Modelling Evolution in Structured Populations Involving Multiplayer Interactions

We consider models of evolution in structured populations involving multiplayer games. Whilst also discussing other models, we focus on the modelling framework developed by Broom and Rychtář (J Theor Biol 302:70–80, 2012) onwards. This includes key progress so far, the main gaps and limitations, the relationship and synergies with other models and a discussion of the direction of future work. In this regard as well as discussing existing work, there is some new research on the applicability and robustness of current models with respect to using them to model real populations. This is an important potential advance, as previously all of the work has been entirely theoretical. In particular, the most complex models will have many parameters, and we concentrate on considering simpler versions with a small number of parameters which still possess the key features which would make them applicable. We find that these models are generally robust, in particular issues that can arise related to small payoff changes at critical values and removal of pivotal vertices would have similar effects on other modelling system including evolutionary graph theory. These often occur where it can be argued that there is a lack of robustness in the real system that the model faithfully picks up, and so is not a problematic feature.


Introduction
Evolutionary game theory is an important methodology in the modelling of biological populations. The original evolutionary game theory models involved well-mixed populations, generally centred upon contests involving pairs of individuals (see [56,57]). Sometimes, there were no such contests at all, and games were "against the field", i.e. involving the whole population [39]. Simple models, for example the sex ratio game [39] and the Hawk-Dove game [57], have been used to explain initially surprising biological phenomena. The most commonly used evolutionary dynamics for well-mixed models are the replicator dynamics [41,93]. Even for very simple models of this type, there is still a lot to discover, see, for example, [9]. For an overview of much of classical game theory and also more modern applications, see [14,20,30].
Real populations, however, have structure with some individuals more likely to encounter others based upon factors such as geographic location, social connections and specific behaviours. Perc et al. [72] is a good review of modelling involving structured populations and evolutionary dynamics, focusing on multiplayer games and in particular public goods games. As explained in [72], how to model structure often reduces to two key questions: (1) How are pairs or groups of individuals formed to determine the payoffs that contribute to the fitnesses of individuals? (2) Given the fitnesses, how are individuals selected for replacement or reproduction? Thus, a number of ways have been developed to incorporate structure into the population.
A natural way to incorporate structure, especially spatial structure, is through regular lattices, e.g. hexagonal (with 6 neighbours per individual) or square (with 4 or 8 neighbours, depending upon whether corner connections are counted). Here, every individual occupies a point on a lattice and interacts with those close by. This can be in pairwise games against each of its neighbours or in a multiplayer game against all of its neighbours. In the latter case, an individual will also be involved in multiplayer games centred upon each of its neighbours, and hence, it will interact with all of its neighbours' neighbours. The importance of group interactions involving individuals that are not directly connected was shown in [91]. That increasing group size did not necessarily tend towards mean field behaviour was shown in [90] (see also [89]), and in general the relationship between the evolution of cooperation, for example, and group size is not necessarily straightforward [42].
In some models, all individuals update simultaneously, with changes depending upon their type and configuration of neighbours [76,100]. This can lead to interesting, but often not biologically realistic, patterns, and in keeping with the focus of this paper, here we consider asynchronous updating only. Often a single individual is selected at random, its fitness calculated by the appropriate games, and then, a random neighbour is selected and its fitness similarly found. The first individual then replaces the second with a direct copy of itself as some increasing function of the difference of their fitnesses [72].
More recently, evolutionary graph theory was developed, which in many ways is similar to the above; for example, each vertex of a graph is still occupied by a single individual, but now a much greater variety of potential structures can be considered. This originated with [49] (see also [5,16,52,54,67]) with good reviews (of work up to the time of publication) in [2,80]. Both the structure of the population and the evolutionary dynamics used can have big effects on evolution, (see [16,49,97,98]). For example, population structure is one of the key factors that can help cooperative behaviour evolve [60]. Heterogeneity can be an important factor allowing clusters of cooperators to form around hubs [48]. Theoretical work for general graphs is hard except for in the case of weak selection, where important simplifications can occur. A powerful recent methodology based upon coalescents has been developed in [3].
Real networks are often highly heterogeneous and follow certain structural patterns. For instance, the degree distribution can often follow a power law structure [29]. There are a number of ways of generating large structures with interesting properties randomly [77], and often real complex networks, for example [51,75,88], which consider networks of dolphins, macaques and zebra, respectively, are used as the basis of models to evolve on in [3].
In some models, the individuals are not fixed on a particular place but can move over the underlying structure. These movements can be random (see, e.g. [28]), or can be affected by their interactions with other individuals, for example [43,104]. They may choose to move with a higher probability for lower payoffs from interactions, or if they have an unfavourable group composition. We can see how a version of such models fits into our framework in Sect. 4.
In other models, the structure itself can evolve. This might be because links between individuals are the focus of strategic interactions, as in the models of [84][85][86] or [11][12][13]27]. In the former, the structure changes as new individuals are born and die, and interesting and complex patterns emerge out of simple rules, reminiscent of the cellular automata discussed earlier. In the latter, the population is fixed but the links change, in a constant battle between individuals trying to achieve a favourable configuration. Other models are co-evolutionary [64,65], and as the composition of the population changes, the structure changes too. For example, if links can be adjusted based upon payoffs received as in [102] or other related factors [103,105,106], this can strongly promote cooperation.
One limitation of evolutionary graph theory (without evolving structure) is that interactions are usually pairwise, although as described above in the case of lattices, multiplayer games can be incorporated by having an individual play a game with all of its neighbours simultaneously (e.g. [47]). The group sizes here are still constant for regular graphs, and there is a constant set of groups for each individual on an irregular graph, and so the process still lacks flexibility and realism. The methodology that we principally describe in this paper was developed by [17] to model multiplayer games with a variable number of players within general structures. We shall go on to describe this methodology in detail.
In Sect. 2, we introduce the framework and describe how population structure is introduced, including the simplest independent model. In Sect. 3, we introduce the evolutionary dynamics and consider the full dynamic model. In Sects. 4 and 5, we consider models more complex than the simplest where specific independence assumptions are removed (history independence in Sect. 4 and individual movement independence, "row independence", in Sect. 5). In Sect. 6, we consider the potential applicability of our framework and its various models to real populations, and in Sect. 7, we consider the robustness of our models to errors, for example in measurement. Finally, Sect. 8 is an overall discussion of work so far.

