Cooperative success in epithelial public goods games
CCooperative success in epithelial public goods games
Jessie Renton ∗ Karen M. Page
Department of Mathematics, University College London,Gower Street, London WC1E 6BT, UK
Abstract
Cancer cells obtain mutations which rely on the production of diffusiblegrowth factors to confer a fitness benefit. These mutations can be consid-ered cooperative, and studied as public goods games within the frameworkof evolutionary game theory. The population structure, benefit function andupdate rule all influence the evolutionary success of cooperators. We modelthe evolution of cooperation in epithelial cells using the Voronoi tessellationmodel. Unlike traditional evolutionary graph theory, this allows us to imple-ment global updating, for which birth and death events are spatially decoupled.We compare, for a sigmoid benefit function, the conditions for cooperation tobe favoured and/or beneficial for well-mixed and structured populations. Wefind that when population structure is combined with global updating, coop-eration is more successful than if there were local updating or the populationwere well-mixed. Interestingly, the qualitative behaviour for the well-mixedpopulation and the Voronoi tessellation model is remarkably similar, but thelatter case requires significantly lower incentives to ensure cooperation.
Keywords— multiplayer games, cooperation, evolutionary game theory, epithelium,Voronoi tessellation, epithelial automata
Oncogenesis is a process of somatic evolution. In order to become cancerous thereare certain key mutations which cells must obtain, corresponding to the hallmarksof cancer [1, 2]. Evolutionary game theory provides a framework for modellingmutations which have a fitness effect beyond the cell itself. For example, certainmutations can be considered cooperative, in that they invoke a cost to the cellwhich is recuperated as a shared benefit. This is evident when the benefit relies onthe production of a diffusible growth factor [3, 4], as is the case for a number ofthe hallmarks of cancer such as self-sufficiency in growth signalling and sustainedangiogesis. The Warburg effect, whereby tumour cells metabolise through glycolysiseven when oxygen is abundant [5], can also be considered cooperative [6]. ∗ Corresponding author: [email protected] a r X i v : . [ q - b i o . P E ] J a n .2 Public goods games Applications of evolutionary game theory to model cancer evolution have mainlyfocussed on two-player games, whereby cells participate in multiple pairwise inter-actions within the population [7–9]. Interactions between cancer cells however, tendto happen in groups. For example, a cell producing a growth factor will providea benefit to other cells within its diffusion range. These types of mutations arethus better represented as multiplayer public goods games (PGGs) [10], played be-tween producer (cooperator) and non-producer (defector) cells. The former producegrowth factor at a fixed cost to their fitness. Both producers and non-producers re-ceive a fitness benefit as a function of the frequency of producers in their interactionneighbourhood.The most common PGG, known as the N-person prisoner’s dilemma (NPD), usesa linear benefit function [11, 12]. However, non-linear benefit functions may be morerealistic, and can lead to much richer dynamics, even for well-mixed populations. Anexample is the volunteer’s dilemma (VD), which defines the benefit as a Heavisidestep function [13–16].A sigmoid benefit function has been proposed as an appropriate model for growthfactor production. Experiments on neuroendocrine pancreatic cancer cells in vitro have found sigmoid dependence of proliferation rates on the concentration of growthfactor IGF-II [17]. Furthermore such a function is relatively general, with both theNPD and VD arising as special cases [18].
Most cancers originate in epithelia. These are tissues formed of sheets of cells, whichare approximately polygonal on their apical surfaces. It is important to take intoaccount this population structure when modelling the evolutionary dynamics. Forboth two-player cooperation games [19, 20] and multiplayer PGGs [21], cooperatorstend to have greater success in structured populations, as compared to well-mixedones, because they are able to form mutually beneficial clusters. Epithelial cellstend to have six neighbours on average, and thus can be represented as a hexagonallattice. Introducing more realistic population structures, with small variation inneighbour number, does not have a significant impact on evolutionary outcomes[22, 23].The success of cooperation is also dependent on the update dynamics. Evolutionon structured populations is usually modelled within the framework of evolutionarygraph theory [24], in which the population is represented as a fixed graph. Thepopulation then evolves according to an update rule. In general update rules canbe divided into two categories: local and global [25].
A local update involves a spatial relationship between birth and death events. Evo-lutionary graph theory usually requires a local update rule in order to maintain thefixed graph structure. Two commonly used local update rules are defined as follows:• birth-death: a cell is selected to divide with probability proportional to fitness;one of its neighbours is chosen to die uniformly at random.• death-birth: a cell is chosen to die uniformly at random; one of its neighboursis selected to divide with probability proportional to fitness.2n both cases the offspring of the dividing cell occupies the empty site left by thedead cell [26]. The choice between these update rules has a substantive effect onevolutionary outcomes. For example, consider a two-player prisoner’s dilemma gameand a population represented by a regular graph. Cooperation can be favoured fora death-birth update rule, so long as the benefit is high enough. For the birth-deathupdate however, as is the case with a well-mixed population, cooperation is onlyfavoured for an infinitely large benefit [19]. Other local update rules, such as thosefor which selection acts on death, can be found in the literature [27].
Under a global update rule there is no spatial dependence between birth and deathevents, thus cells are selected to reproduce and die from the population as a whole.Global updating is generally seen for well-mixed populations, or when populationsare organised in phenotype space [28] or by sets [29].Within evolutionary graph theory the shift update rule is an example of globalupdating. In this case a cell is chosen to divide with probability proportional tofitness, and a second cell is chosen to die. A path is then selected on the graphwhich connects the two. Cells are shifted along this path until there is an emptynode next to the dividing cell for its progeny to occupy. This kind of update workswell on a one-dimensional lattice [30], and promotes cooperation, even compared tothe death-birth update. However, it becomes more complex in two-dimensions [31].Alternatively, a previous study has used the Voronoi tessellation model [32] toincorporate global updating consistent with epithelial structure [23]. Unlike tradi-tional evolutionary graph theory, the structure is dynamic, thus cells are able todivide and die independently. Using a two-player prisoner’s dilemma game, cooper-ation was found to be more successful in this context, than for a death-birth updaterule.
For stochastic evolutionary games without mutation, we can compare the successof different strategies by calculating fixation probabilities. Here we consider thedynamics of two cell types: A and B . The fixation probability ρ X is then definedas the probability that a single initial mutant X will eventually take over the entirepopulation. We consider two measures for the success of an A mutant [26, 33]:• A is a beneficial mutation when ρ A > ρ . Here ρ = 1 /Z is the fixationprobability for a neutral mutant and Z is the population size.• A is favoured by selection, or has an evolutionary advantage, when ρ A > ρ B .This is equivalent to the condition that the equilibrium frequency of A isgreater than a half when mutation is allowed ( A is the dominant strategy).In general, these conditions are not equivalent, thus it is possible for a mutationto be beneficial but not favoured, or vice versa. One or the other condition mightbe more relevant to quantifying mutant success depending on the circumstances.Furthermore, under certain conditions these two conditions are equivalent [33].The remainder of this paper explores conditions under which a mutation is ben-eficial and/or favoured. We begin in Section 2 by setting out the mathematical for-malism for multiplayer evolutionary games, focussing particularly on PGGs playedbetween cooperators and defectors. Section 3 then introduces the σ -rule, which is3sed to determine whether a strategy is favoured. We outline several known re-sults on graphs with local update rules, as well as deriving results for a birth-deathand shift update rule on a cycle. We then derive the conditions for favourabilityon a general population structure with global updating. In Section 4 we derive asimilar rule, but for a strategy to be beneficial. In Section 5 we apply this theoryto consider conditions for cooperator success in an epithelium, using spatial statis-tics calculated through simulation of the Voronoi tessellation model. Finally, inSection 6, we discuss the implications of our work for the evolution of cooperativepublic goods in epithelia and make some remarks on the different significance ofbeneficial and favourable mutants. We consider an arbitrary multiplayer game with two strategies, A and B . Playersinteract in groups of size N = k + 1 , and obtain payoffs a j,k and b j,k respectively,where j is the number of A co-players and k is the total number of co-players. Fora graph-structured population, the co-players are direct neighbours. The fitness ofeach individual is then defined as δa j,k or δb j,k , where δ is the selectionstrength parameter.The population evolves according to a Moran process [34], i.e. at each time-stepone individual dies and another reproduces, thus keeping the population size, Z ,constant. How these individuals are chosen is determined by the update rule. Weconsider cases where reproduction, but not death, is dependent on fitness.Many of the results we derive in the following sections are for general games,however we are focussed on PGGs played between producer/cooperator cells ( C )and non-producer/defector cells ( D ). These games are defined by a benefit function b · β ( x ) and a cost function which we take to be constant c , with b > c . Here x is the proportion of cooperators in a cell’s interaction group. Thus cooperator anddefector payoffs are defined respectively as a j,k = b · β (cid:18) j + 1 k + 1 (cid:19) − c , b j,k = b · β (cid:18) jk + 1 (cid:19) . (1)In order to ensure that the payoff is higher when all players cooperate than whenno players cooperate we enforce the condition b · β ( N ) − c ≥ b · β (0) . Often this isdone by setting c = 1 , b > , β ( N ) = 1 and β (0) = 0 .The NPD and VD can both be defined in this form by specifying the benefitfunctions: β ( x ) = x (NPD) (2) β ( x ) = Θ( x − ˜ x ) , (VD) (3)where Θ( x ) is the Heaviside step function and ˜ x is the minimum proportion ofcooperators required to obtain the benefit. Furthermore we can define a generalsigmoid benefit function: β ( x ) = α ( x ) − α (0) α (1) − α (0) , (4)where α ( x ) = 11 + e s ( h − x ) (5)is the logistic function, s is the steepness and h is the inflection point. We can regainthe NPD and VD by taking the limits s → and s → ∞ respectively (see Figure 1).4 x ( x ) x Figure 1: Logistic benefit function.
Left panel: h = 0 . ; s = 1 (dash-dot), s = 10 (dash)and s = 1000 (solid). We can regain the limiting cases by letting s → (NPD) or s → ∞ (VD). Right panel: s = 10 ; h = 0 . (dash-dot), h = 0 . (dash), h = 0 . (solid). σ -rule: conditions for cooperation to be favoured For a particular update rule and population structure, the σ -rule allows us to deter-mine which is the favoured strategy [35]. We recall from Section 1.4, that a strategy A is favoured over B , when ρ A > ρ B .The σ -rule states that ρ A > ρ B ⇐⇒ k (cid:88) j =0 σ j ( a j − b k − j ) > , (6)where σ j are the structure coefficients. It is assumed that the group size, N = k + 1 , is fixed, thus we have let a j,k = a j and b j,k = b j . The structure coefficientsare dependent on the population structure and update rule, but not the payoffs.Therefore if we calculate σ j for a given population structure and update rule, wecan determine the favoured strategy for any game.For certain population structures, such as the well-mixed population and thecycle graph, the state is fully described by the number of A -players, n . Thus we candefine the ratio of fixation probabilities as ρ A ρ B = Z − (cid:89) n =1 T + n T − n , (7)where T ± n are the transition probabilities to go from n → n ± A -type individuals[36]. This does not hold in general, as the transition probabilities in more complexpopulation structures will depend on the spatial configurations of different cell types,and thus are not uniquely defined by n . However, it is still possible to utilisethis equation, as we see in Sections 3.4 and 4, by averaging over possible states toapproximate T ± n .In the following we consider various cases where the structure coefficients can becalculated from transition probabilities in the weak selection limit, i.e. when δ (cid:28) . The structure coefficients for a well-mixed population are given by [37] σ j = (cid:40) , if ≤ j ≤ N − Z − NZ , if j = N − (8)5see also Section 3.4). Thus we can obtain the condition for ρ A > ρ B by pluggingthese into Equation (6). For a PGG defined by Equation (1) the condition thatcooperators are favoured is Z − NZ b [ β (1) − β (0)] > N − (cid:88) j =0 σ j c . (9)This becomes bc > N ( Z − Z − N , (10)when we set β (1) = 1 and β (0) = 0 . Clearly the shape of the benefit function doesnot impact whether cooperation is favoured. For a large population this conditionbecomes b/c > N . We can obtain exact expressions for the structure coefficients of the cycle graph,in the weak selection limit. The cycle is a one-dimensional lattice with periodicboundary conditions. Individuals interact with their two nearest-neighbours, thuswe have group size N = 3 . The structure coefficients for the death-birth update rule are derived in [21]. Theyare given by σ = 1 , σ = Z − , σ = Z − . (11)From Equation (6) we obtain the condition for cooperation to be favoured under anNPD, defined by Equation (2): bc > Z − Z − , (12)which for Z → ∞ becomes b/c > / . These conditions are lower than thoseobtained for a well-mixed population. For a general PGG defined by Equation (1)we can write down the condition bc > Z − Z − β (1) + β (2 / − β (1 / − β (0)] . (13) We calculate the structure coefficients for the birth-update rule, following the methodused in [21]. For the cycle, the transition probabilities are uniquely defined by thenumber of A -players in the population, n . Thus we can write down the ratio oftransition probabilities for each n . For a birth-death update rule these are T + n T − n = (1 + δa ) / (1 + δb ) , if n = 1(1 + δa ) / (1 + δb ) , if < n < Z − δa ) / (1 + δb ) , if n = Z − . (14)Substituting these into Equation (7), and taking the limit, δ (cid:28) , we obtain ρ A ρ B ≈ δ [ a − b + ( Z − a − b )] . (15)6n order that ρ A > ρ B , the second term must be positive. Thus, comparing thiscondition with Equation (6), we find the structure coefficients σ = 1 , σ = Z − , σ = 0 . (16)For the NPD, cooperation is favoured when bc > Z − Z − , (17)which becomes b/c > in the large population limit, Z → ∞ . These conditionsare equivalent to those obtained for the well-mixed population. For a general PGGdefined by Equation (1) the condition is bc > Z − Z − β (2 / − β (1 / . (18) Again we follow the same procedure to derive the structure coefficients for the shiftupdate rule. This time the ratio of transition probabilities is given by T + n T − n = ( Z − δa )2(1+ δb )+( Z − δb ) , if n = 1 ( Z − n )(2(1+ δa )+( n − δa )) n (2(1+ δb )+( Z − n − δb )) , if < n < Z − δa )+( Z − δa )( Z − δb ) , if n = Z − . (19)In the weak selection limit, δ (cid:28) , Equation (7) becomes ρ A ρ B ≈ δ [( a − b ) + 2( H Z − − a − b )) + ( Z − H Z − )( a − b )] , (20)where H m is the m -th harmonic number: H m = m (cid:88) n =1 n . (21)Thus the structure coefficients are given by σ = 1 , σ = 2( H Z − − , σ = ( Z − H Z − ) . (22)The condition for cooperation to be favoured in the NPD is bc > Z − Z − − H Z − . (23)In the large population limit this becomes b/c > . As this condition is required inthe definition of the NPD, we can state that cooperation is always favoured in thelarge population limit for a shift update under weak selection.In fact, if we consider a general cooperation game as defined by Equation (1) weobtain the condition bc > β (1) − β (0) (24)in the large population limit, Z → ∞ . Letting β (1) = 1 and β (0) = 0 , we regain thecondition b/c > . Thus for the shift update on the cycle, as with the well-mixedpopulation, the condition for cooperation to be favoured is not dependent on theshape of the benefit function (although in this case we required the large populationlimit). Furthermore cooperation is favoured on the cycle with shift update for allPGGs, as defined by Equation (1), given that the population is sufficiently large.7 .3 Approximate results for k -regular graphs In higher dimensions the transition probabilities are no longer uniquely defined bythe number of A -players in the population, but depend also on their configuration.Peña et al [21] have derived expressions for the structure coefficients of regular graphsof degree k ≥ , with death-birth updating, using pair approximation and diffusionapproximation [19]. They compared theoretical predictions with simulation resultsfor the case of a volunteer’s dilemma game. They find a good fit for random regulargraphs, but that the approximations underestimate the critical benefit-to-cost ratiofor lattices.We do not state the full expressions here which are non-trivial functions of k .The condition for cooperation to be favoured with the NPD in the large populationlimit ( Z (cid:29) k ) is given by bc > k + 12 . (25) In the following we derive a general expression for the structure coefficients undera global update rule. The proceeding sections have considered games played on afixed graph or well-mixed population, within groups of fixed size, N . For well-mixedpopulations we were free to choose N (although some results required N (cid:28) Z ),while for regular graphs we set N = k + 1 , where k is the degree of the graph. Herewe relax this condition and allow for variable group size.We make the assumption that there is a fixed distribution, f A/Bj ( n, k ) , definingthe probability that an A/B -player interacts with j co-players of type A , given ithas k co-players in total and there are n players of type A in the population. Ifthe population were defined on a graph, this would be the probability of an A/B -player having j A -type neighbours, given k total neighbours. This assumption istrue for a well-mixed population or cycle graph, but not necessarily for other popu-lation structures where f A/Bj ( n, k ) depends on the specific configuration of players.The frequency of individuals with k neighbours is given by g k . We make the fur-ther assumptions that this distribution is fixed, and does not depend on type. SeeAppendix B for a discussion of the validity of this assumption for the VT model.In general, for a global update rule, we can define the transition probabilities T + n = Z − nZ nF A nF A + ( Z − n ) F B T − n = nZ ( Z − n ) F B nF A + ( Z − n ) F B , (26)where F A = 1 + δ Z − (cid:88) k =1 k (cid:88) j =0 f Aj ( n, k ) g k a j,k (27) F B = 1 + δ Z − (cid:88) k =1 k (cid:88) j =0 f Bj ( n, k ) g k b j,k (28)are the population averaged fitnesses. The payoffs a j,k and b j,k depend explicitly onthe number of neighbours k .Substituting Equations (26) to (28) into Equation (7), and taking the weakselection limit we obtain ρ A ρ B ≈ δ Z − (cid:88) n =1 Z − (cid:88) k =1 k (cid:88) j =0 g k [ f Aj ( n, k ) a j,k − f Bj ( n, k ) b j,k ] (cid:124) (cid:123)(cid:122) (cid:125) Γ . (29)8hus ρ A > ρ B when Γ > . In the weak selection limit, f Aj ( n, k ) = f Bk − j ( Z − n, k ) (30)must hold by symmetry, and thus Z − (cid:88) n =1 f Aj ( n, k ) = Z − (cid:88) n =1 f Bk − j ( n, k ) . (31)Therefore we have Γ = Z − (cid:88) k =1 k (cid:88) j =0 Z − (cid:88) n =1 g k f Aj ( n, k )( a j,k − b k − j,k ) . (32)The condition for A to be favoured over B is thus given by Z − (cid:88) k =1 k (cid:88) j =0 σ j,k ( a j,k − b k − j,k ) > , (33)where σ j,k = g k Z − (cid:88) n =1 f Aj ( n, k ) (34)are the structure coefficients. For a fixed group size, N = k + 1 , this reduces toEquation (6), with σ j = Z − (cid:88) n =1 f Aj ( n ) , (35)where we have dropped the explicit dependence on k .Recall that this derivation is based on the assumption that g k and f Aj ( n, k ) arefixed. While this is not true in most cases, we can obtain an approximation for thestructure coefficients by averaging over a large ensemble of population configurations,i.e. letting f Aj ( n ) = (cid:104) f Aj ( n ) (cid:105) . Here, (cid:104) . (cid:105) represents the mean taken over possibleconfigurations and the indicates that these are obtained in the neutral selectionlimit, i.e. δ = 0 .The well-mixed population is an example where f Aj ( n ) is fixed. It is defined bya hypergeometric distribution: f Aj ( n ) = (cid:18) Z − k (cid:19) − (cid:18) n − j (cid:19)(cid:18) Z − nk − j (cid:19) . (36)We can therefore find the structure coefficients [37] by substituting this expressionfor f Aj ( n ) into Equation (35): σ j = (cid:18) Z − k (cid:19) − Z − (cid:88) n =1 (cid:18) n − j (cid:19)(cid:18) Z − nk − j (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) S . (37)The sum can be written S = Z − (cid:88) m =0 (cid:18) mj (cid:19)(cid:18) Z − m − k − j (cid:19) − (cid:18) Z − j (cid:19)(cid:18) k − j (cid:19) . (38)9or ≤ j < k the second term vanishes and we obtain S = (cid:0) Zk +1 (cid:1) , using theChu–Vandermonde identity. For j = k we have S = (cid:0) Zk +1 (cid:1) − (cid:0) Z − j (cid:1) . Thus thestructure coefficients are given by σ j = (cid:40) Zk +1 if ≤ j < k Z − k − k +1 if j = k . (39)These are equivalent to Equation (8) up to a constant factor. The cycle graph alsohas a fixed distribution, f Aj ( n ) , thus the structure coefficients for the shift updaterule can also be obtained exactly using Equation (35).For a variable group size the structure coefficients for the well-mixed populationare given by σ j,k = g k σ j ( k ) = g k (cid:40) Zk +1 if ≤ j < k Z − k − k +1 if j = k, (40)where σ j ( k ) are defined in Equation (39).As we have seen in previous sections, once the structure coefficients have beendetermined, we can use Equation (6) or Equation (33) to find the condition underwhich cooperation is favoured. For a PGG as defined by Equation (1) this is givenby bc > Z − (cid:80) Z − k =1 (cid:80) kj =0 σ j,k (cid:2) β (cid:0) j +1 k +1 (cid:1) − β (cid:0) k − jk +1 (cid:1)(cid:3) . (41) Thus far we have considered conditions where the fixation probability for a single A -mutant is greater than the fixation probability for a single B -mutant, i.e. ρ A >ρ B . Given this condition holds, we can say that A is favoured by selection, or hasan evolutionary advantage over B . It has also been shown previously that wheremutation is allowed, this condition corresponds to the point where A is more frequentin the population than B [38].An alternative measure of a mutant’s success can be obtained by comparing it toa neutral mutant, i.e. considering the condition ρ A > ρ . If this condition is satisfiedwe say the mutation is beneficial. These conditions do not always coincide, thus itis possible for cooperation to be beneficial but not favoured or vice versa. In thefollowing we derive the condition for an A -mutant to be beneficial.As in the previous section, we make the assumption that the distributions g k and f A/Bj ( n, k ) are fixed. Thus the population averaged fitnesses of A and B players aredefined by Equations (27) and (28) and the transition probabilities by Equation (26).The fixation probability for a single A -mutant [36] is then given by ρ A = (cid:34) Z − (cid:88) m =1 m (cid:89) n =1 T − n T + n (cid:35) − . (42)Substituting in the transition probabilities and taking the weak selection limit δ (cid:28) we obtain ρ A = 1 Z + δZ Z − (cid:88) k =1 k (cid:88) j =0 (cid:0) θ Aj,k a j,k − θ Bj,k b j,k (cid:1) + O ( δ ) , (43)10here we have defined θ Aj,k = g k Z − (cid:88) m =1 m (cid:88) n =1 f Aj ( n, k ) (44) θ Bj,k = g k Z − (cid:88) m =1 m (cid:88) n =1 f Bj ( n, k ) = g k Z − (cid:88) m =1 m (cid:88) n =1 f Ak − j ( Z − n, k ) . (45)The final equality is obtained by symmetry arguments in the weak selection limit.The condition for A to be a beneficial mutation, ρ A > /Z , is therefore given by Z − (cid:88) k =1 k (cid:88) j =0 (cid:0) θ Aj,k a j,k − θ Bj,k b j,k (cid:1) > . (46)If we consider a PGG as defined by Equation (1), then cooperation is beneficialwhen bc > Z ( Z − (cid:80) Z − k =1 (cid:80) kj =0 (cid:2) θ Aj,k β (cid:0) j +1 k +1 (cid:1) − θ Bj,k β (cid:0) jk +1 (cid:1)(cid:3) . (47)For a fixed group size N = k + 1 these conditions simplify to k (cid:88) j =0 (cid:0) θ Aj a j − θ Bj b j (cid:1) > (48)and bc > Z ( Z − (cid:80) kj =0 (cid:2) θ Aj β (cid:0) j +1 k +1 (cid:1) − θ Bj β (cid:0) jk +1 (cid:1)(cid:3) , (49)where θ A/Bj = Z − (cid:88) m =1 m (cid:88) n =1 f A/Bj ( n, k ) . (50) A number of studies have considered the evolutionary dynamics of sigmoid PGGsin epithelia, representing the tissue either as a well-mixed population [39], or a fixedgraph structure with various local update rules [22, 40]. Here we use the frameworkintroduced in [23] to incorporate explicit tissue dynamics, using the Voronoi tes-sellation (VT) model, with a spatially decoupled (global) update rule. This meansthat when the population is updated, a division and death occur simultaneously,but there is no spatial dependence between the two events.In this section we briefly introduce the VT model, before calculating conditionsunder which cooperation is favoured and beneficial for a sigmoid PGG. We verifytheoretical results by running simulations in various parameter regimes. We alsocompute the gradient of selection in order to obtain a fuller picture of the dynamics.In all cases we compare VT model results with the well-mixed population.
The VT model represents a tissue as a set of points, corresponding to cell centres[32, 41]. The shape of each cell, as well as its neighbour connections, is determinedby performing a Voronoi tessellation. Cells move subject to spring-like forces which11hey exert on their neighbours. A full description of the VT model used is given inAppendix A.We obtain g k and f Aj ( n, k ) by averaging over a large ensemble of possible statesin the weak selection limit. We then make the assumption that variation around thismean can be neglected. Figures 2 and 3 show the distributions g k and f Aj ( n, k ) for theVT model under neutral selection, calculated by averaging over simulations, eachof which starts with a single neutral mutant and is run to fixation. See Appendix Bfor further discussion on neighbour distributions in the VT model and the validityof assuming g k is independent of n and cell type. k F r e q u e n c y , g k Figure 2: Degree distribution for the Voronoi tessellation model. Error bars show standarddeviation. Data is obtained from simulations with population size, Z = 100 . The condition for cooperation to be favoured can be approximated by calculatingthe structure coefficients using Equation (34). Figure 4 plots the VT structurecoefficients with those for a well-mixed population as defined by Equation (39).Using the structure coefficients we can derive the condition for cooperation to befavoured for an arbitrary PGG, as defined by Equation (1). We define the criticalbenefit-to-cost ratio ( b/c ) ∗ , such that ρ C > ρ D when b/c > ( b/c ) ∗ . Thus fromEquation (41) we can write (cid:18) bc (cid:19) ∗ = Z − (cid:80) Z − k =1 (cid:80) kj =0 σ j,k (cid:2) β (cid:0) j +1 k +1 (cid:1) − β (cid:0) k − jk +1 (cid:1)(cid:3) . (51)For an NPD, defined by Equation (2), this becomes (cid:18) bc (cid:19) ∗ = Z − (cid:80) Z − k =1 (cid:80) kj =0 σ j,k j +1 − kk +1 . (52)Substituting in the structure coefficients we obtain ( b/c ) ∗ ≈ . for the VTmodel with decoupled update rule and population size Z = 100 . For a well-mixedpopulation with the same group size distribution we obtain ( b/c ) ∗ ≈ . . As wewould expect there is a significant increase in the success of cooperative mutantsunder the VT model. This is due to the high level of assortment in the VT model,which means cooperators are likely to have more cooperator neighbours than defec-tors. 12 f A j ( n , k ) k = 5 0 100nk = 6 0 100nk = 70 100n0.00.20.40.60.81.0 f B j ( n , k ) k = 5 0 100nk = 6 0 100nk = 7j=0 7 Figure 3: Frequency distributions f Aj ( n, k ) and f Bj ( n, k ) for Z = 100 . These define theprobability that a cell of type A/B has j neighbours of type A , given k neighbours totaland n cells of type A in the population. The lower panel compares values of f Bj ( n, k ) calculated directly through simulation (dashed) with values obtained from the simulateddata for A cells defined by f Bj ( n, k ) = f Ak − j ( Z − n, k ) . j , k k = 5 0 5jk = 6 0 5jk = 7 VTWM Figure 4: Comparing the structure coefficients for the Voronoi tessellation model withdecoupled update (VT) and a well-mixed (WM) population. Both have the same neighbourdistribution seen in Figure 2.
13n average, cells have six neighbours, thus the mean group size is seven. Wecan therefore compare the critical benefit-to-cost ratio for a well-mixed populationwith variable group size, given above, to that of a well-mixed population with fixedgroup size, N = 7 . The latter is given by Equation (10) to be ( b/c ) ∗ = 7 . . Clearly,incorporating variation in group size has a negligible impact on whether cooperationis favoured in this case. ( b / c ) * s = 1 s = 50 1h2.55.07.5 ( b / c ) * s = 10 0 1hs = 20 VTHLWM Figure 5: Comparing the critical benefit-to-cost ratio, ( b/c ) ∗ at which ρ C > ρ D , for alogistic benefit function. For a well-mixed population with N = 7 (WM), ( b/c ) ∗ is highest,and independent of the inflection point, h , and steepness, s . For the Voronoi tessellationmodel with decoupled update (VT) and fixed hexagonal lattice with death-birth update(HL), ( b/c ) ∗ varies with h and s . For small s the benefit function approaches linearity andwe regain the results for an NPD. We can also use Equation (51) to determine ( b/c ) ∗ for a sigmoid benefit function,defined by Equation (4). Recall that the logistic function has two parameters: thesteepness, s , and the inflection point, h . Figure 5 compares the predicted valuesof ( b/c ) ∗ for the VT model, with those for a well-mixed (WM) population withgroup size 7, and hexagonal lattice (HL) with death-birth update rule. These areobtained from Equation (51) by using the relevant structure coefficients in each case(structure coefficients for death-birth update on regular graphs are derived in [21]).Values of ( b/c ) ∗ are symmetric across h = 0 . for all three cases, and minimisedat h = 0 . for the hexagonal lattice and VT model. In Appendix C we show that ( b/c ) ∗ is in fact minimised at h = 0 . , so long as the structure coefficients increasewith j for ≤ j < k . For the well-mixed population ( b/c ) ∗ does not vary witheither s or h . Furthermore, it is clear for all population types, that as the NPD isapproached ( s → ), ( b/c ) ∗ becomes independent of h .In all parameter regimes, ( b/c ) ∗ is highest for the well-mixed population. Boththe VT model with decoupled update and HL with death-birth update show similarvariation with s and h , however ( b/c ) ∗ is always lower for the VT model. There-fore in terms of thresholds for favourability, we can determine that cooperation ismost successful in the VT model with decoupled update, followed by the hexagonallattice with death-birth update. Cooperation does least well in the well-mixed pop-ulation. This suggests that both population structure and global updating promotecooperation. 14igure 6 (right panel) shows the variation of ( b/c ) ∗ with h and s for the VTmodel. As we have discussed, these results are based on the approximation that f Aj ( n, k ) and g k are fixed. In order to verify the accuracy of this approximationwe compare Equation (51) with simulation results in Figure 7. Simulated values of ( b/c ) ∗ were obtained for each parameter set ( s, h ) as follows. We calculated ρ C/D forvarious b/c values, by running simulations of the VT model to fixation, startingwith a single C/D mutant and population size Z = 100 . In all simulations we usesmall selection strength ( δ = 0 . ) and set c = 1 . Thus ( b/c ) ∗ is determined by thepoint at which ρ C = ρ D . s ( b / c ) *0 ( b / c ) *1 Figure 6: Critical benefit-to-cost ratios for the VT model with decoupled update. Theseare given by Equations (51) and (53) for a PGG with logistic benefit function, definedby Equation (4). Parameters s and h correspond to the steepness and inflection point ofthe benefit function, respectively. Cooperation is beneficial when b/c > ( b/c ) ∗ (left) andfavoured when b/c > ( b/c ) ∗ (right). ( b / c ) * s = 5 0 1hs = 10 Figure 7: Critical benefit-to-cost ratio, ( b/c ) ∗ , above which ρ C > ρ D , for a logistic benefitfunction. The solid line plots Equation (51) and circles are simulation data. For both s = 5 and s = 10 there is symmetry across h = 0 . , at which point ( b/c ) ∗ is minimised. There is a decent fit between simulation and theory. It is possible this could beimproved by running larger numbers of simulations, however the model is computa-tionally expensive. In any case the qualitative behaviour is consistent. For a fixedsteepness, s , ( b/c ) ∗ is minimised at h = 0 . and (near) symmetric across this value.The values of ( b/c ) ∗ are highest when h = 0 and h = 1 , where the benefit functionprovides diminishing returns or increasing returns respectively.15 .3 Beneficial cooperation Thus far we have considered conditions for cooperation to be favoured, i.e. where ρ C > ρ D . We can also define the critical benefit-to-cost ratio ( b/c ) ∗ above whichcooperation is beneficial, i.e. ρ C > ρ . From Equation (47) this is given by (cid:18) bc (cid:19) ∗ = Z ( Z − (cid:80) Z − k =1 (cid:80) kj =0 (cid:2) θ Aj,k β (cid:0) j +1 k +1 (cid:1) − θ Bj,k β (cid:0) jk +1 (cid:1)(cid:3) , (53)where θ A/Bj,k are calculated from the distributions f A/Bj ( n, k ) and g k according toEquations (44) and (45).Figure 6 (left panel) plots ( b/c ) ∗ against s and h . We can see that for large s , ( b/c ) ∗ is maximised at h = 1 and has a minimum at h ≈ . . For smaller s this minimum moves towards h = 0 . As s decreases further, the logistic functionapproaches linear and there is negligible variation in ( b/c ) ∗ with h . In the limit s → the game becomes an NPD, with ( b/c ) ∗ = ( b/c ) ∗ ≈ . . Figure 8 comparessimulated values of ( b/c ) ∗ with the theoretical prediction, finding good agreementbetween the two for a range of s and h values. ( b / c ) * s = 1 0 1hs = 5 0 1hs = 10 Figure 8: Critical benefit-to-cost ratio, ( b/c ) ∗ , above which ρ C > ρ , for a logistic benefitfunction. The solid line plots Equation (53) and circles are simulation data. For small s the logistic benefit function becomes near linear and the game approaches an N P D , thusthere is little variation in ( b/c ) ∗ . For larger s there is strong dependence on the inflectionpoint, h , particularly for h > . . We saw in Figure 5 that the critical benefit-to-cost ratios for cooperation to befavoured, ( b/c ) ∗ , are lower in the VT model compared to the well-mixed population.Figure 9 plots ( b/c ) ∗ for a well-mixed population with N = 7 and the VT modelwith decoupled update, showing clearly that the critical benefit-to-cost ratios forcooperation to be beneficial are also lower for the VT model. Thus under bothmeasures, cooperation is promoted by the VT model. In contrast to ( b/c ) ∗ , whichwas independent of the shape of the benefit function for the well-mixed population, ( b/c ) ∗ is an increasing function of h , so long as s is sufficiently large.In general, conditions for cooperation to be beneficial are not equivalent to con-ditions for cooperation to be favoured. This is evident from Figure 10, where weplot ( b/c ) ∗ and ( b/c ) ∗ against h for various values of s . The parameter space canbe divided into regions where cooperation is both favoured and beneficial, favouredbut not beneficial, beneficial but not favoured, and neither favoured nor beneficial.From Figure 10 we can see that ( b/c ) ∗ = ( b/c ) ∗ when h = 0 . . Furthermore,as s decreases, the regions of parameter space where cooperation is beneficial, butnot favoured, or favoured, but not beneficial, get smaller. For sufficiently small s we obtain ( b/c ) ∗ ≈ ( b/c ) ∗ . For both the VT model and well-mixed populations16 .0 0.5 1.0h5101520 ( b / c ) * s = 1 0.0 0.5 1.0hs = 5 0.0 0.5 1.0hs = 10 VTWM Figure 9: Comparing the critical benefit-to-cost ratio, ( b/c ) ∗ at which ρ C > ρ , for aPGG with logistic benefit function. The critical ratio is always higher for the well-mixedpopulation with N = 7 (WM), than for the Voronoi tessellation model with decoupledupdate (VT). For small s the benefit function becomes near linear and variation of ( b/c ) ∗ with h is small. For WM, ( b/c ) ∗ increases with h , taking its minimum value at h = 0 . Bycontrast, for VT, there is a minimum of ( b/c ) ∗ at h ≈ . when s is sufficiently large. Forboth WM and VT, ( b/c ) ∗ is maximised at h = 1 , for any given s . it is clear that behaviour where cooperation is beneficial but not favoured, is onlypossible when h < . . Conversely behaviour where cooperation is favoured but notbeneficial occurs only when h > . .We can understand this intuitively by considering the extreme cases ( h = 0 , )of the VD game, obtained by letting s → ∞ . When h = 0 , a cooperator alwaysreceives the full benefit, even if it has no cooperator neighbours. Defectors requirea single cooperator neighbour to obtain the benefit. Thus both cooperators anddefectors have higher than average fitness early on in the invasion process, whenthey are most vulnerable to extinction. It is therefore possible, depending on thebenefit-to-cost ratio, that both perform better than a neutral invader, and thereforeboth are beneficial mutations. However, one can still be favoured over the other ifits fixation probability is higher.The converse is true when h = 1 : defectors will never receive any benefit, andcooperators only obtain the benefit when surrounded by other cooperators. Thuswhen the number of cooperators/defectors is small, they have lower than average fit-ness, and there is a high chance they die out early in the invasion process. Therefore,it is possible that neither performs better than a neutral invader. We can obtain more insight into what is happening in the different parameter regionsby looking at the gradient of selection, G ( n ) = T + ( n ) − T − ( n ) . The transitionprobabilities are defined by Equation (26), thus in the weak selection limit, δ (cid:28) ,the gradient of selection becomes G ( n ) ≈ Z − nZ nZ δ (cid:40) Z − (cid:88) k =1 k (cid:88) j =0 g k ( f Aj ( n, k ) a j,k − f Bj ( n, k ) b j,k ) (cid:41) . (54)The sum essentially gives the difference in expected payoffs of A and B players.Thus, the right-hand side is identical to the replicator equation which describes thedeterministic dynamics which can be obtained in the large-population limit. Thesign of G ( n ) indicates the direction of selection, and we can consider the roots of G ( n ) as ‘fixed points’. Of course, for a finite population there are only two absorbing17 b / c VT s = 1 s = 50 1h246 b / c s = 10 0 1hs = 20 61626 b / c WM s = 1 s = 50 1h61626 b / c s = 10 0 1hs = 20 Figure 10: Success of cooperator mutants in the VT model (left) and well-mixed population(right), for a PGG with logistic benefit function. The solid line corresponds to ( b/c ) ∗ , where ρ C = ρ . The dashed line corresponds to ( b/c ) ∗ , where ρ C = ρ D . Blue region (top): C is beneficial and favoured ( ρ C > ρ D and ρ C > ρ ). Green region (left): C is beneficialbut not favoured ( ρ D > ρ C > ρ ). Pink region (right): C is favoured but not beneficial( ρ > ρ C > ρ D ). Orange region (bottom): C is neither beneficial not favoured ( ρ C < ρ D and ρ C < ρ ). states, n = 0 and n = Z , however the location of fixed points is still important. Inparticular the system may remain for a long time near a stable fixed point. Wecan classify the behaviour of the system in different parameter regions based on thefixed points of the gradient of selection.Figure 11 plots G ( n ) for a PGG with various values of h , s and b/c , both for theVT model and well-mixed population. There are four dynamical regimes, consistentwith the deterministic results for PGG in a well-mixed population in [39]:(i) Dominance: there are only two fixed points at n = 0 and n = Z . Defectiondominates if the n = 0 fixed point is stable, while cooperation dominates ifthe n = Z fixed point is stable.(ii) Coexistence: there is an internal stable fixed point, n R , along with two unstablefixed points at n = 0 and n = Z . Selection pushes the system towards thestable fixed point, thus it can take a long time to reach one of the absorbingstates.(iii) Coordination: there is an internal unstable fixed point, n L , along with twostable fixed points at n = 0 and n = Z .(iv) Coexistence & coordination:
In addition to the fixed points at n = 0 (stable)and n = Z (unstable), there is both an unstable internal fixed point on theleft, n L , and a stable internal fixed point on the right, n R . Thus it resemblescoexistence, in that there is a stable mixed state; and coordination in thatthere are two stable fixed points.These regimes are all familiar in the evolutionary game theory literature for well-mixed populations. The first three correspond to the behaviour of two-player ma-trix games in well-mixed populations [36]: (i) prisoner’s dilemma/harmony game,18ii) snowdrift game, and (iii) stag-hunt game. The final type, coexistence & coor-dination, arises for both the N-player stag-hunt [42] and N-player snowdrift games[43].For the well-mixed population we see dominance when s is sufficiently small, andthus the PGG is approximating an NPD. As expected, cooperation is dominant when b/c is sufficiently high. For higher values of s there is a wide range of behaviour. Ina region around h = 0 . , if b/c is large enough, there are coexistence & coordinationdynamics. There is a large basin of attraction for n R and if the system reaches thisfixed point it will spend a long time in the vicinity. However, a single mutant invadermust cross n L to reach this, against the selection pressure. As b/c is increased, n L and n R move further apart, increasing the size of the basin of attraction for n R . For h = 0 . , the gradient of selection is symmetric ( n L = Z − n R ).Decreasing h from . , causes n L and n R to move to the left, eventually enteringthe coexistence regime. The basin of attraction for the internal stable fixed pointis now < n < Z . The system may spend a large amount of time near this point,although it will ultimately end up in one of the absorbing states. In the coexistenceregime, as we discussed in Section 5.3 for the VD game with h = 0 , cooperators anddefectors have a selective advantage when they are in sufficiently small numbers.This can lead to the case where both are beneficial mutants, and thus cooperationis able to be beneficial but not favoured.Conversely as h is increased from . , n L and n R move to the right and we enterthe coordination regime. This corresponds to the region in Figure 10 where veryhigh values of the benefit-to-cost ratio are required for cooperation to be beneficial,even when cooperation is favoured. In Section 5.3 we argued, for the VD game with h = 0 . , that this is due to the fact that both cooperators and defectors are ata disadvantage when in small numbers. Indeed this is the defining feature of thecoordination regime, that n = 0 and n = 1 are stable fixed points. Thus any invaderis at a disadvantage initially, as selection pushes it towards dying out. Therefore itis possible that defectors and cooperators can be at an evolutionary disadvantagecompared to a neutral mutant.The VT model behaviour is qualitatively very similar to that of the well-mixedpopulation. The major difference is that the full spectrum of behaviour is availablefor a much smaller range of b/c values for the VT model. This means that coopera-tion is successful at smaller benefit-to-cost ratios, as is consistent with our previousfindings. It should be noted however, that these classifications are often approxi-mate for the VT model. We observe, in a number of cases, additional fixed pointsvery close to n = 0 and n = Z . It is also clear from Figure 11 that the coexistence& coordination behaviour is much less pronounced that it is for the well-mixed case,with the internal fixed points much closer to the boundaries. Population structure has long been established as a mechanism for promoting theevolution of cooperation [44]. If interactions are limited to an individual’s neigh-bourhood, then cooperators can form mutually beneficial clusters. However, thesuccess of cooperation is also influenced by the update rule. Results for the cyclegraph in Section 3 demonstrate that a global update rule, such as the shift update,can lead to less stringent conditions for cooperation to be favoured when comparedto local update rules. Within the local update rules there are also clear differences:cooperation tends to fare better with a death-birth update rule than a birth-death.In fact for the birth-death update on a cycle, the condition for cooperation to be19 .0000.005 G ( n ) VT h = 0.0 h = 0.2 h = 0.5 h = 0.8 s = h = 1.00.0000.005 G ( n ) s = G ( n ) s = b=1.53.50.0000.005 G ( n ) WM h = 0.0 h = 0.2 h = 0.5 h = 0.8 s = h = 1.00.0000.005 G ( n ) s = G ( n ) s = b=311 Figure 11: Gradient of selection, G ( n ) for a PGG with logistic benefit function in a well-mixed population (WM) and Voronoi tessellation model with decoupled update (VT). Thequalitative behaviour is very similar between the two, however occurs at different values ofthe benefit, b/c . Where G ( n ) > , selection is working to increase n , and vice versa. Theroots of G ( n ) = 0 can be considered as fixed points, and we can use these to classify thebehaviour in different parameter regions. s and h of the logistic benefit function,as well as the benefit-to-cost ratio. In general, a lower benefit-to-cost ratio is re-quired for cooperative success for the VT model, than the well-mixed population.In other words cells need a lower incentive to cooperate. This is consistent with ourexpectations: both models use global updating, however the population structurein the VT model allows for positive assortment of cooperators.Although cooperation is more successful in the VT model, than the well-mixedpopulation, the qualitative behaviour is very similar. We have characterised theevolutionary dynamics by considering conditions for which cooperation is beneficialand/or favourable, as well as calculating the gradient of selection.As long as the steepness, s , is large enough, we tend to see coexistence behaviourwhen the inflection point, h , is less than a half and coordination behaviour whenit is greater. These regimes are characterised by the fixed points of the gradient ofselection. They also correspond to the regions in parameter space where cooperationis beneficial, but not favourable (coexistence), and favourable, but not beneficial(coordination). For small steepness, the game approaches an NPD and there isdominance behaviour. In this regime, conditions for cooperation to be beneficialand favoured coincide.Examining the gradient of selection enables us to identify an additional dynam-ical regime: mixed coexistence & coordination, which occurs in a region around h = 0 . , as long as s and b/c are sufficiently large. This regime is characterised bytwo stable fixed points, one corresponding to all-defection, and the other to a het-erogenous, majority-cooperator state. This dynamic has been identified previouslyfor both well-mixed populations [39] and graph-structured populations with localupdating [22]. We have shown that it can also occur for the VT model, however theinternal fixed points tend to be much closer to the boundaries.It is beyond the scope of this paper to consider the full dynamics for an ep-ithelial population structure with local update rules. However, we have consideredconditions for cooperation to be favourable on a hexagonal lattice with death-birthupdate rule using results from [21]. We found the critical benefit-to-cost ratios tobe intermediate between the well-mixed population and VT model. This is consis-tent with previous results for two-player games [23]. Taken together, these resultssuggest not only that population structure promotes cooperation, but that globalupdating also plays a crucial role.It is worth taking a moment to consider the implications of beneficial andfavourable mutations for invasion, and how we distinguish between the two con-cepts. Whether or not a mutation is beneficial is perhaps the most relevant measurefor a single invasion event. It essentially tells us that the mutated cell has a higherprobability of invasion in a wild-type population than a wild-type cell would have,and therefore it has an evolutionary advantage. The significance of a mutation beingfavourable is a little less clear, as it compares two different invasion processes: theprobability of invasion of a mutated cell in a wild-type population is higher than the21onverse scenario, where a wild-type cell invades a population of mutants. However,the condition for a mutant to be favoured is also equivalent to the condition forcooperation to dominate, if mutation is allowed.These results suggest that there are a number of important considerations whichneed to be taken into account when modelling the evolution of cooperation in ep-ithelia: population structure, update rule and benefit function. In particular, thesuccess of cooperation could easily be underestimated, either by ignoring populationstructure, or by including population structure, but considering only local updaterules. If birth and death are decoupled in real epithelia, our results suggest thatproducers of growth factor can have an evolutionary advantage, even when the costto themselves is relatively high compared to the benefit.These local and global update rules represent two extremes: absolute spatialcoupling of birth and death, or absolute decoupling, respectively. In reality, therelationship between cell division (birth) and apoptosis/differentiation (death) inan epithelium, is likely somewhere in between the two. Regulatory processes suchas contact inhibition [45] and crowding-induced extrusion [46, 47] could lead to aweaker form of spatial coupling, somewhere in between local and global updating.Furthermore, the form of coupling could look more like a death-birth process [48] ora birth-death process, depending on the specific dynamics of the epithelium. Under-standing the nature of this coupling is therefore crucial for predicting evolutionaryoutcomes in epithelia. Data accessibility.
The code and data can be accessed at https://github.com/jessierenton/pgg-epithelium
Authors’ contributions.
JR and KMP designed the research. JR carried out the re-search and wrote the paper. KMP edited the paper.
Acknowledgments.
This research was funded by an EPSRC studentship held by JR.
Competing interests.
We declare we have no competing interests.
References [1] D. Hanahan and R. A. Weinberg, “The hallmarks of cancer,”
Cell , vol. 100,no. 1, pp. 57–70, 2000.[2] D. Hanahan and R. A. Weinberg, “Hallmarks of cancer: The next generation,”
Cell , vol. 144, no. 5, pp. 646–674, 2011.[3] J. Jouanneau, G. Moens, Y. Bourgeois, M. F. Poupon, and J. P. Thiery, “Aminority of carcinoma cells producing acidic fibroblast growth factor induces acommunity effect for tumor progression.,”
Proceedings of the National Academyof Sciences of the United States of America , vol. 91, pp. 286–90, jan 1994.[4] R. Axelrod, D. E. Axelrod, and K. J. Pienta, “Evolution of cooperation amongtumor cells,”
Proceedings of the National Academy of Sciences of the UnitedStates of America , vol. 103, pp. 13474–9, sep 2006.[5] O. Warburg, “On the origin of cancer cells,”
Science , vol. 123, no. 3191, pp. 309–314, 1956.[6] M. Archetti, “Evolutionary dynamics of the Warburg effect: Glycolysis as acollective action problem among cancer cells,”
Journal of Theoretical Biology ,vol. 341, pp. 1–8, jan 2014. 227] I. P. M. Tomlinson, “Game-theory models of interactions between tumour cells,”
European Journal of Cancer , vol. 33, pp. 1495–1500, aug 1997.[8] D. Basanta and A. Deutsch, “A game theoretical perspective on the somatic evo-lution of cancer,” in
Selected Topics in Cancer Modeling , pp. 1–16, BirkhäuserBoston, 2008.[9] S. Hummert, K. Bohl, D. Basanta, A. Deutsch, S. Werner, G. Theißen,A. Schroeter, and S. Schuster, “Evolutionary game theory: cells as players,”
Molecular BioSystems , vol. 10, pp. 3044–3065, oct 2014.[10] M. Archetti and I. Scheuring, “Review: Game theory of public goods in one-shotsocial dilemmas without assortment,”
Journal of Theoretical Biology , vol. 299,pp. 9–20, apr 2012.[11] C. Hauert, S. De Monte, J. Hofbauer, and K. Sigmund, “Volunteering as RedQueen mechanism for cooperation in public goods games,”
Science , vol. 296,pp. 1129–1132, may 2002.[12] F. C. Santos, M. D. Santos, and J. M. Pacheco, “Social diversity promotes theemergence of cooperation in public goods games,”
Nature , vol. 454, no. 7201,pp. 213–216, 2008.[13] L. A. Bach, S. M. Bentzen, J. Alsner, and F. B. Christiansen, “An evolutionary-game model of tumour–cell interactions: possible relevance to gene therapy,”
European Journal of Cancer , vol. 37, pp. 2116–2120, nov 2001.[14] L. A. Bach, T. Helvik, and F. B. Christiansen, “The evolution of n-playercooperation - Threshold games and ESS bifurcations,”
Journal of TheoreticalBiology , vol. 238, pp. 426–434, jan 2006.[15] M. Archetti, “The volunteer’s dilemma and the optimal size of a social group,”
Journal of Theoretical Biology , vol. 261, pp. 475–480, dec 2009.[16] M. Archetti, “Cooperation as a volunteer’s dilemma and the strategy of con-flict in public goods games,”
Journal of Evolutionary Biology , vol. 22, no. 11,pp. 2192–2200, 2009.[17] M. Archetti, D. A. Ferraro, and G. Christofori, “Heterogeneity for IGF-II pro-duction maintained by public goods dynamics in neuroendocrine pancreaticcancer,”
Proceedings of the National Academy of Sciences of the United Statesof America , vol. 112, pp. 1833–1838, feb 2015.[18] M. Archetti and I. Scheuring, “Coexistence of cooperation and defection inpublic goods games,”
Evolution , vol. 65, no. 4, pp. 1140–1148, 2011.[19] H. Ohtsuki, C. Hauert, E. Lieberman, and M. A. Nowak, “A simple rule forthe evolution of cooperation on graphs and social networks,”
Nature , vol. 441,no. 7092, pp. 502–505, 2006.[20] M. A. Nowak, C. E. Tarnita, and T. Antal, “Evolutionary dynamics in struc-tured populations,”
Philosophical Transactions of the Royal Society B: Biolog-ical Sciences , vol. 365, no. 1537, pp. 19–30, 2010.[21] J. Peña, B. Wu, J. Arranz, and A. Traulsen, “Evolutionary Games of Multi-player Cooperation on Graphs,”
PLoS Computational Biology , vol. 12, no. 8,pp. 1–15, 2016. 2322] M. Archetti, “Cooperation among cancer cells as public goods games on Voronoinetworks,”
Journal of Theoretical Biology , vol. 396, pp. 191–203, may 2016.[23] J. Renton and K. M. Page, “Evolution of cooperation in an epithelium,”
Journalof the Royal Society Interface , vol. 16, p. 20180918, 2019.[24] E. Lieberman, C. Hauert, and M. A. Nowak, “Evolutionary dynamics ongraphs,”
Nature , vol. 433, no. 7023, pp. 312–316, 2005.[25] C. G. Nathanson, C. E. Tarnita, and M. A. Nowak, “Calculating EvolutionaryDynamics in Structured Populations,”
PLoS Computational Biology , vol. 5,p. e1000615, dec 2009.[26] J. Zukewich, V. Kurella, M. Doebeli, and C. Hauert, “Consolidating birth-deathand death-birth processes in structured populations,”
PLoS ONE , vol. 8, no. 1,p. e54639, 2013.[27] N. Masuda, “Directionality of contact networks suppresses selection pressure inevolutionary dynamics,”
Journal of Theoretical Biology , vol. 258, pp. 323–334,may 2009.[28] T. Antal, H. Ohtsuki, J. Wakeley, P. D. Taylor, and M. A. Nowak, “Evolutionof cooperation by phenotypic similarity,”
Proceedings of the National Academyof Sciences of the United States of America , vol. 106, pp. 8597–8600, may 2009.[29] C. E. Tarnita, T. Antal, H. Ohtsuki, and M. A. Nowak, “Evolutionary dynamicsin set structured populations.,”
Proceedings of the National Academy of Sciencesof the United States of America , vol. 106, pp. 8601–8604, may 2009.[30] B. Allen and M. A. Nowak, “Evolutionary shift dynamics on a cycle,”
Journalof Theoretical Biology , vol. 311, pp. 28–39, oct 2012.[31] A. Pavlogiannis, K. Chatterjee, B. Adlam, and M. A. Nowak, “Cellular coop-eration with shift updating and repulsion,”
Scientific Reports , vol. 5, dec 2015.[32] F. A. Meineke, C. S. Potten, and M. Loeffler, “Cell migration and organizationin the intestinal crypt using a lattice-free model,”
Cell Proliferation , vol. 34,no. 4, pp. 253–266, 2001.[33] W. Maciejewski, F. Fu, and C. Hauert, “Evolutionary Game Dynamics in Pop-ulations with Heterogenous Structures,”
PLoS Computational Biology , vol. 10,no. 4, p. e1003567, 2014.[34] P. A. P. Moran, “Random processes in genetics,”
Mathematical Proceedings ofthe Cambridge Philosophical Society , vol. 54, pp. 60–71, jan 1958.[35] B. Wu, A. Traulsen, and C. Gokhale, “Dynamic Properties of EvolutionaryMulti-player Games in Finite Populations,”
Games , vol. 4, pp. 182–199, may2013.[36] A. Traulsen, C. Hauert, and H.-G. Schuster, “Stochastic evolutionary gamedynamics,” in
Reviews of Nonlinear Dynamics and Complexity (H.-G. Schuster,ed.), vol. 2, ch. 1, Wiley-VCH, 2009.[37] C. S. Gokhale and A. Traulsen, “Evolutionary games in the multiverse,”
Pro-ceedings of the National Academy of Sciences of the United States of America ,vol. 107, pp. 5500–5504, mar 2010. 2438] C. E. Tarnita and P. D. Taylor, “Measures of relative fitness of social behaviorsin finite structured population models,”
American Naturalist , vol. 184, no. 4,pp. 477–488, 2014.[39] M. Archetti, “Evolutionary game theory of growth factor production: implica-tions for tumour heterogeneity and resistance to therapies,”
British Journal ofCancer , vol. 109, pp. 1056–1062, aug 2013.[40] M. Archetti, “Dynamics of growth factor production in monolayers of cancercells and evolution of resistance to anticancer therapies,”
Evolutionary Appli-cations , vol. 6, pp. 1146–1159, dec 2013.[41] I. M. M. Van Leeuwen, G. R. Mirams, A. Walter, A. Fletcher, P. Murray, J. Os-borne, S. Varma, S. J. Young, J. Cooper, B. Doyle, J. Pitt-Francis, L. Momta-han, P. Pathmanathan, J. P. Whiteley, S. J. Chapman, D. J. Gavaghan, O. E.Jensen, J. R. King, P. K. Maini, S. L. Waters, and H. M. Byrne, “An integrativecomputational model for intestinal tissue renewal,”
Cell Proliferation , vol. 42,no. 5, pp. 617–636, 2009.[42] J. M. Pacheco, F. C. Santos, M. O. Souza, and B. Skyrms, “Evolutionarydynamics of collective action in N-person stag hunt dilemmas,”
Proceedingsof the Royal Society B: Biological Sciences , vol. 276, no. 1655, pp. 315–321,2009.[43] M. O. Souza, J. M. Pacheco, and F. C. Santos, “Evolution of cooperation underN-person snowdrift games,”
Journal of Theoretical Biology , vol. 260, pp. 581–588, oct 2009.[44] M. A. Nowak, “Five rules for the evolution of cooperation,”
Science , vol. 314,no. 5805, pp. 1560–1563, 2006.[45] A. I. Mcclatchey and A. S. Yap, “Contact inhibition (of proliferation) redux,”
Current Opinion in Cell Biology , vol. 24, no. 5, pp. 685–694, 2012.[46] G. T. Eisenhoffer, P. D. Loftus, M. Yoshigi, H. Otsuna, C.-b. Chien, P. A.Morcos, and J. Rosenblatt, “Crowding induces live cell extrusion to maintainhomeostatic cell numbers in epithelia,”
Nature , vol. 484, no. 7395, pp. 546–549,2012.[47] R. Fernandez-Gonzalez and J. A. Zallen, “Feeling the squeeze: Live-cell extru-sion limits cell density in epithelia,” may 2012.[48] K. R. Mesa, K. Kawaguchi, K. Cockburn, D. Gonzalez, J. Boucher, T. Xin,A. M. Klein, and V. Greco, “Homeostatic Epidermal Stem Cell Self-Renewal IsDriven by Local Differentiation,”
Cell Stem Cell , vol. 23, no. 5, pp. 677–686.e4,2018.[49] J. M. Osborne, A. G. Fletcher, J. M. Pitt-Francis, P. K. Maini, and D. J.Gavaghan, “Comparing individual-based approaches to modelling the self-organization of multicellular tissues,”
PLOS Computational Biology , vol. 13,no. 2, p. e1005387, 2017. 25 ppendix A Voronoi tessellation model
The Voronoi tessellation (VT) model was introduced in [32, 41]. We use the versionand parameter values from [23] in all simulations in this paper. Parameters are givenin Table 1 and the model is defined as follows.
Table 1: Table of parameters used in the Voronoi tessellation model [23, 49].
Parameter Description Value µ spring constant 50 s natural cell separation 1 (cid:15) initial sibling cell separation 0.1 η drag coefficient 1 ∆ t time-step (hours) 0.005 λ rate of division/death (hours − ) − The VT model represents a tissue as a set of points points in a fixed domain withperiodic boundary conditions. Each point corresponds to a cell-centre and movessubject to spring-like forces cells exert on their neighbours.We define the force cell j exerts on its neighbour i to be F ij ( t ) = − µ ˆ r ij ( t )( | r ij ( t ) | − s ij ( t )) , (55)where µ is the spring constant, r ij = r i − r j is the vector pointing from cell j to cell i , and s ij is the natural separation between cells i and j . This is set be aconstant s ij = s , with the exception that newborn sibling cells have a separation (cid:15) immediately after division. For these cells, s ij grows linearly over the course of anhour to reach s . The total force on a cell i is given by F i ( t ) = (cid:88) j ∈N i ( t ) F ij , (56)where N i ( t ) denotes the set of cells neighbouring i .It is assumed that motion is over-damped and thus the equation of motion for i is the first order differential equation η d r i dt = F i ( t ) , (57)where η is the damping constant. This is solved numerically using r i ( t + ∆ t ) = r i ( t ) + ∆ tη F i , (58)where the time-step, ∆ t , is sufficiently small to ensure numerical stability.At each time-step, after cells have moved, a Voronoi tessellation is performed.This partitions the domain into polygonal regions, each corresponding to the shape ofa cell. It also defines the neighbourhood connections, which are needed to determinethe forces cells exert on one another, as well as cell fitnesses.The population evolves through a process of sequential cell divisions and deaths.In order to keep population size constant, cell divisions and deaths are constrainedto occur simultaneously. These update events occur at rate λ , according to a con-tinuous time Moran process. When an update occurs, a cell is chosen to divide with26robability proportional to fitness. This parent cell is removed from the tissue andreplaced with two identical progeny cells, separated by a distance (cid:15) , across a randomaxis. Simultaneously, a cell is chosen to die uniformly at random, and is removedfrom the tissue. After an update event the neighbourhood structure is defined onceagain by performing a Voronoi tessellation. Appendix B Neighbour distributions in the VT model
In Sections 3 and 4 we derived conditions under which cooperation is favoured andbeneficial, given by Equations (41) and (47) respectively. These derivations arebased on the assumption that the frequency of cells with k neighbours is a fixeddistribution, g k , independent of the cell type or the number of cooperators in thepopulation, n . Figure 12: Neighbour distributions in the VT model for cooperators (C) and defectors (D),for varying cooperator population size, n . Data is generated from simulations with totalpopulation size Z = 100 in the neutral selection limit, δ = 0 . Figure 12 plots neighbour distributions from simulations of the VT model forcooperators and defectors at different values of n . It is clear from the plot thatthe assumption is a reasonable one. The neighbour distributions are approximatelyequal for different values of n and for the two cell types. The exception is whenthere are either very few cooperators or very few defectors, i.e. near n = 1 and n = 99 respectively. In the case where there is only one or very few cooperators,the cooperator neighbour distribution becomes slightly more narrow. The converseis true when there are few defectors. Appendix C Minimising the critical benefit-to-costratio at which cooperation is favoured
In Section 4 we considered conditions for cooperative success for a sigmoid benefitfunction, as defined by Equation (4). It is clear from Figure 5 that the critical27enefit-to-cost ratios, ( b/c ) ∗ , at which ρ C = ρ D , are minimised at h = 0 . , andsymmetric across that point. This appears to hold for both the Voronoi tessellationmodel with decoupled update, and for the death-birth update on a fixed hexagonallattice. In the following we show that this is indeed true for any system where < s < ∞ and the structure coefficients, σ j , are increasing for ≤ j < k .We rewrite Equation (51), defining ( c/b ) ∗ , such that cooperation is favoured for c/b < ( c/b ) ∗ (cid:16) cb (cid:17) ∗ = 1 Z − k (cid:88) j =0 σ j (cid:20) β (cid:18) j + 1 k + 1 (cid:19) − β (cid:18) k − jk + 1 (cid:19)(cid:21) . (59)We have assumed that the number of neighbours, k , is fixed, however the results areeasily generalisable to variable k . Defining γ j = β (cid:18) j + 1 k + 1 (cid:19) − β (cid:18) k − jk + 1 (cid:19) (60)we obtain (cid:16) cb (cid:17) ∗ = 1 Z − σ k + (cid:88) k>j ≥ k/ ( σ j − σ k − j − ) γ j . (61)By taking derivatives with respect to h we show that for k/ ≥ j < k , γ j ( h ) ismaximised when h = 0 . . In order that this corresponds to a unique maximum of ( c/b ) ∗ , and thus a minimum of the critical benefit-to-cost ratio, certain conditionson σ j must be satisfied.First we show that γ j has one extremum at h = 0 . for < s < ∞ . Wesubstitute Equation (4) into Equation (60) and take the first derivative with respectto h , letting r = j +1 k +1 . Thus we obtain dγ j dh = ddh (cid:20) (1 + e s ( h − r ) ) − − (1 + e s ( h + r − ) − (1 + e s ( h − ) − − (1 + e sh ) − (cid:21) (62) = ddh (cid:20) e s ( r − − e − sr − e − s · (1 + e s ( h − )(1 + e sh )(1 + e s ( h − r ) )(1 + e s ( h + r − ) (cid:21) (63) = s · e s ( r − − e − sr − e − s · e sh (1 + e − s − e − sr − e s ( r − )(1 − e s (2 h − )(1 + e s ( h + r − ) (1 + e s ( h − r ) ) . (64)Setting dγ j /dh = 0 , gives one root at h = 0 . , for < s < ∞ . This is a uniquestationary point of ( c/b ) ∗ so long as there is at least one value of j ∈ [ k/ , k ) forwhich ( σ j − σ k − j − ) (cid:54) = 0 . We can show that this is a maximum by considering thesecond derivative at h = 0 . d γ j dh (cid:12)(cid:12)(cid:12)(cid:12) h = = − s · e s/ (1 + e − s − e − sr − e s ( r − )( e s ( r − − e − sr )(1 − e − s )(1 + e s ( r − / ) (1 + e − s ( r − / ) (65)which is negative given that / < r < . This corresponds to ( k − / < j < k ,encompassing all the values of j which we sum over in Equation (61). Therefore, inorder that ( c/b ) ∗ is maximised when h = 0 . , we require that ( σ j − σ k − j − ) ≥ for k/ ≤ j < k and non-zero for at least one value of j in the range. This condition isguaranteed if σ j is an increasing, but not constant, function for ≤ j < k .It is clear from Figure 4 that σ j +1 ,k > σ j,k ∀ j, k for the VT model with decoupledupdate, therefore h = 0 . maximises ( c/b ) ∗ in this case. For k -regular graphs withdeath-birth update rule, we can verify whether this is true by using the approximateexpressions for the structure coefficients derived in [21]. These are plotted for various28 values in Figure 13. For smaller values of k , we can see that σ j is strictly increasingfor ≤ j < k . However, as k increases, a growing region appears for which σ j isconstant. So long as there is at least one value of j < k for which ( σ j − σ k − j − ) > , ( c/b ) ∗ is maximised at h = 0 . . However, as k → ∞ , we approach the case where σ j are constant for j < k , and we regain the well-mixed population result that ( c/b ) ∗ are independent of h . j k 510204080 Figure 13: Structure coefficients, σ j , for k -regular graphs with death-birth update rule[21]. It is clear that σ j is increasing (or constant) for j < k . Thus far we have limited ourselves to the case where < s < ∞ . In the limit s → , we obtain an NPD game with a linear benefit function which is independentof h . The value of ( c/b ) ∗ therefore does not depend on h either, as can be seen inFigure 5. In the limit s → ∞ , the VD game is approached and the benefit functionceases to be continuous. In this case the unique maximum at h = 0 . is maintainedonly if σ j are strictly increasing, and therefore ( σ j − σ k − j − ) > . This is true for theVT model with decoupled update and for k -regular graphs with death-birth update,if k is sufficiently small. On the other hand, if ( σ j − σ k − j − ) = 0 for some values of j ∈ [ k/ , k ) , h = 0 . ceases to be an isolated maximum, and there is a region of h values, around h = 0 . , which maximise ( c/b ) ∗1