Evolutionary dynamics and the evolution of multiplayer cooperation in a subdivided population

The classical models of evolution have been developed to incorporate structured populations using evolutionary graph theory and, more recently, a new framework has been developed to allow for more ﬂexible population structures which potentially change through time and can accommodate multiplayer games with variable group sizes. In this paper we extend this work in three key ways. Firstly by developing a complete set of evolutionary dynamics so that the range of dynamic processes used in classical evolutionary graph theory can be applied. Secondly, by building upon previous models to allow for a general subpopulation structure, where all subpopulation members have a common movement distribution. Subpopulations can have varying levels of stability, represented by the proportion of interactions occurring between subpopulation members; in our representation of the population all subpopulation members are represented by a single vertex. In conjunction with this we extend the important concept of temperature (the temperature of a vertex is the sum of all the weights coming into that vertex; generally, the higher the temperature, the higher the rate of turnover of individuals at a vertex). Finally, we have used these new developments to consider the evolution of cooperation in a class of populations which possess this subpopulation structure using a multiplayer public goods game. We show that cooperation can evolve providing that subpopulations are suﬃciently stable, with the smaller the subpopulations the easier it is for cooperation to evolve. We introduce a new concept of temperature, namely “subgroup temper-ature”, which can be used to explain our results.


Introduction
Evolutionary game theory has proved to be a very successful way of modelling the evolution of, and behaviour within, populations.The classical models Email addresses: karan.pattni.1@city.ac.uk (Karan Pattni), mark.broom@city.ac.uk (Mark Broom), rychtar@uncg.edu(Jan Rychtář) mainly focused on well-mixed populations playing two player games [31,30], or alternatively playing games against the entire population [30].Simple models such as the Hawk-Dove game [29] and the sex ratio game [20] have been used to explain important biological phenomena.These models were developed to consider finite populations explicitly [34, (although see [32,33] for important earlier non-game theoretic work) and structured populations using the now widespread methodology of evolutionary graph theory originated in [26] (see also [3,9,52,27], and [1,44] for reviews).Such population structures can have a profound effect on the result of the evolutionary process even when individuals have a fixed fitness [26,28,40].
Further, even for a given structure, the rules of the evolutionary dynamics have a significant effect on the evolution of the population.
Previous work has investigated a number of important questions, the most widely considered being how cooperation can evolve.The evolution of cooperation, where individuals make sacrifices to help others, can seem paradoxical within the context of natural selection, especially amongst unrelated individuals.There are a number of ways that mathematical modelling has demonstrated that cooperation can occur [35]; one key way is through the presence of population structure, which can mean that cooperative individuals are more likely to interact with other cooperators, which makes them resistant to exploitation by defectors [36,42].In particular, this is true for structures where individuals are heterogeneous [43] allowing hubs or clusters of cooperators to form.The dynamics that one uses are also important; for example [36] showed that death-birth or birth-death dynamics with selection on the second event promotes cooperation but not when selection happens in the first event.
One limitation of evolutionary graph theory is that it naturally lends itself to pairwise games, whereas real populations can often involve the simultaneous interaction of many individuals [45,15].Multiplayer games, whilst more common in economic modelling [21,6], have become used in increasing frequency within evolutionary games starting with [38,7] (see also [14,18]) and it is important to incorporate these too into the modelling of structured populations.
A multiplayer public goods game [4,5,19,54], (and this type of game is central to our paper too, see Section 2.2) has been used in evolutionary graph theory [25,51,24,41,56], but this typically involves forming an individual and all of its neighbours into a group and allowing them to play a game.Although this is convenient, it is not really natural because there is no mechanism for deciding how individuals spend their time, and so how they share that time with others, either singly or in groups.
More recently a general framework has been developed [10,13,8,11] which considers the interaction of populations in a more flexible way, where groups of any size can form, with different propensity potentially depending upon a number of factors, including the history of the process.Crucially, the key elements of evolutionary graph theory of population structure, game and evolutionary dynamics occur for this new framework too; this makes it capable of analysing different spatial structures whilst providing the flexibility for different multiplayer interactions.Prior to the current paper, the actual applications of the The fully independent model from [10].There are N individuals who are distributed over M places such that In visits place Pm with probability pnm.Individuals interact with one another when they meet, for example, I 1 and I 2 can interact with one another when they meet in P 1 .
above framework have been limited.In particular only a single evolutionary dynamics (the BDB dynamics from the current paper) has been used, and only relatively simple populations, which resembled those in evolutionary graph theory (the population consisting of individuals each resident at a unique graph vertex) have been considered.
In this paper we further develop the general theory of the framework originated in [10].We first show how to represent subpopulations using a reduced graphical representation within our structure, which will then allow us to potentially consider larger populations with a richer structure than previously.We then demonstrate how to apply a standard set of evolutionary dynamics to consider a range of evolutionary processes.This is vital since, as mentioned above, dynamics can have a big effect on the outcome of evolution within other models, including evolutionary graph theory, and as we will see, this is certainly also true for our work.Finally we use these new tools to consider the evolution of cooperation using a multiplayer public goods game [51,48,49,4] and show that cooperation can occur when both the structure and evolutionary dynamics act together in favour of the cooperators.
The paper is structured as follows: in Section 2 the model framework is described, including how to incorporate subpopulations.In Section 3 a standard set of evolutionary dynamics to be used with our model are defined.In Section 4 we introduce and discuss the important concepts of fixation probability and temperature.In Section 5 we study the evolution of cooperation in our model with subpopulations.Section 6 is then a general discussion.

A framework for modelling evolution in structured populations
A framework for modelling the movement of individuals was presented in [10].This is a very general and flexible methodology, the details of which are not necessary for the current paper.Below we describe the fully independent version of this framework in which individuals move independently of each other and Payoff to cooperator in a group (including itself) of c cooperators and d defectors.
State in which all individuals are cooperators.
Mean fixation probability of a cooperator.
Weighted adjacency matrix that represents an evolutionary graph.v n Vertex n of an evolutionary graph.independently of the population's history (any past movements), and a version of the fully independent model called the territorial raider model as introduced in [10] and further developed in [8].We then develop a generalization of this model, which then forms the basis of much of the work in this paper, although we note that Section 3 in particular is more general.Important terms used in the current paper are given in Table 1.

The population structure
We begin by introducing the fully independent model.Consider a population made up of N individuals I 1 , . . ., I N who can move around M places P 1 , . . ., P M .
The probability of individual I n being at place P m is denoted by p nm ; see Figure 1 for a visual representation using a bi-partite graph.When individuals move around they form groups.Let G denote any group of individuals, then the probability χ(m, G) that group G forms in place P m is given by We can show from equation (2.1) that This follows intuitively from the fact that individual I n has to be present in some place P m in some group G at any given time.The mean size of an individual's group (see also [13]) is given by where the simplification of the denominator follows from equation (2.2).
When a group of individuals is formed they will then interact with one another.In particular, individual I n will receive a payoff that depends upon the group G it is present in and the place P m occupied by this group.This is denoted as R n,m,G and was referred to in [10] as a direct group interaction payoff because individual I n only interacts with other individuals with whom it is directly present ( [10] allowed for a more general class of payoff but this is the only type we will consider, and hence will just refer to it as the payoff).
Individual I n 's fitness is then calculated by averaging its payoffs over all possible groups and places that these groups can form as follows: (2.4) We now move on to consider the territorial raider model.In the territorial raider model, each individual I n has its own place P n with no unoccupied places and, therefore, there is a one-to-one correspondence between individuals and places.A graph is used to represent the structure of the population where each vertex represents an individual and its corresponding home such that two connected individuals can raid each other's home places (see Figure 2).The probability of raiding another's home place is governed by a common movement parameter called home fidelity, h, that measures an individuals' preference for their home place.In particular, an individual with d neighbours would stay on their home place with probability h/(h + d) or raid any one of its neighbours' home places with an equal probability of 1/(h + d) (see Figure 2).We now generalise the territorial raider model to include subpopulations, based upon their movement distributions.We will see that individuals within a given subpopulation are more likely to interact with each other than with members of other subpopulations, and this will affect the success of their strategies.
Consider the fully independent model.We define a subpopulation of individuals as a division of individuals from the main population that is well-mixed [10], which simply means that all of these individuals have an identical distribution over the places.In particular, for a subpopulation Q we have that p im = p jm ∀ i, j ∈ Q and m = 1, . . ., M .This can be visualised in terms of a bipartite graph as in Figure 1 where the vertices are now occupied by subpopulations rather than individuals.This subpopulation structure is thus a special case of the fully independent model.
For simplicity we will assume that individuals move as they do in the territorial raider model; thus our model is a generalization of the territorial raider We will assume that individuals in subpopulation Q m treat place P m as their home place, so that there is a one-to-one correspondence between subpopulations and places.However, because we allow subpopulations to be empty, we can have places in which no individuals reside.As before, the movement probabilities of the individuals is governed by the home fidelity h.In particular, a subpopulation Q m that can visit d neighbouring places will stay in home place P m with probability h/(h + d) or move to one of its neighbouring places with probability 1/(h + d).Note that when there is one individual in each subpopulation, that is we recover the original territorial raider model.This information can be visually represented in two different ways as shown in Figure 3, which includes a graph whose vertices represent both subpopulations and places.This generalized territorial raider model will be the basis of our detailed investigation of the evolution of cooperation in Section 5.

A multiplayer public goods game
A multiplayer Hawk-Dove game [46] and a public goods game were considered in [8], though there are other games that can be considered like the multiplayer stag hunt game [37].
In this paper we focus only on the multiplayer public goods game based on the game defined by [51], where an individual's payoff is an average of two player public goods games (just a version of the standard prisoner's dilemma) played with each of its group mates.Players can either cooperate (C) or defect (D).
A cooperator always pays a cost 1 so that the other player receives a reward v and a defector pays no cost but only receives a reward when present with a cooperator.Note that the cost is set to 1 because scaling all the payoffs by some other cost value does not affect the outcome of the game and, therefore, the reward v is a multiple of the cost.The payoff matrix is thus given by In [51] and most models involving public goods games, individuals are never alone, and so what happens in the case they are alone is not considered.However, in our case it is possible for an individual to be alone, for example, an individual could remain on its home place and not be raided.As in [8], we will assume that a lone cooperator still pays a cost but does not receive a reward and lone defectors receive nothing.There are other ways that we can allocate rewards to lone individuals; for example, in [22] there is a specific strategy, the loner strategy, where cooperators choose to be alone and not pay a cost.Our Table 2: Dynamics defined using the replacement weight as in [40].In each case, B (D) is appended to the name of the dynamics if selection happens in the birth (death) event.
choice seems a natural generalisation of the prisoners dilemma model [51], where individuals pay a cost but do not benefit from their own contributions.We note that our version makes cooperation harder to evolve than the alternatives.Thus if cooperators thrive in a population using our model, this can be thought of as strong support for the evolution of cooperation.
In the multiplayer public goods game, the payoffs to cooperators and defectors playing within a group of c cooperators and d defectors (including themselves) is then respectively given by where r is a background payoff, which is also a multiple of the cost, that every individual receives, representing the contribution from activities that are not related to the games.Generally, the effect of selection is weaker the larger the value of r (for example, see [12], Chapter 2).The payoff is then given by ) when I n is a cooperator (defector) and |G| = c+d, which can then be substituted into Equation 2.4 to find the individual's fitness.Note that here the payoffs do not depend upon the place occupied by the individuals, that is, R n,m,G ≡ R n,G .