The Framework of Broom and Rychtář: Population Structure
Broom and Rychtář [17] introduced a new framework for modelling the evolution of a population in a structured environment. This framework was based upon three core components: (1) the population (environment) structure; (2) the interactions between the individuals (evolutionary games); and (3) the evolutionary dynamics of the population. The framework can accommodate a variety of structures and evolutionary games; moreover, the evolutionary games are multiplayer games, and they may include interactions in groups of varying size, including lone individuals. The initial paper [17] defined the modelling framework and discussed several example structures, but it did not include the evolutionary dynamics. In this section, we shall introduce the population structure and (multiplayer) interactions among individuals in a structured environment and compare several sample structures.

Basic Setup
Broom and Rychtář [17] considered a population of N individuals I 1 , . . . , I N who can move to M distinct places P 1 , . . . , P M . The distribution of the individuals over the places at time t was defined by a binary matrix X (t) = (X n,m (t)) of size N × M given by if the individual I n is at the place P m at time t, 0, otherwise, and the probability distribution of the population at time t, dependent on the entire history of the system x <t = (x 1 , x 2 , . . . , x t−1 ), was written as where x and x 1 , x 2 , . . . , x t−1 are any particular values that the matrix X may take. Similarly, denoted the probability of the individual I n being at the place P m at time t given the history x <t . The population distribution thus followed a random process, which could depend upon its entire history. This allowed a vast range of possible models and therefore [17] defined several specific scenarios, each of which included some simplifications.
If the population distribution is independent of the history of the process: the model was called history independent. A history-independent process was called homo- If the probability of an individual visiting a place does not depend upon the simultaneous movements of other individuals (although it may depend upon the history): then the model was called row independent.
If the model is history independent, homogeneous, and row independent, then we have the following independent movement probabilities: p n,m,t (x <t ) = p n,m for all n, m, t, x <t , and in this case, the model was called simply independent. Figure 1 provides a sample illustration of an independent process. Broom and Rychtář [17] introduced one special case of an independent model, the territorial raider model. In this model, individuals I 1 , . . . , I N each occupy their own place P 1 , . . . , P N , and the places are arranged as vertices of a connected graph. Individuals may either stay at their own place, or move to one of the neighbouring places (although returning to their own place for the start of the next move). Each individual follows its own distribution. The movement patterns within this model may still be complex, and [23] introduced a single Fig. 1 Illustration of an independent process. The node I i is connected to the node P j only if the probability p i, j of the individual I i visiting the place P j is not 0 Illustration of a territorial raider model. The territory is divided into four places P 1 , . . . , P 4 . Each place P i is home to an individual I i . The individual I i may move freely within the place P i and any of its neighbouring places. This model can also be described by a graph where the vertices are the places and the neighbouring places are connected by the edges. With the home fidelity parameter h, the movement probabilities are given by h+2 . The territorial raider model is an independent process, and it represents a particular case of the scenario in Fig. 1 home fidelity parameter h, which governed movement in the entire population. The home fidelity parameter could take any nonnegative value, and any individual was h times as likely to stay at its own place than to move to any of the neighbouring places. Figure 2 provides a sample illustration of a territorial raider model. Broom and Rychtář [19] considered a variant of an independent model involving strategic movements, the territorial movement game. As in the standard setup, the game was played by N individuals I 1 , . . . , I N that could move between M distinct places P 1 , . . . , P M . A collection of all places P m that could be visited by an individual I n was called the territory of the individual. Each individual chose a place in its territory to move to. These choices were made simultaneously (i.e. not knowing what others chose to do). Once an individual made its choice and visited the place, it received a payoff based on the total number of occupants in that particular place. The fitness of any individual was simply given by this payoff, and the solution to the territorial movement game was such a movement that any individual would have received a lower payoff if it unilaterally changed its position (given the fixed positions of other individuals). This solution was called the stable occupancy function.
There were two main classes of the territorial movement game considered: (1) a habitat selection game, which is a generalisation of the classical ideal free distribution [31,35]; and (2) a predator dilution game. In the habitat selection game, the payoffs to individuals decreased with the number of the occupants at their place. It was shown that there was a unique stable occupancy function over the places for the unstructured population (where all individuals could visit all places), but that this was not true in general for structured populations. In the predator dilution game, the payoffs to individuals increased with the number of the occupants at their place. In general, there could be many stable occupancy functions for this game, including in the unstructured case. A general method for finding all such solutions was given.

Fitness
An important component of any evolutionary model is the fitness of the individuals within the population. The reward for an individual I n at time t was denoted R(n, x, t, x <t ). Assuming only the current distribution affects the reward, the mean reward can be considered as the fitness where R(n, x) is the reward to I n given the population distribution x.
The reward to an individual at any given time depends upon the others that it is interacting with and the place where this interaction occurs. The group G of individuals meeting at P m is G = {I j ; x j,m = 1}. We note that in the territorial movement game model of [19] described in the previous section, this fitness was greatly simplified by the fact that it only depended upon the number of individuals at a place (movements were strategic, but no games were played at the place) and all individuals adopting pure strategies led to a particular group distribution with probability 1 for each strategy combination.
Let P(X •,m = χ G )(x <t ) be the probability of a group G meeting at the place P m at time t given the history x <t . For the independent model, this probability can be computed as Often the reward to an individual will only depend upon the place that it occupies and the group of individuals at that place; in this case, R(n, x) = R(n, m, χ G ).

Multiplayer Games
We present three multiplayer games that were considered in [17] and later papers. Each game has two strategies, which we label A and B, respectively.
The multiplayer Hawk-Dove game Here, A represents the Hawk strategy and B the Dove strategy. All individuals that meet at a place play a multiplayer Hawk-Dove game, where they compete for a single reward of value V . If all individuals are Doves, they share the reward. If there are one or more Hawks, the Doves flee and the Hawks fight, with the winner getting the reward and the rest paying a cost C. All individuals also receive a background reward of value R, representing the reward gained from other activities. Thus in a group of a Hawks and b Doves, the average payoffs for the Hawks and the Doves, respectively, are The public goods game Here, A represents Cooperate and B represents Defect. Each cooperator pays a cost C to provide a benefit V , which is shared by all other individuals in the group (but not the providing cooperator); note that the cost is paid even by lone cooperators. All players also receive a background payoff R. Thus, for a group of a Cooperators and b Defectors, the average payoffs to Cooperators and Defectors, respectively, are We should also note that the public goods game presented above is just one version of a social dilemma, and there are many other types of games that can involve multiple (and indeed a variable number of) players, which pick cooperative versus non-cooperative strategies. These have been treated in a variety of papers (e.g. see [7]), and considering the variability of the number of players in particular was an important focus of [71] and [24].
The fixed fitness game Here, individuals receive a constant payoff which only depends on their own type (i.e. there is no real interaction and this does not depend upon the group composition), with type A individuals receiving reward V and the background reward R, and type B only receiving R. For a group of a type As and b type Bs, we thus have

Comparison Between Structures
The payoffs to individuals engaged in multiplayer interactions usually depend on the size of the groups where such interactions take place. The group sizes, in turn, depend on the underlying structure. For example, in the territorial raider model on a complete graph any group size from 1 to N is possible, but on a circle graph only groups of size 1, 2, and 3 may form. Consequently, the population structure crucially influences the payoffs in the multiplayer game. To see whether the population structure can influence the outcome of a multiplayer game in other ways besides possible group sizes, Broom and Rychtář [17] defined the notion of a completely mixed population where p n,m = p for all n and m. Every individual is equally likely to meet any other individual in a completely mixed population, and the mean group size in such a population is 1 + (N − 1) p. If a population structure has a mean group size g, one can compare it to a completely mixed population with parameter p = (g − 1)/(N − 1). This makes the comparison fair as both cases have the same mean group size (but may differ in higher moments).
As an example, it was shown in [17] that for the territorial raider model with the Hawk-Dove game, evolution would lead to a lower proportion of Hawks on the star graph than for a completely mixed population. This analysis was then extended to other graph structures in [26]. For example, we should expect a higher proportion of Hawks on the circle graph than for the completely mixed population.

The Framework of Broom and Rychtář: The Full Dynamic Model
The evolutionary dynamics of finite (structured) populations are based on a Moran process to keep the population size constant, and six main updating mechanisms are used in evolutionary graph theory [54,67,68]. Each evolutionary dynamics is based on a series of updating events. Each updating event consists of one individual being selected for birth (reproduction) and one individual for death (replacement); the offspring of the individual selected for birth replaces the individual selected for death. (In classical evolutionary graph theory, an individual selected for birth may replace only an individual to whom it is connected by an edge.) The six evolutionary dynamics differ in the order in which the birth and death events take place and the stage at which selection occurs. Several dynamics may also be combined into a single mixed dynamics, see, for example, [94,107]. Different updating mechanisms can yield different evolutionary outcomes, see, for example, [101].
One of the most commonly studied evolutionary dynamics is the BDB dynamics, also known as an invasion process [6,49,59,60,83]. In the BDB dynamics, first an individual is selected for birth, and then, it copies itself into a randomly selected neighbour; individuals' fitness is used to weigh the birth process. Broom et al. [23] adapted this dynamics into the framework of [17] by considering different structures for population interactions and population dynamics; such a decoupling of these structures was previously considered in [61,62]. In these earlier works, the population and replacement structures were completely unrelated. But in our framework, the replacement structure is determined by the population structure and composition. The probability that an individual I i replaces an individual I j follows a distribution that depends upon the frequency of interactions between these individuals. Specifically, the probability of I j being replaced by a copy of is the probability of the individual I i being selected for birth (weighted by the fitness), and d i j = w i j /( k w ik ) is the probability of the individual I j being selected for death given that I i was selected for birth. Here, w i j are the replacement weights, which can be thought of as the expected proportion of time the individual I i spends interacting with the individual I j . The replacement weights are computed as where χ(m, G) is the probability of a group G meeting at the place P m . In particular, an individual selected for reproduction may replace itself with its own offspring; the population composition does not change in this case. We have k w ik = 1, and hence for the BDB updating mechanism, where fitness plays no role in the death event, In [67,68], all six evolutionary dynamics were considered, and they are summarised in Table 1.
With at most two types of individuals (A or B) present in the population at any time, the state of the population is uniquely determined by listing the positions of individuals of one type. Let S denote the set of positions (graph vertices) occupied by type A individuals in the population and P SS denote the transition probability from the state S to the state S in the Table 1 Summary of six evolutionary dynamics used in evolutionary graph theory. Here, F j is the fitness on an individual I j and w i j are the replacement weights. The fitness derives from a multiplayer game, and the values of w i j come from (7) Dynamics dynamic process of the game. Then, for the BDB updating mechanism, we have when S = S and Because the population is finite, the evolutionary dynamics ultimately results in fixation of either type A or type B individuals in the population. The fixation probability of A individuals occupying a set S, A S , is a solution to the equation with the boundary conditions One of the key results of [67] was that the temperature T n = i =n w in , 1 ≤ n ≤ N , was a good predictor of the average fixation probability of a single type A individual (mutant) in a population of N − 1 type B individuals (residents). High (low) temperatures generally enhanced (suppressed) the strength of selection. This observation agrees with other results of evolutionary graph theory, see, for example, [21, 22,58] In particular, finding these temperatures provides a computationally efficient method of estimating the mutant fixation probabilities [79], and it may potentially lead to some general analytical insight. The mean temperature was also a good predictor of the fixation probabilities for a territorial raider model in [23]. But this paper considered only populations on small networks (size 4), and hence, the general nature of these predictions remained an open question. Larger populations and a wider variety of network structures were considered in [77]. Whilst the general thesis of the earlier investigations was confirmed, this more detailed study showed that the relationship between the mean temperature and fixation probability is complex, and in particular, it depends on the game being played.
Additionally, Pattni et al. [68] considered an extension of the territorial raider model where each place was home to a subpopulation of individuals rather than a lone individual. The DBB and DBD dynamics were used for the evolution of this subdivided population playing the public goods game. This led to a natural extension of the notion of temperature, the subpopulation temperature, which was a better predictor than the standard temperature for this type of population. The conclusion was that for the DBB dynamics low subpopulation temperatures were conducive to the evolution of cooperation in the public goods game, whilst cooperation could never evolve for the DBD dynamics. Moreover, small subpopulations (but not subpopulations of single individuals) were particularly helpful for cooperation to evolve.