Evolutionary dynamics
In this section we revisit the standard dynamics of evolutionary graph theory, before demonstrating how we can adapt each of them to our framework.For the current work there will actually only be two distinct dynamics, but for more general cases each will be distinct, and so it is important to consider them all.
We start by recalling the dynamics from evolutionary graph theory.

Evolutionary dynamics in evolutionary graph theory
An evolutionary graph [26,40]  For each of these dynamics, natural selection can influence the birth ('B' appended to name) or death ('D' appended to name) event.We use the definitions of [28] who extensively studied a set of each of these dynamics.In terms of the exact formulae of the transition probabilities, we use those of [40] as summarised in Table 2.In these definitions, the dynamics are a function of the replacement structure W and the fitnesses of the individuals such that the individual on vertex v n has fitness F n .

Evolutionary dynamics in our framework
In [8] a birth-death dynamics was defined to be used with the territorial raider model.In this section we shall develop a consistent set of dynamics for our framework.In particular, we will show that we can adapt the above dynamics widely used in evolutionary graph theory.
To consider the evolution of the population it is useful to think of the individuals in the population in an abstract way.In particular, individuals in the population change through time and, therefore, it is better to think of I i as a position that an individual can occupy.These positions are referred to as I-vertices in [8] and have a particular relationship to the places, although as the population evolves the actual individual, and in particular the type of individual, occupying the position may change.We will generally simply refer to these I-vertices as "individuals" but make the distinction where necessary.
This leads to a natural way to create evolutionary dynamics for our framework; namely, by mapping each individual I i to vertex v i , we can incorporate the replacement weights of different interaction methods straight into the formulae from Table 2.All that remains is to choose the replacement weights appropriately.
The replacement weights used here are based on the assumption that an offspring of individual I i is likely to replace another individual I j proportional to the time I i and I j spend together.The offspring of I i can also replace I i itself and it does this proportional to the time I i spends alone.Therefore, when i = j, the probability that I i and I j meet is given by summing χ(m, G) over all m such that i, j ∈ G.When they meet, we assume that I i will spend an equal amount of time with each other individual in group G and, therefore, weight χ(m, G) with 1/(|G| − 1) since there are |G| − 1 other individuals (an alternative weighting could be 1/|G| that allows interaction within groups larger than one to contribute to the probability of I i 's offspring replacing itself).Note that this is consistent with the payoffs from our public goods game, where each pairwise payoff equally contributes to the total payoff an individual receives.On the other hand, when i = j, we sum χ(m, G) over all m such that G = {i}.Here there is no need to weight χ(m, G) because I i is alone.
The replacement weights are therefore calculated as follows Thus we have a new set of evolutionary dynamics which can be applied to our framework in a wide variety of situations (including those that we consider later in this paper).Note that the dynamics used in [8] is the BDB dynamics defined from the above process.
By our definition W is symmetric, that is w ij = w ji ∀i, j, because the probability of I i meeting I j within any given group is clearly the same as that of I j meeting I i .We also have that W is doubly stochastic, that is 1 = j w ij = i w ij for all i, j, because w ij is the proportion of time I i spends with I j (with w ii the proportion of time it spends alone), and it is always in precisely one of these N categories.In this case, W is referred to as being isothermal [26,40].
We note that the results above hold because of the particular weights w ij that we have chosen.Although these are natural, they are not the only possibility.
In particular we could have alternative weights where w ij and w ji are not in general equal and/or where W is not isothermal.

The fixation probability
The (mean) fixation probability ρ C (ρ D ) is the probability that the offspring of a randomly placed mutant cooperator (defector) eventually replaces the entire population.This can be uniformly at random as in [26]; alternatively, one can use the mutant appearance distribution as described in [2].[8] used a version of this where they weighted the fixation probabilities using the mean temperature.
For this current work we use the arithmetic mean, as the difference between these two approaches is negligible here, with the arithmetic mean being greater than or equal to the weighted mean [2].For more details on how the fixation probability is calculated, see the Appendix.
As in [50], we will use the neutral fixation probability 1/N as a benchmark when comparing cooperators and defectors using their fixation probabilities.In particular, [50] say that selection opposes D replacing C when ρ C < 1/N and

Concepts of temperature
In [26] the in temperature (or just the temperature) of a vertex of an evolutionary graph was introduced to measure how likely an individual occupying a particular vertex is to be replaced by another individual's offspring.[28] extended this definition and introduced the out temperature of a vertex of an evolutionary graph to measure how likely the offpsring of the individual occupying that vertex will replace another individual.These definitions of the in and out temperatures of individual I n for an evolutionary graph W are respectively defined as follows In general, the in and out temperatures can be different.However, in our case, W is doubly stochastic and symmetric and, therefore, the in and out temperatures are identical.We therefore work with the definition of only in temperature and simply refer to it as the temperature.
An alternative version of the definition of temperature (used in [8]) is the strict temperature that measures how often an individual is likely to be replaced by other individuals excluding itself.Since W is doubly stochastic, the strict temperature of individual I n for an evolutionary graph W is given by The definition of strict temperature can be extended to subpopulations to give the strict subpopulation temperature.This measures how likely an individual in subpopulation Q m is to be replaced by an individual in another subpopulation.Clearly all individuals in a subpopulation have the same temperature (for any of our temperature definitions), since they all have the same movement distribution.The strict subpopulation temperature is calculated by summing all weights w ij such that I i is not part of subpopulation Q m and I j is part of subpopulation Q m giving  This means that if there is only one subpopulation then its strict subpopulation temperature is 0 by definition, that is, We note that a strategy introduced in one subpopulation can spread throughout the population because W is strongly connected.This implies that if there is more that one non-empty subpopulation then the strict subpopulation temperature is non-zero for all non-empty subpopulations, that is, To measure the connectedness of the subpopulations, that is how often the different subpopulations interact with one another, we use the mean strict subpopulation temperature that is defined as follows

Cooperation in generalized territorial raider models
In this section we study the effect that different model parameters have on the evolution of cooperation.For models investigating the evolution of cooperation using evolutionary graph theory, both the evolution and interaction As in Figure 4, cooperators evolve better in small groups (larger than 1), namely groups of size two and three, with a small advantage for groups of size four.
of individuals are dictated by a fixed structure, following games with a fixed number of players (almost always two).In our model the replacement structure emerges from the interactions between individuals, involving games with a varying number of players, and therefore give us a different perspective on the evolution of cooperation.
We note that no simulations were run to calculate the fixation probabilities in this paper, rather, all the states of the population were explicitly calculated following the procedure described in the Appendix.

The effect of the dynamics
As we mentioned in Section 1, for evolutionary graph theory models, cooperation is favoured when using DBB or BDD dynamics, but not DBD or BDB dynamics, if the structure allows a cluster of cooperators to form (also see [36]).This is consistent with [8] where we studied the effect of the BDB dynamics on the public goods game and cooperators generally performed poorly.It was shown that defectors dominate regardless of the structure of the population and the game parameters.We are now in a position to revisit the public goods game with more flexibility both in terms of the dynamics and the structure of the population.In terms of the dynamics, the results for BDB and DBD are identical (as are those for BDD and DBB), because the replacement structure W is symmetric and doubly stochastic, so whether birth or death occurs first (but not whether selection occurs in the first or second position) is irrelevant, see Table 2. Furthermore, the LB and LD dynamics are equivalent to the BDB and DBD dynamics, respectively, because W is isothermal.This can be shown for LB dynamics (and similarly for LD dynamics) as follows Thus in what follows, we only mention one dynamics from each pair, in each case the DB dynamics.
For DBD dynamics, the defectors do better than cooperators regardless of the population structure.However, for DBB dynamics, cooperators are favoured over defectors for certain population structures.In particular, these structures that favour cooperators contain small subpopulations, ideally of two individuals.
We can see this in Figure 4, where the fixation probability is plotted against different complete population structures for the DBD (Figure 4a) and DBB (Figure 4b) dynamics (as explained in the caption, for each population, each number in its representation corresponds to a subpopulation of that size).For example, for the complete structure 222 where there are 3 subpopulations of size 2, the cooperators outperform defectors by a large amount.
To understand why this is the case, consider a population of two individuals where one individual is a cooperator and the other a defector.Within such a population, the cooperator will be less fit than the defector.For DBD dynamics, the least fit individual is most likely to be chosen for death and the fixation probability is proportional to the fitness of the individual.This means that a cooperator has a low fixation probability compared to a defector.However, when using DBB dynamics, one of the two individuals in randomly chosen for death and immediately replaced by the offspring of the other individual.This means that regardless of the fitness of the individual, each type will fixate with probability 1/2.For sufficiently high home fidelity parameter h, individuals primarily interact with their members of their own subpopulation.Therefore, in such a population where there exists a subpopulation of two individuals, a cluster of two cooperators is more likely to form when using DBB dynamics.
This cluster of cooperators has a fitness larger than that of a cluster of defectors, provided that v > 1, thereby establishing a stronghold against defectors.In fact, a subpopulation of sufficiently small size (but greater than one) can establish a stronghold against defectors as shown in Figure 5.Here the fixation probability is plotted against a complete structure with four subpopulations that each have size ranging from 1 to 6. Subpopulations of size two are best for cooperation, with their advantage over defectors declining as the size of the subpopulation increases.Given the parameters used, subpopulations of two to four cooperators can successfully resist invasion, but larger subpopulations cannot.

The effect of the temperature
In [8] the strict temperature and mean group size were both shown to be strongly correlated with the fixation probability, with the effect of the former shown to be stronger.We therefore focus on the temperature, namely the strict subpopulation temperature.Note that in [8] there is one-to-one correspondence between individuals and places, which implies that the strict temperature and strict subpopulation temperature are identical, but this is not the case here.
The individual temperature is a measure of how often an individual interacts with other individuals including those who are part of the same subpopulation; thus an individual may have a high temperature but that does not mean it is interacting with individuals from other subpopulations.In particular whenever  In particular, we notice that the fixation probability of the cooperators is decreasing with the mean subpopulation temperature.
( individuals are not alone very often, this temperature does not vary so much between different individuals, and so is not a useful concept when there are nontrivial subgroups.The strict subpopulation temperature, on the other hand, considers interactions with individuals only from other subpopulations, and thus can be very variable.We shall see that this temperature is a good predictor of important population properties.
The mean strict subpopulation temperature decreases when home fidelity increases as shown in Figure 6a.This is because the individuals are more likely to remain on their home place than visit another place as home fidelity increases, therefore reducing interactions with other subpopulations, and in particular the probability that a member of one subpopulation replaces a member of another at any given time.
In [8] it was shown that for BDB dynamics for structures where each subpopulation is of size one, there was a linear relationship between the strict (subpopulation) temperature and the fixation probability, with the higher the temperature, the stronger the effect of selection.We investigated this for DBB dynamics, and found an opposite linear effect, which is consistent with [28] who showed that the DBB dynamics suppresses the effect of selection the most for the complete graph.We note that this relationship only holds for relatively weak selection, and we can reverse the relationship (and make it non-linear) by increasing the value of the reward.
To promote cooperation we need a structure involving a subpopulation of size at least two.However, whether these structures promote cooperation or not also depends upon the base fitness and reward, and so we assume that the base fitness and reward are sufficiently large for this to be the case, see Section 5.4.In this case, decreasing the temperature by increasing the home fidelity promotes cooperation.In particular, the relationship between the mean fixation probability of cooperators and the mean strict subpopulation temperature is negative and nonlinear as shown in Figure 6b.The nonlinearity arises not only from the nonlinear payoff function of the public good game, but also from the fact that there exists a subpopulation that has size at least two.For cooperators, the mean fixation probability is negatively correlated with the mean strict subpopulation temperature because the mean strict subpopulation temperature is highest when home fidelity is lowest, i.e. when cooperators cannot separate themselves from the population and form clusters, consequently defection evolves.On the other hand, for low mean strict subpopulation temperature, and so high home fidelity, it is easier to form clusters of cooperators that allows cooperation to evolve.This kind of behaviour is also evident in Figures 4 and     7.

The effect of the number of places
In [8] each individual had their home place and there were no empty places (non home places) that individuals could visit.In our case, individuals can visit non home places and we therefore investigate what effect this has on the evolution of cooperation.shows the effect of compensating for empty places by increasing the home fidelity such that the probability of staying in their home place, pnn, remains the same.We start at h = 30 for the 33 and 222 structures.As an empty place is added, h is increased so that pnn = 30/31 for the 330,. . .,330000 structures and pnn = 30/32 for 2220,. . .,222000 structures.In all cases r = 30 and v = 10.We can see that after compensating for the above effect, the influence of introducing empty places is both reversed and weakened.Figure (b) shows the mean strict subpopulation temperature dropping off when we compensate for the empty places by increasing the home fidelity such that pnn remains the same.
As seen in Figure 4, increasing the number of empty places that subpopulations can visit, whilst keeping all other parameters constant, makes it more difficult for cooperation to evolve.In particular, this effect is prominent for structures where cooperators were initially doing well.For example, for the structure 222 where the cooperators do best, increasing the number of places significantly reduces their fixation probability whilst increasing that of the defectors.Here increasing the number of places acts in the same way as decreasing the home fidelity, i.e. as decreasing the amount of time an individual spends in its home place with members of its subpopulation.Thus the amount of time an individual spends alone or with individuals not from its subpopulation increases, so that the overall fitness of a cooperative subpopulation will decrease (they still pay a cost but do not receive a benefit when alone).In terms of the dynamics, spending more time alone would increase the effect of selection in DBB dynamics because an individual with higher fitness that is randomly chosen for death is more likely to be replaced by its own offspring, which affects the cooperators adversely.A cooperative subpopulation will also have lower fitness because its members are more likely to interact with individuals from other subpopulations, therefore exposing them to defectors.The increased interaction between individuals will also increase the effect of selection in DBB dynamics because an individual with higher fitness that is randomly chosen for death is less likely to be replaced by an individual with lower fitness in the same subpopulation.
The increase in the number of places can be compensated for by increasing the home fidelity, so that individuals stay in their home place with the same probability.This has the effect of decreasing the mean strict subpopulation temperature as individuals are more likely to spend time with members of their subpopulation.This is shown in Figure 8, where we can see that the effect of adding empty places is now reversed, although the strength of this reverse effect is weak.
and is approximately 0 otherwise.Thus the fitness of an individual can be evaluated assuming that we have a group containing precisely all individuals from its subpopulation with probability 1. Due to the symmetric nature of our population, the weights for any two individuals in the same subpopulation will be the same, as will the weights for any two members of different subpopulations.
Denoting the latter as w O , which will be small, we have w ij = w O when I i and I j are not in the same subpopulation, and otherwise, with the probability of self-replacement negligible.
It follows that only replacements within subpopulations will happen, except very rarely.Thus we can assume that the battle within any mixed subpopulation of cooperator (C) and defector (D) individuals will be resolved with fixation of one type or the other before any new mixed subpopulation appears.
We thus consider a two stage process.Firstly, a new mixed group appears.
This occurs rarely, through the invasion of a cooperator into a defector subpopulation, or a defector into a cooperator subpopulation.Assuming that there are currently M C (M D = M − M C ) cooperator (defector) subpopulations, such a transition happens with probability of a cooperator into a defector subpopulation, or we obtain that the ratio of the two expressions in equations (5.1) and (5.2), and thus the relative frequency that the new invasions happen, is thus for large v and r.
The second process considers fixation within a well-mixed group of size L.
Following [23] we obtain the formula for the fixation probability of i cooperators within a population of size L.Here β k (δ k ) is the probability that the next event is the replacement of a defector (cooperator) by a cooperator (defector), when the number of cooperators is k.
We have here ) For sufficiently large r, we obtain where (5.8) The fixation probability of a single cooperator in a group of defectors is given by ρ C,L = x 1 , and the fixation probability of a single defector in a group of cooperators is ρ D,L = 1 − x L−1 .We thus have This implies that (5.10) Following our assumptions, the population evolves following a succession of invasions of subpopulations either of cooperators by defectors or of defectors by cooperators.The probability that the next such event will be the invasion of a subpopulation of defectors by a cooperator is simply where r S = p CI ρ C,L /p DI ρ D,L is the forward bias [40] of cooperative groups within our population.For a cooperator to fixate in the population it must first fixate within its group with probability ρ C,L , after which, there is a competition between groups proceeding precisely as in a Moran process, so that we have with the equivalent expression for ρ D , It is clear from equation (5.10) that r S > 1, so that ρ C is greater than ρ C,L (1 − 1/r S ) for any M .Letting M become large means that 1/N = 1/M L will be less than ρ C , but larger than ρ D , so that inequality (4.1) holds.This means that for sufficiently large h, r and v, we have that cooperation evolves for any given subpopulation size L. Thus cooperation can potentially evolve for arbitrarily large subpopulations, although as we have seen previously, it is easier for smaller subpopulations.

Discussion
In [10] a new framework for the flexible modelling of structured populations using multiplayer interactions was introduced, see also [8,13,11].This work built on classical evolutionary graph theory, but was limited in terms of the dynamics used.In this paper we have developed this framework further.Most importantly we have developed a full range of dynamics to apply in the framework, which will allow us to consider many different evolutionary scenarios.In particular these can be applied for the fully independent model in general, not just the examples considered here, enabling us to use a fuller range of the possibilities that our flexible framework allows.Thus this paper can be thought to complete the basic development phase of our work.
We have then developed the fully independent model to incorporate subpopulations and in particular consider a generalized version of the territorial raider model introduced in [8].This is beneficial because previously the fully independent model, represented in the bipartite graph in Figure 1, would require a vertex for every individual as well as an additional vertex for every available place.Now we just need a vertex per subpopulation, potentially allowing a small number of very large subpopulations to be considered, which would not have been possible previously.Thus this generalization allows us to look at much larger populations, which most real populations are, but still be able to use some analytical methods.The fact that larger populations can be considered without increasing complexity in turn allows us to incorporate other features that will enable our model to be applied more widely, as discussed below relating to mobile populations.This type of structure has been considered in a slightly different context, for example, the island-or community-structured populations of [53].In this model interactions occur at multiple levels, interactions between community members being more common than those with non-community members where interaction occurs at multiple levels.Members of one community first play a public goods game and then join the members of another community and play a public goods game such that, at the highest level, the entire population plays a public goods game.This is in contrast to our case, where individuals only play a game if they are present in the same place at the same time.They showed that cooperation can evolve when DBB dynamics are used and selection is weak within communities, which is consistent with our results.
We note that the framework of [8] is capable of modelling far wider behaviour than that developed here, in particular it is able to consider dynamic populations whose distributions continuously change due to their history, and the interactions that they have.Thus it can incorporate the type of situations with mobile populations modelled in [55,47].In particular, movement can follow a stochastic process in which the individuals move depending upon their current state as in [16].This is an important step in the development of realistic population models, for example related to territorial behaviour where animals can cover long distances, or movement behaviour varies throughout the year as seen in, for example, African wild dogs that live in packs [17].In a recently submitted paper [39] we have developed a Markov chain version of our model, and again consider a combination of theoretical developments and the specific application of the evolution of cooperation.This is our first step in the type of history-dependent analysis described above.
We then applied our new methodology to an example, considering the evolution of cooperation within a population involving subpopulations.We saw as in evolutionary graph theory that the choice of dynamics is crucial, and that DBD (and BDB) dynamics would not allow cooperation to evolve, but that DBB (and BDD) would, which is consistent with [36].Further, using the latter dynamics, the size and the level of isolation of the subpopulations is important, with the smaller the subpopulations and the greater the isolation, the greater the chance for cooperation to evolve.Unsurprisingly, the larger the level of reward v, the better the cooperators do.In particular, the larger the subpopulations, the larger the reward v required for cooperation to evolve; note that this is similar to the requirement that the benefit-to-cost ratio exceeds the average number of neighbours an individual has from [36].
We see from Figure 6 that our new idea of strict subgroup temperature is important in explaining the level of cooperation that evolves.Low (high) temperature helps promote the invasion of cooperators (defectors).In particular, higher temperatures allow cooperators to cluster more strongly and benefit more from cooperating with one another.We note that this raises a more general question about temperature.Within subpopulation temperature includes replacement weights between pairs of individuals from different subpopulations, but excludes weights between pairs from within the same subpopulation.What if two individuals have very similar, but not identical, movement distributions (and thus whilst formally not within the same subpopulation, for practical purposes they might as well be)?Under the current definition no distinction is made between this and two individuals whose distributions are completely different.

Figure 1 :
Figure1: The fully independent model from[10].There are N individuals who are distributed over M places such that In visits place Pm with probability pnm.Individuals interact with one another when they meet, for example, I 1 and I 2 can interact with one another when they meet in P 1 .

Figure 2 :
Figure 2: The territorial raider model of [10, 8].(a) Population structure represented using a graph where vertices represent individuals and places.Individual In lives in place Pn and can visit any neighbouring places.For example, the home place of I 1 is place P 1 but it can visit places P 2 , P 3 and P 4 .(b) An alternative visualization on a bi-partite graph where individuals and places are clearly separated.

Figure 3 :
Figure 3: The generalized territorial raider model.(a) Individuals that are members of subpopulation Qm live in place Pm but can visit neighbouring places.The territory of subpopulation {I 1 , I 2 } consists of places P 1 and P 2 , the territory of subpopulation {I 3 , I 4 } consists of places P 1 , P 2 and P 3 , the territory of subpopulation {I 5 } consists of P 2 and P 3 .(b) An alternative visualization as multiplayer interactions on a bi-partite graph where individuals and places are clearly separated.
is a graph represented by a weighted adjacency matrix W = (w ij ) where w ij ∈ [0, ∞) is referred to as the replacement weight.Each vertex v n of the evolutionary graph is occupied by one individual and if w ij > 0 then the individual on v i can place a copy of itself in v j by replacing the individual there.It is assumed that the weights are chosen so that the evolutionary graph is strongly connected, which means that there is a route of finite length between any pair of vertices v i and v j .The weighted adjacency matrix W is therefore said to define the replacement structure.Assuming that there is only one replacement per update event, there are several different ways to calculate the probability of a replacement event r ij where a copy of the individual on v i replaces the individual on v j .In particular, we can broadly classify these in terms of the order in which v i and v j are picked.For birth-death dynamics (BD) the birth event happens first where the individual on v i is chosen for birth with probability b i .The individual on v j is then chosen for death conditional on the individual on v i giving birth with probability d ij , thus we have the replacement probability r ij = b i d ij .For death-birth dynamics (DB) the death event happens first where the individual on v j is chosen for death with probability d i .The individual on v i is then chosen for birth conditional on the death of individual on v j with probability b ij , thus r ij = d i b ij .For link dynamics (L) both birth and death events happen simultaneously and therefore r ij cannot be decomposed.

Figure 4 :
Figure 4: Comparing average fixation probability for different complete structures where figure (a) uses DBD dynamics and figure (b) uses DBB dynamics.Each number indicates a subpopulation of a certain density.For example 60 is a complete structure with 2 subpopulations of size 6 and 0 respectively; 2220 has three subpopulations of size 2 and one of size 0.In each case the public goods game parameters are r = 30, v = 10 and movement parameter is h = 30.We see that in figure (a) for the DBD dynamics, cooperators perform poorly in all cases.In figure (b), cooperators do better for small groups (greater than one).Increasing the number of empty places is beneficial for defectors.

Figure 5 :
Figure 5: Comparing average fixation probability for different δ that is the size (or density) of each subpopulation in a complete graph with 4 subpopulations.The public goods game parameters are set to r = 30, v = 11, the movement parameters are set to h = 30 and dynamics used are DBB.As in Figure4, cooperators evolve better in small groups (larger than 1), namely groups of size two and three, with a small advantage for groups of size four.

Figure 6 :
Figure 6: Figure (a) plots the mean subpopulation temperature against the home fidelity h for a complete population structure with 3 subpopulations of size 2 each.Figure (b) then plots the fixation probabilities against these values of the mean subpopulation temperature where r = 30 and v = 10 for the public goods game, and the dynamics used are DBB.In particular, we notice that the fixation probability of the cooperators is decreasing with the mean subpopulation temperature.

Figure 7 :
Figure 7: Comparing different population structures for the public goods game with various complete graphs for a population size of 12 where (1,12) means there is 1 subpopulation with 12 individuals, (2,6) means there are 2 subpopulations with 6 individuals and so on.We have set r = 30 and v = 10, home fidelity h = 30 and the dynamics used is DBB.

Figure 8 :
Figure 8: Figure (a)shows the effect of compensating for empty places by increasing the home fidelity such that the probability of staying in their home place, pnn, remains the same.We start at h = 30 for the 33 and 222 structures.As an empty place is added, h is increased so that pnn = 30/31 for the 330,. . .,330000 structures and pnn = 30/32 for 2220,. . .,222000 structures.In all cases r = 30 and v = 10.We can see that after compensating for the above effect, the influence of introducing empty places is both reversed and weakened.Figure(b)   shows the mean strict subpopulation temperature dropping off when we compensate for the empty places by increasing the home fidelity such that pnn remains the same.

5. 4 .
The effect of a large home fidelity Consider a well-mixed population of M subpopulations each containing L individuals, so that N = M L, as described in Section 2.1, where h is very large.Consequently from equation (3.1), χ(m, G) is approximately 1

. 2 )
of a defector into a cooperator subpopulation.The terms F L (C) and F L (D) are the fitnesses of cooperator and defector individuals within their own subpopulations, and are obtained directly from equations (2.4) and (2.6), and the terms O(w O ) are of the order of w O , and very small.Further denoting x = v/[r(L−1)]

Table 1 :
Notation used in the paper.