Breaking History Independence: The Markov Model
Pattni et al. [69] extended the territorial raider model to the case where individuals explored the environment following a Markov process. Each vertex of a graph was home to exactly one individual. Individuals explored the environment through a fixed number of movement steps, which comprised an exploration phase. At each movement step, every individual decided strategically whether to stay at the current location or to move to one of the neighbouring locations at random. These decisions were made independently of the choices of the other individuals. Individuals that ended up in the same location after any given movement step played one round of the public goods game, which was described in Sect. 2.3. The multiplayer game payoffs accumulated over the entire exploration phase were added up to inform the individuals' fitness. At the end of the exploration phase, all individuals instantaneously returned home to reproduce. After the reproduction event (using the BDB dynamics), which might have resulted in a different population composition, the fitness of all individuals was reset, and a new exploration phase began.
We will now explain how the individuals explored the environment and the Markov nature of the process. Each individual in the population had two independent traits that constituted its type: (1) the interactive strategy (cooperate or defect), which determined the individual's strategy in the public goods game, and (2) the exploration strategy (staying propensity), which corresponded to the probability that an individual was going to stay at the current location if the individual was alone. The actual decision whether an individual was going to stay or move depended on its intrinsic exploration strategy and the current state of the environment; whence the Markov property.
Given a current state of the environment during an exploration phase, each individual evaluated the composition of the group it was currently part of. If the individual found the group attractive, then it was more likely to stay, and if it found the group unattractive, then it was more likely to move to another location. More specifically, let G n (m t−1 ) denote the set of individuals that are located at the same place as the individual I n at time t − 1; we call G n (m t−1 ) the group of the individual I n at time t − 1. The attractiveness of the group to the individual I n (who is a member of the group) was computed as where β j is the attractiveness of a group member I j . This attractiveness depended on the interactive strategy of the individual, and it was defined as and the values β C = 1, β D = −1 were used. Letting α n be the staying propensity of the individual I n , the probability h n (G n (m t−1 )) that an individual I n stayed at time t − 1 at its current location, where it was part of the group G n (m t−1 ), was computed using the sigmoid function where 0 < σ < 1 was the sensitivity to the group composition parameter. The smaller the value of σ , the more sensitive individuals are to the group composition, that is, the more likely they are to react to the composition of their group. In particular, it follows from (17) that h n (G n (m t−1 )) = α n if the individual I n is alone. Notice that the staying propensity α n in this model is similar to the home fidelity parameter h in the basic territorial raider model presented in Sect. 2.1. However, in this model the staying propensity was one of the traits of the individuals, rather than a parameter of the model, and it could evolve. Each movement incurred a fixed cost λ, which was subtracted from the payoff in the corresponding round of the public goods game [given by (3) and (4)]. Consequently, higher movement costs resulted in higher optimal staying propensities in the population.
Pattni et al. [69] demonstrated that lower movement cost, longer exploration time, and larger population size helped cooperation to evolve and resist replacement by defectors on the complete graph. One notable observation was that the type of evolutionary dynamics used to evolve the population did not have much effect on the outcome. In particular, cooperation could evolve under the BDB updating mechanism, which is not the case in classical evolutionary graph theory [60].
This work was extended in [33] to larger populations and different network structures: the complete, circle, and star graphs. The patterns observed in [69] for the complete graph held for larger population sizes, but the overall outcomes varied among different network structures. It turned out that the stability of the population of defectors was determined mainly by the movement cost: mutant cooperators could replace resident defectors only if the movement cost λ was less than a threshold value, which varied only slightly between different networks and population sizes. But the stability of the population of cooperators was determined mainly by the network structure. Cooperators resisted replacement by mutant defectors for all movement costs on the complete graph as long as the population was sufficiently large, but mutant defectors always replaced cooperators on the star graph. On the circle graph, there was a threshold intermediate value of the movement cost so that cooperators dominated for all lower values, and defectors dominated for all higher values. [33] conjectured that two of the network topology characteristics-the clustering coefficient and degree centralisationwere partially responsible for these outcomes, and a road map for further investigation of this phenomenon was proposed.

Breaking Row Independence: The Correlated Movement Model
Broom et al. [25] considered row dependence and developed a general method for generating correlated movement in a history-independent process. In particular a suite of methodologies for generating animal dispersal patterns was developed. The primary aim of this paper was to model animal movement in terms of specific mechanisms and to consider the relationships of such mechanisms with the resulting placement distributions. We can predict the outcome if the dispersal mechanism is known, and conversely we can test what mechanisms could result in a given outcome.
The secondary aim of the paper was to be able to generate distributions with certain desired properties, such as the probability distribution over the distinct places that could be visited, and a given level of coordination among the animals. We developed a range of measures of coordinated behaviour, so that we could create models with sufficient flexibility, and investigated the properties of all of the models that we developed.
The paper considered a range of models where individuals arrived sequentially, and chose their patch according to predefined rules. A certain utility was usually assigned to each patch, which was a function of both the intrinsic value of the patch and the number of individuals on it. Thus, denoting the value of a patch P m by V m and the current number of individuals on the patch by Y m , the utility of the patch was given by U m (V m , Y m ).
One class of model was purely deterministic, where individuals selected the patch with the currently highest utility. An example of a deterministic model is the ideal free distribution discussed in the previous section. For example, following Parker's matching principle [66], we have Here, the more individuals are on a patch, the lower its utility is, and hence, the individuals will spread over the patches proportionally to the intrinsic patch values. There could also be situations where the utility of a patch increases with the number of occupants, and then, the individuals will tend to form a single herd of the entire population. Utility functions could also be combined to form other alternatives (e.g. a weighted average of a number of functions). On the other hand, individuals do not need to choose the place with the best utility with probability 1. We considered a second class of models where the places were chosen randomly weighted by the utility. For a utility function independent of Y m , this yielded an i.i.d. distribution, which was the baseline of comparison for other models. Broom et al. [25] also considered different urn models in this context, for example the Polya urn [44]. Furthermore, we devised a procedure, called the wheel, that could generate movement distributions with prescribed probabilities to visit specific places and with varying degrees of coordination between the individuals. Figure 3 illustrates possible outcomes in two probabilistic models: (1) the random distribution model, where the utility U m of a patch is equal to its intrinsic value V m , and (2) the Polya urn model with 10 balls. We also created hybrid models, which combined basic models to create models with intermediate properties. There were two classes of hybrid models: (1) the whole population used the same model, which was selected at random, and (2) each individual independently chose one of the available models according to a given probability distribution.
We considered measures of aggregation within populations. In particular, the measure T was defined as the probability that a randomly selected pair of individuals is on the same patch. Denoting the final number of individuals on the patch P m by X m , the value of this measure is The theoretical maximum of T is 1, and we found an expression for the lower bound on T , conditional on a given distribution of a randomly chosen individual over the patches. We also showed that all possible intermediate cases could occur, giving us a very flexible methodology.

Framework Applicability
We have considered a number of models above. To what extent are they applicable to real populations? To be applicable, they need to be complex enough to explain real features, but also sufficiently simple so that parameters can be estimated to prevent model overfitting. Analytical models often employ too few parameters, whilst models based on simulations usually employ too many parameters because it is tempting to include anything that might be relevant.
In each of the model types above, we can identify the number of degrees of freedom necessary to fully describe the model, which corresponds to the maximum number of parameters that are potentially needed to be estimated in order to use the model. But this number is often too large for practical purposes. It is therefore important to identify simpler versions of the models with fewer parameters. For instance, the territorial raider model may contain up to M − 1 movement parameters for each individual (to determine the probability p m of moving to place P m , for m = 1, . . . , M − 1; the last probability being given by 1 − M−1 m=1 p m ), but the notion of home fidelity reduces this number to a single parameter for the whole population as long as the underlying population structure is fixed. We will indeed be assuming below that the population structure (network) is fixed. In general, the population structure also contributes model parameters. A graph with M vertices (places in the environment) can have up to M(M − 1)/2 edges and hence the population structure could be thought to have that many parameters, although we could alternatively list all possible graphs of size M and then use the position of the graph in the list as a single structure parameter.

Independent Model
Here, N individuals may move to a maximum of M places. Each individual goes to exactly one place, and hence, the movement of each individual is governed by a probability distribution with M components, resulting in M − 1 free parameters. Thus, there are N (M − 1) degrees of freedom. A well-mixed population considered in [17] had only one parameter. For the territorial raider model, we have N = M and hence there are N (N − 1) degrees of freedom. [23] developed a version of the territorial raider model with a single parameter, the home fidelity parameter (as described in Sect. 2.1), which works for any structure.

Individual Markov Movement Model
Here, each individual moves conditional upon its own most recent position only. It has a probability distribution over M places, with M − 1 free parameters as above, but this distribution can be different for each potential place it currently occupies. Thus, we have N M(M − 1) degrees of freedom.
The simplest model should include at least two parameters, for example, a baseline staying propensity parameter [69] and an adjustment parameter for the current position of an individual. An adjustment parameter could be realized as a multiplier, which depends on the degree of the vertex an individual currently occupies. As an example, if an individual with staying propensity α is located at a vertex of degree k, then it stays at its current place with probability αβ k /(αβ k + k) and moves to any neighbouring place with probability 1/(αβ k + k).

General Markov Model
Here, the distribution of each individual could depend upon the specific combination of the current places of all N individuals. There are M N such combinations, so we have N M N (M − 1) degrees of freedom.
The model in [69], discussed in Sect. 4, contained five movement parameters: one governing the cost of movement, one governing the exploration time before a reset of the population, and three governing the preference for staying or moving to a neighbouring place at each time point. In this case, the role of the adjustment parameter β was played by the attractiveness of the individual group members (β for cooperators and −β for defectors), but another parameter was sensitivity to the group composition.

History-Independent Correlated Movement Model
Here, we consider the placement of N individuals over M places as a single outcome. If each individual is unique, we would have M N possible outcomes (i.e. M N distributions of individuals over the places). However, if the individuals are identical, then we are interested in how many (rather than which) individuals occupy a given place, that is, we are interested in the number of M-tuples (n 1 , n 2 , . . . , n M ) with n m ≥ 0 and m n m = N . This is a variant of a stars and bars problem [34], and it is well known that there are N +M−1 N of such Mtuples (see also [20,Exercise 9.1]). One of the possible outcomes must happen, hence we have N +M−1 N − 1 degrees of freedom, and we can choose many parameters to describe the probability of a given distribution of N identical individuals over M places.
There are a number of different ways to construct models with fewer parameters. One natural way is to consider M − 1 parameters e 1 , e 2  Another natural way to construct such models is to specify the quality of each place (i.e. M parameters in total); this is done, for example, in the ideal free distribution model [31,35] or, more generally, in habitat selection games [19,25].
In each of the above models for arbitrarily different places, all M − 1 parameters e 1 , e 2 , . . . , e M−1 are needed, but these (or some underlying place quality) values might be specified in advance, and so they are not free model parameters (analogously to specifying an underlying graph structure). For example, we might want to investigate the case where all places have the same expected number of individuals.
Of course, even the simpler models previously described can have a relatively large number of degrees of freedom if more sophisticated modelling is required. For example, in general urn models [44] we may need to describe whether or not a ball will be returned to the urn after the draw (since it may depend on the ball number, there are M degrees of freedom) and also how many balls (and with what numbers) will be added to the urn after the draw. (This means we need to specify up to M numbers for each of the M possible draws for M 2 degrees of freedom.) For the wheel, the coordination between individuals may be more complex. The relative positions of the individuals can be given by random distributions, and then, for each of the individuals I 2 , . . . , I N , we would have to specify the distribution and all its parameters resulting in a potentially unbounded number of degrees of freedom. However, one often considers only a single and simple random distribution, such as a normal distribution, and we thus need only 2(N − 1) parameters to describe the mean and standard deviation of the position of each individual if each individual has its own distribution; this number may be considerably smaller otherwise. For the ideal free distribution or the habitat selection game in general, one should potentially specify the payoff to an individual at any of the M places when in the group of size 1, 2, . . . , N , resulting in up to M N degrees of freedom.
Combining pairs of distinct deterministic models through either of the hybrid models requires only a single parameter to give the weighting of each deterministic model used, which then varies smoothly between each of the models in the extreme. However, we should note that different models will require different numbers of parameters, so the combined model will require the sum of the numbers of parameters of each model plus one.
In general, we need to carefully think about the situation that we wish to model and see what number and type of parameters are necessary to build a sufficiently realistic model.

Markov Correlated Movement Model
If we had a Markov correlated movement model, then the whole population would move over the possible N +M−1 N combinations, conditional upon the current combination. Thus, we would have N +M−1 N N +M−1 N − 1 parameters. The minimum parameter model for a specified target would require two parameters: one that governed herding as discussed in Sect. 6.4 (e.g. B from the Polya urn) and another that governed the history-dependent aspects as discussed in Sect. 6.2 (e.g. staying propensity α).

Game and Dynamics Parameters
We assume that every player in the game can follow one of S strategies and that the payoff to an individual depends only on how many individuals are together and what those individuals play (i.e. not on who exactly is with the individual and not on what happens on other places). The payoff will depend on the place (M possibilities), the strategy of the focal player (S possibilities), and how many of the remaining N − 1 players are with the focal player and what strategies they choose. To evaluate the latter, notice that we have to assign N − 1 players to one of the following S + 1 roles: be with the focal player and play one of the S strategies, or be somewhere else (in which case their strategy does not matter). This gives N −1+S S possibilities. Overall, there are SM N −1+S S − 1 degrees of freedom ("−1" since some combination must occur). We note that as S is unbounded, so is the potential number of degrees of freedom.
If we consider games with two strategies only, which are commonly used in structured evolutionary models, there are then 2M N +1 2 − 1 = M N(N + 1) − 1 degrees of freedom. Of course in practice, much simpler games with relatively few parameters are used. Here, we focus only on the three multiplayer games that we considered in Sect. 2.3. The background reward R and the benefit V are two parameters in all of the three games (Hawk-Dove, public goods, and fixed fitness games). Moreover, the cost C was an additional parameter for the first two games. However, often it is only the relative size of the parameters that matter and so one can fix V = 1 and specify only the other two (or one).
We note that the dynamics are governed completely by the values of w i j and F i , which in turn depend upon just the structure parameters (in the case of w i j ) or the structure and the game parameters (for F i ). Thus, once we have selected our dynamics, there are no additional contributions not already involved. We note that in the games, one of the parameters is the background fitness parameter, so this could be replaced with a strength of selection parameter within the dynamics, but this will leave the number of parameters involved unchanged.
The above discussion assumes that one specific dynamics was chosen, similarly to fixing the underlying structure as discussed at the beginning of Sect. 6. But there are several different dynamics, and they could be mixed in some ways. These considerations may result in more parameters that govern dynamics selection.

Framework Robustness
Modelling a biological system requires estimating the parameters of the model. Any potential errors in such estimates may affect the outcomes of the model. We therefore investigate the robustness of our framework with respect to small changes in the parameter values. We start with simpler systems related to our framework, from which it is partially derived, and explore if and when small changes in parameters can result in different outcomes. We then discuss the meaning of our results in general terms, and the implications both for modelling and for a more general interpretation of biological systems.
For evolutionary games in finite populations, small changes in parameter values for game payoffs will mostly have a correspondingly small effect on the fixation probability. However, sometimes they can have a huge influence on it. For a large population with fixed fitnesses (1 for residents and r for mutants), if r changes from slightly larger than 1 to slightly lower than 1, then the fixation probability can change from an appreciable value to near 0 if the population starts with a number of mutants. In the usual starting situation of a single mutant, the fixation probability for r > 1 and large N is just over 1 − 1/r (and exactly this in the infinite population limit), so that for instance a fitness reduction from 1.1 to 1.05 would reduce the fixation probability almost in half (from 0.091 to 0.047), and a further reduction of the same size would take it to (just over) 0. For games in a finite population behaviour is more complicated, but similar threshold cases arise.
The above is just a special case of games on graphs in general, and hence small changes in fitness can yield significantly different outcomes for games on graphs too. However, in a structured environment, we also need to estimate the underlying population structure. An example of a small change in the population structure is a graph rewiring, where an edge is removed and replaced with a different one. In general, such simple rewirings do not have a significant effect on the fixation probability, as demonstrated in [22]. Thus, from this point of view it seems that measuring the payoff function accurately might be more important than getting the underlying population structure exactly right.
Adlam and Nowak [1] found a robustness result for the isothermal theorem related to random graphs and essentially showed that rewiring has little effect in this case. Hindersin and Traulsen [40] studied the effect of edge removal on the fixation time, which was robust to the removal even of a large number of edges. Möller et al. [58] carried out a numerical study which also supported the robustness of evolutionary graph properties to structural change. We note that there will be some situations where changes to structure can have a big effect, for example specific graphs like the star that are subject to vertex or edge removal, as this will make the graph unconnected. Similarly, there are some very special directed graphs which act as amplifies (or suppressors) of selection [49] where the structure is critical and small structural changes would make big differences. Amplifying graphs and their properties are discussed in detail in [4,94].
In general, vertex removal (especially targeted rather than random removal) for graphs with hubs can have a large effect. However, it seems that changes due to errors in observation of real structures which are otherwise well connected would not lead to large changes.
The most basic model in the framework of [17] is the independent model and particularly the territorial raider model. Even in this simplest model there could be issues with payoffs near certain critical points. Aside from these exceptional cases, Broom et al. [23] showed that the fixation probability was fairly robust to changes in the parameter values. Even changes between graph types seemed to make little difference for the small graphs considered in that paper. A much larger class of graphs was covered in [77], and in particular, they considered some important graph properties in assessing the fixation probability. The temperature and mean group size seemed good predictors in general, but there could be big differences in the fixation probabilities of graphs with the same temperature, especially for certain types of graph. This was most pronounced for the more extreme cases of the Barabasi-Albert graphs [8], which are scale-free networks constructed using the mechanism of preferential attachment, and preferential attachment leads to high variability in the degrees of the vertices. Thus on occasion, knowledge of a small number of graph properties (at least those we are currently aware of) is not always sufficient to predict the fixation probability, and this is true for evolutionary graph theory too. To this end, [58] investigate all small graphs considering fixation probabilities and fixation times. Rychtář and Stadler [74] showed that gradually changing the underlying graph from completely regular (a square lattice) to a completely random one increases the fixation probability with the randomness. Nevertheless, the increase in fixation probability is not significant, and thus, the fixation probability could be estimated by the known formulae for regular graphs.
A version of the territorial raider model was developed in [36,37]. Here, individuals live on the graph vertices, each vertex having a value 1. Each individual decided whether to stay at home, or raid a neighbouring place. There were some advantages in contests for the defender. The game resulted in strict Nash equilibria if and only if all individuals raided someone else's place. This is possible for some graphs and not others, and one can convert a graph from one type to the "opposite" type via a simple rewiring. Thus in this case, the solution can be broken by a small change in the graph structure. This is an example of non-robustness. We note this is somewhat different from finding the fixation probability.
Small changes-smaller than a change in the edges involved-can often radically change Nash equilibrium (NE) results. NEs are based upon deterministic calculations where all other individuals behave rationally, typically without stochasticity. NEs on games in sequential form can be sensitive to the choices that would be made on parts of the game tree that are never visited, for example, and there is much subtlety (and many NE concepts) involved in their analysis, see [96], and for other interesting complications involving sequential strategies, see [18].
For the Markov movement models discussed in Sect. 4, changes in parameter values (except for the same critical cases as described above) in general did not make significant differences in [69], which considered the complete graph for relatively small populations. Different graphs were considered in [33], and here, there was one case where the graph structure led to significant changes in outcome for small parameter changes. For "large" populations on the circle graph, a small change in movement cost (from 0.1 to 0.2) leads to the optimal choice of the staying propensity for defector mutants jumping from the smallest possible value (0.01) to the largest possible one (0.99). However, this did not lead to a nontrivial change in the fixation probability. Also, this jump phenomenon in the best response can occur even in simple matrix games (and even two-strategy matrix games), where a small parameter change can switch the best strategy from choosing the first strategy of two with probability 1 to probability 0.
We have hypothesised that the behaviour of the population in this model depends upon topological characteristics of the underlying networks such as the clustering coefficient and degree centralisation. The three networks we have considered-complete, circle, and star graphs-have extreme values of these characteristics: 1 (0), 0 (0) and 0 (1), respectively, for the clustering coefficient (degree centralisation). We observed completely different outcomes for the stability of the population of cooperators on these three networks. We will test this hypothesis for other networks with intermediate values of the clustering coefficient and degree centralisation in future work, which is in its early stages.
We have not yet incorporated the correlated movement models from Sect. 5 into evolutionary processes. But we can already see some interesting cases of non-robustness. We considered three distinct mechanisms whereby individuals formed a single herd in a deterministic optimisation model. These mechanisms were "follow the leader" (the first individual is denoted the leader and all others follow this individual), "follow the majority", and "follow the predecessor". In a hybrid type I model, combined with some other distribution, various types of aggregation could be achieved, but all three of these would give the same result. However, for hybrid type II they were all different. In particular hybridising with the Random model, with a small probability of picking Random, is the same as the selected deterministic model with rare random errors. Even small errors multiply leaving all three cases very different. We note though that this does not mean that the process per se is not robust, only that it is possible to generate a non-robust model if sufficient care is not taken in properly defining the mechanism.
In many cases where a model is not robust, it may not be that the model is a poor representation of reality. Here, the reality may not be "robust" in some sense too. Thus small changes in the real world can lead to huge effects. A classic example of this is the "butterfly effect", where small errors lead to large long-term ones, see, e.g. [32,50]. Similar chaotic behaviour emerges even from the simple logistic map [55], which relates to the population dynamics of fruit flies. The development of animal colouring patterns such as spots and stripes can emerge through sensitive dynamical changes, so-called Turing patterns [53,78,95].
Real situations are impossible to model accurately at such critical junctures, as the smallest error in measurement would lead to a completely wrong prediction, just as the smallest change in the real situation would. Thus, we should not necessarily be so concerned when our models have problems of this kind, which may simply be faithfully reflecting reality. In general, a feature of a good model is that if a model characteristic is (not) robust and predicts some real phenomena of our interest, then that phenomena should (not) be robust in reality.
We have found situations where our framework is robust to small changes, and situations where it is not. The non-robust cases linked to the structure are generally due to some extreme properties of the structure, such as the Barabasi-Albert graphs from [77], which may not often be relevant to real populations. Other features, to do with payoffs, can be critical in general, but these are at certain threshold cases where similar phenomena occur for even the simplest evolutionary models.

Discussion
In this paper, we have looked in detail at the modelling of evolution in structured populations, with a focus on the modelling framework developed by Broom and Rychtář [17]. This has involved a review of previous work and also new work considering the robustness and applicability of the different models.
There are key elements in common with other models of evolution in structured populations. The key modelling systems, the most prominent being evolutionary graph theory [16,49,60], involve a structure, types of interaction (typically games) and an evolutionary dynamics for how the population changes. Whilst these are often interlinked, with weightings used in the evolutionary dynamics arising from the structure on which the games are played, this does not have to be so. The evolutionary dynamics can be completely separated arising from a totally different underlying structure, as in [61,62]. This is realistic if individuals spend much of their lives away from their breeding sites, for example. The framework from [17] and evolutionary graph theory can both easily involve such uncoupled models where relevant, though work using the same structure for both interactions and dynamics is much more common. The coalescent-based methodology of [3] in particular may potentially be applicable in many situations. One distinction of the methodology in [17] is that multiplayer interactions of varying sizes arise naturally, rather than having only pairwise contests, or other imposed fixed group sizes.
We note that an alternative method of incorporating such structures was developed in [92], where there is a range of available sets to be members of, and every individual in the population is a member of its own collection of sets. Here, variable group sizes occur and can change throughout the evolutionary process, often due to individual strategic choices of updating their set membership following payoffs received. The work in [92] focused on a situation with cooperating and defecting individuals, and cooperators only cooperated with those who shared sufficiently many common sets. The methodology is potentially much more general, and there is considerable scope for further exploration.
Another way of building upon network structure to make it more complex is by extension to a multilayer network. A good review of this modelling approach is given in [99]. Here, the standard graph structure is extended to contain many levels, each of which has its own set of graph properties, with connections between the levels also able to vary accordingly. They classify them into three main types: multiplex networks, interdependent networks, and interconnected networks. The latter is a way of more efficiently analysing what at first sight might appear a complex graph, and this is related to the idea of sub-populations in the [17] framework, as introduced in [68]. Unlike for interconnected networks, interdependent networks have separate network structures, but the individuals within them are connected by some non-network property, for example if those in one provide (or compete for) resources for (with) another. Multiplex networks on the other hand contain all vertices in each layer, examples being metabolic networks or transport networks covering different systems between the same cities. The situation where the evolutionary dynamics occurs in a different relationship to the interaction networks as described above is also an example of a multiplex network [62]. Wang et al. [99] then go on to consider evolutionary games on multilayer networks.
In Sect. 2.3, we discussed a simple multi-player public goods game which has been used as the basis of investigating the evolution of cooperation on the various structures involved. In general, the evolution of cooperation is very heavily studied, and there are many different types of such multiplayer models [24]. These models are extended further to consider many different aspects of behaviour involved in cooperation, including punishment, pooled punishment etc. These are considered in great detail in [73], who then go on to consider the evolution of such processes on a structure, principally a rectangular lattice. Important insights into key features of populations, for example the conditions for cooperation to evolve, can be found using weak selection versions of the models, as we discussed in the introduction. Whilst solving more complex models in the strong selection case is difficult, and coming up with general analytical solutions is only possible for very special cases, they can be tackled in other ways. The most obvious way is by direct simulation, but there are also approximation methods that can lead to more general insights. Such methods for approximating the evolutionary process were developed in [63], using existing methods from epidemic modelling as a basis. There are in fact a number of such methods from statistical physics used for modelling epidemics [10,45,46,70,81,82], and there is certainly more potential to apply such methods to evolutionary processes in the future. See also [38], which applied similar methods to the modelling of kleptoparasitism [15,87].
In general, we believe that the methodology we have discussed in this paper is a promising way to model complex populations where multiplayer interactions of various types and with varying group sizes take place. There are also other interesting methods, as discussed above. Such work is certainly still in its early stages, with much opportunity for future exploration.