Bipartite Graphs as Models of Population Structures in Evolutionary Multiplayer Games

By combining evolutionary game theory and graph theory, “games on graphs” study the evolutionary dynamics of frequency-dependent selection in population structures modeled as geographical or social networks. Networks are usually represented by means of unipartite graphs, and social interactions by two-person games such as the famous prisoner’s dilemma. Unipartite graphs have also been used for modeling interactions going beyond pairwise interactions. In this paper, we argue that bipartite graphs are a better alternative to unipartite graphs for describing population structures in evolutionary multiplayer games. To illustrate this point, we make use of bipartite graphs to investigate, by means of computer simulations, the evolution of cooperation under the conventional and the distributed N-person prisoner’s dilemma. We show that several implicit assumptions arising from the standard approach based on unipartite graphs (such as the definition of replacement neighborhoods, the intertwining of individual and group diversity, and the large overlap of interaction neighborhoods) can have a large impact on the resulting evolutionary dynamics. Our work provides a clear example of the importance of construction procedures in games on graphs, of the suitability of bigraphs and hypergraphs for computational modeling, and of the importance of concepts from social network analysis such as centrality, centralization and bipartite clustering for the understanding of dynamical processes occurring on networked population structures.


Introduction
Since the pioneering work of Maynard Smith and Price [1], evolutionary game theory [2] has become a valuable tool to describe and study evolutionary dynamics when fitness is frequency-dependent. Evolutionary game theory builds on the theory of games [3] by considering populations of individuals whose success or fitness depends on the outcome of social interactions. Behavioral strategies are genetically or culturally inherited, so that the relative abundance of fitter strategies increases over time due to natural selection or social learning. When populations are assumed to be infinite and well-mixed, the replicator dynamics [4,5] offers a deterministic and exact account of the evolutionary dynamics.
In spite of the importance of the replicator dynamics as a mathematical tool for investigating evolutionary dynamics, it is obvious that real populations are never infinite nor perfectly wellmixed. Games on graphs (see Refs. [6] and [7] for reviews) go beyond these two simplifying assumptions by considering finitesized populations embedded in graphs representing geographical isolation or social networks. A graph G~(V ,E) consists of a set V of vertices and a set E of edges connecting pairs of vertices. In general models of games on graphs, individuals are placed on two graphs with the same set of vertices [8]: the interaction graph G~(V ,E G ) and the replacement graph H~(V ,E H ). Evolution-ary dynamics are specified so that, first, individuals play twoperson games with their neighbors in the interaction graph G, and second, strategy updating takes place along the edges of the replacement graph H. Although the set of edges of the replacement graph may differ from the set of edges of the interaction graph, it is usually assumed that E H~EG so that G and H effectively coincide.
A perusal of the vast literature on games on graphs highlights the importance of network structure in the evolutionary dynamics of different games. Two particular results are worth mentioning. First, although unconditional cooperation under the one-shot prisoner's dilemma (PD) is not evolutionarily stable in infinite and well-mixed populations, it can be viable in sparse homogeneous networks [9][10][11]. ''Spatial reciprocity'' [12], ''network reciprocity'' [13] and ''graph selection'' [14] are different labels that have been coined in order to contrast such effect with other cooperationpromoting mechanisms (see, however, Refs. [11,[15][16][17] for the close connections between network reciprocity and kin selection via limited dispersal, and between games on graphs and inclusive fitness theory). Second, heterogeneous population structures such as scale-free networks [18] can significantly promote cooperation under the PD and other social dilemmas [19,20], although such promotion strongly depends on several details of the network, the payoff functions and the updating rules [7,[21][22][23][24][25][26].
Notwithstanding the importance of pairwise social interactions, many situations in real social systems require the collective action of groups comprised by more than two individuals. Moreover, interactions within these larger groups can not always be represented as disjointed collections of two-person games [27]. Public goods games (PGGs) are paradigmatic among such nondecomposable multiplayer games. PGGs are models of situations where individuals face the dilemma of providing and/or maintaining a public good: a common resource that is both nonexcludable (no individual can be excluded from its consumption) and non-rivalrous (one individual's use of the public good does not diminish its availability to another individual) [28]. Digestive enzymes in yeast [29], ATP in heterotrophic microorganisms [30], webs in social spiders [31], alarm calls in meerkats [32], collective hunting in lions [33], and open-source software in contemporary humans [34] are typical examples of public goods whose abusive exploitation by non-contributing individuals may lead to the socalled tragedy of the commons [35]: a situation in which nobody contributes and therefore no public good is produced or maintained.
By far, the most well known PGG is the N-person prisoner's dilemma (NPD) [36]. In this game, each individual in a group of size N has to decide whether to cooperate (by contributing to a common pot) or to defect (by refraining from contributing). The sum of the individual contributions is multiplied by a factor r and then equally distributed among all players, including those who did not contribute. No matter the decisions taken by the other players, it is always better to defect if 1vrvN. In infinite wellmixed populations and in the absence of cooperation-promoting mechanisms, defection is evolutionarily stable and the replicator dynamics predicts the ultimate extinction of Cs. However, as it is also the case in the two-person PD, cooperation in the NPD can be sustained in structured populations under particular life-cycle assumptions. Hauert et al. [37] studied a spatial NPD resulting from placing the individuals in the nodes of a two-dimensional lattice and restricting interactions to nearest neighbors. In this model and for large values of r (but still for rvN) Cs are able to survive by minimizing interactions with Ds through cluster formation. Hauert et al.'s model has been extended by Santos et al. [38], who used scale-free networks instead of regular lattices as population structures. The highly heterogeneous degree distributions of scale-free networks introduce social diversity both at the individual level (players vary greatly with respect to the number of games they take part in) and at the group level (different games are played by different numbers of players). Social diversity brings up a moderate promotion of cooperation when Cs pay a fixed cost c per game, but a significant boost when Cs pay a fixed cost c for all the games they play. In the following we use the terminology introduced in Ref. [39] and call conventional NPD the former case and distributed NPD the latter.
The way networks are constructed is a common feature of Refs. [37,38] and many other papers dealing with evolutionary multiplayer games on networks [40][41][42][43][44][45][46][47][48][49][50][51][52]. We refer to this construction procedure as the graph approach. According to this framework, nodes of a graph G~(V ,E) define both individuals playing a game and games being played by the focal individual plus its direct neighbors, so that an individual with z neighbors takes part in zz1 games: the one centered on itself plus z games, each centered on one of its neighbors. Fitness or social success is given by the sum of payoffs collected in these zz1 games, and competition or imitation takes place along the edges of the graph.
An alternative way of looking at the population structure resulting from the graph approach is realizing that while the replacement graph is the original graph G, the interaction graph is actually a hypergraph (or a bipartite graph) in which hyperedges (or top vertices) correspond to closed neighborhoods of G (see Figure 1). A hypergraph is the generalization of a graph for the case where edges (called in this case hyperedges) can connect arbitrarily many vertices. A bipartite graph B~(T,\,E), also called a bigraph, consists of two disjoint sets of vertices, T (top vertices) and \ (bottom vertices), and a set of edges, E. The difference between bipartite graphs and standard unipartite graphs is that edges in a bigraph only connect vertices of different kinds. Undirected hypergraphs and bigraphs are mathematically equivalent, but bigraphs are usually easier to implement and to work with. Many real biological and social networks display a natural bipartite structure and can be represented as bigraphs in a straightforward manner. Food webs [53] and metabolic networks [54] are well known biological examples; social examples include affiliation [55] or collaboration networks [56], such as those connecting co-owners of companies [57], film actors [58] and scientists [56]. In this paper, we represent groups/games as top vertices and individuals/players as bottom vertices. Three network statistics will be particularly important. First, the top degree distribution gives the distribution of the number of games being played by a given individual. Second, the bottom degree distribution gives the distribution of the number of individuals playing a given game. Finally, the bipartite clustering coefficient captures correlations between the neighborhoods of bottom vertices, i.e. the degree to which groups overlap.
With the previous definitions, the graph approach can be interpreted as one in which (1) the replacement graph is defined, and (2) the interaction bigraph is constructed from the replacement graph. This approach has become the de facto standard for modeling population structures in multiplayer games on graphs. However, it has many important limitations. First, since both players and games are identified with the same set of vertices, the numbers of games and players are exactly the same, i.e. DTD~D\D in the resulting interaction bigraph. Second, and for the same reason, the top degree distribution and the bottom degree distribution coincide. In real systems, however, these distributions are usually very different. In collaboration networks, for example, the number of papers per author has been shown to follow a power-law distribution while the number of authors per paper generally follows an exponential distribution [56]. Third, the graph approach automatically leads to a relatively large bipartite clustering coefficient. Although such large coefficient seems to be an intrinsic property of many social and biological networks [59], its presence by default in models of games on graphs can be a drawback if the goal is to build null models of connectivity patterns or to study the effects of bipartite clustering. Fourth, while each individual effectively interacts with second-order neighbors in the original graph, strategy updating is posited to occur only between first-order neighbors. As an example of this, consider the graphs depicted in Panels A and B of Figure 1. Note that individual A plays with C and D the game centered on B, but A is not connected to C nor D in the replacement graph. Finally, replacement graphs in the graph approach do not reflect encounter rates between two individuals but rather assume that all neighbors in the replacement graph are equally important. Consider again the graphs depicted in Panels A and B of Figure 1. On the one hand, individual B plays thrice with C (the games centered on B, C and D) but only twice with A (the games centered on A and B). On the other hand, individual B plays (on average) games of smaller size with A (one two-person game centered on A and one four-person game centered on B) and games of larger size with C (two four-person games centered at B and D, and one three-person game centered on C). In any case, the replacement graph posited by the graph approach fails to take into account these heterogeneities, since the connections between B and A and between B and C are equally (un)weighted in this graph.
A new modeling framework for studying networked multiplayer games, recently proposed by Gómez Gardeñ es and co-workers [60,61] and further generalized here, is free of these limitations. We call this framework the bigraph approach. It consists in (1) defining the interaction bigraph B~(T,\,E) so that top vertices correspond to games and bottom vertices to players, and (2) deriving the replacement graph H~(\,E H ) as the (bottom) projection of B (see Panels C and D of Figure 1). The bottom projection of a bigraph B is a graph H~(\,E \ ) so that (i,j)[E \ if and only if i[\ and j[\ are connected at least once to the same top vertex. In addition to the ''unweighted projection'' (UP) considered in Ref. [60], we also consider two weighted projections: the ''unnormalized weighted projection'' (UWP) and the ''normalized weighted projection'' (NWP). With the UWP method, the weight of the link between two players is proportional to the number of common games played by those players; with the NWP method, the sizes of such groups are taken into account when calculating the weights, so that interactions in smaller groups contribute more to the total weight of the link than interactions in larger groups (see the sec:methods section for details).
The bigraph approach circumvents all of the limitations associated with the graph approach we mentioned above. Since the interaction bigraph is defined at the outset, it can have arbitrary numbers of games and players, different degree distributions for games and players and (if required) relatively low bipartite clustering coefficient. In addition, since the replacement graph is obtained as the bottom projection of the interaction graph, individuals playing together at least one game will be connected in the replacement graph. Hence, the neighborhood of player A in the replacement graph shown in Panel D of Figure 1 comprises all the individuals A interacts with, i.e. B,C,D f g . Finally, weighted projections take into account differences in the interaction patterns of players and reflect such differences in the resulting replacement graph. For instance, in the graph shown in Panel D of Figure 1 the weights of the links between players B and A and between players B and C are respectively given by 2 and 3, indicating the number of common games between each pair of players (in Panel D of Figure 1, weights are derived using the UWP method).
In this paper, we make use of the bigraph approach to explore the influence of different topological properties of network structures on the evolutionary dynamics of multiplayer games. We focus on the conventional and the distributed versions of the NPD, as these are among the most studied evolutionary multiplayer games on graphs. Specifically, we investigate the effects of different assumptions on the way of specifying replacement graphs, different top and bottom degree distributions, and different amounts of bipartite clustering. We build interaction bigraphs either from prescribed simple graphs using the graph approach, or from given degree distributions using the configuration model procedure (see the sec:methods section). We denote these interaction bigraphs respectively by the labels fromgraph-X and config-Y-Z, where X stands for the simple graph from which the bigraph is constructed, and Y and Z stand for the degree distributions of the bottom and the top vertices, respectively. Replacement graphs are given either by the graph approach or by the bigraph approach.

Replacement Graphs
In the graph approach, the original graph from which the interaction bigraph is constructed automatically determines the replacement graph. As a result, the subset of players involved in imitation/competition with a given individual is generally smaller than the subset of players with whom such individual interacts. This implicit assumption is in stark contrast with most models of two-person games on graphs, where interaction and replacement neighborhoods perfectly overlap. In the following, we present the results of making interaction and replacement neighborhoods coincide in otherwise standard models of evolutionary multiplayer games on graphs. Figure 2 depicts the results of the evolution of cooperation in the conventional and the distributed NPD for population structures with the same interaction bigraph but different replacement graphs. We plot the cooperation level (the average fraction of Cs for 2000 additional generations after an initial transient of 10 5 generations) as a function of the normalized enhancement factor g~r=n, where n is the average degree of the top nodes of the interaction bigraphs, i.e. n is the mean number of players per game in the population. In each case, the population structure is built (1) by defining a graph G of order Z (i.e. G has Z nodes) and mean degree SzT, and (2) by constructing the interaction bigraph B from In the graph approach the replacement graph H is assumed to be equal to the original graph G, so that interactions take place along the hyperedges of the hypergraph, but strategy updating occurs along the edges of the original graph. The alternative bigraph approach consists in first defining the interaction bigraph B (Panel C) and then obtaining the replacement graph H (Panel D) as the bottom projection of the interaction bigraph. Weights can be attached to the links of the replacement graph according to different heuristics (here, the ''unnormalized weighted projection'' method is used; the width of the links is proportional to the links' weights). The interaction bigraph can be constructed from a bipartite graph model or following the graph approach from a simple graph G (Panel A). In this last case, the replacement graphs due to the graph approach (the original graph shown in Panel A) and to the bigraph approach (the projection of the interaction bigraph shown in Panel D) differ, the latter being denser. doi:10.1371/journal.pone.0044514.g001 G using the graph approach. Hence, n~m~SzTz1, where m is the average degree of the bottom nodes in the interaction bigraphs, i.e. the mean number of games per player in the population. The replacement graph H is given either by G itself (graph approach) or by the projection of B (bigraph approach). In this last case, weights are assigned to the edges of H according to one of three methods: UP, UWP and NWP. In any case, individuals engage in a given number of multiplayer games (according to their connectivity in the interaction bigraph) and accumulate payoffs. The accumulated payoff of each player is then associated with its fitness/success, and competition/imitation is implemented by using a finite population analogue of the replicator dynamics [38,62]: in social learning terms, each individual randomly chooses a neighbor in the replacement graph and, if the neighbor's success is greater than its own success, it imitates the neighbor's strategy with a probability proportional to the success difference (see the sec:methods section for details).
Panel A of Figure 2 shows the results for the case where G is a ring of degree z~4. We refer to the resulting bigraph B as fromgraph-ring. Referring to the original graph G, individuals interact with both their first-order neighbors and their secondorder neighbors. If the replacement graph is given by the graph approach, only first-order neighbors in G are considered for competition/imitation. If the replacement graph is given by the bigraph approach, both first-order and second-order neighbors in G are considered for competition/imitation, possibly with a probability depending on the number of common games (UWP). Note that the larger replacement neighborhoods due to the bigraph approach favors cooperation slightly, but systematically. A detailed analysis of the origin of such promotion, considering the case of two contiguous clusters of Ds and Cs in a ring of degree z~4, can be found in section 1 of Text S1 and Figure S1.
While the larger replacement neighborhoods brought about by the bigraph approach are beneficial to cooperation in bigraphs constructed from rings, they are detrimental to cooperation in bigraphs constructed from Barabási-Albert (BA) scale-free networks, which we call fromgraph-ba. Indeed, as evidenced in Panels B and C of Figure 2, in this case there is systematically less cooperation if replacement neighborhoods coincide with interaction neighborhoods (bigraph approach) than if the original graph is taken as the replacement graph (graph approach). Additionally, in the former case the assignment of weights to the edges of the replacement graph plays a key role in BA networks, as it is evident from the ordering of the curves, with NWP leading to more cooperation than UWP, and UWP to more cooperation than UP.
In order to explain these results, let us briefly recall the mechanism responsible for the promotion of cooperation in the distributed NPD when the interaction and replacement graphs are derived from scale-free networks using the graph approach [38]. Scale-free networks are characterized by the co-existence of few hubs (very well connected individuals) with a vast majority of leaves (poorly connected individuals). Due to their large connectivity, hubs not only take part in many games, consequently accumulating high payoffs, but are also often targeted for competition/imitation by their neighbors. As a result of these two factors, C-hubs and D-hubs easily spread their strategies to their less connected neighbors. However, while C-hubs are favored by a positive feedback mechanism (the more they are imitated, the more Cs in their neighborhoods, and the more their own accumulated payoffs increase) D-hubs are penalized by a negative feedback mechanism (the more they are imitated, the more Ds in his neighborhood, and the more their own accumulated payoffs decrease) that eventually leads to their own demise. Hubs' inherent success along with the feedback mechanisms favoring Cs in interhub competition have been studied using star and double-star graphs as simple models of connectivity patterns in scale-free networks [38,39].
If the replacement graph H is no longer the original graph G (graph approach) but it is rather assumed to be the projection of the interaction bigraph B (bigraph approach), many additional links are present in H that were not in G. Indeed, since each top node of degree z induces a clique consisting of z(z{1)=2 edges, the projection of B is a relatively dense graph, particularly if the top degree distribution is highly heterogeneous [63]. This higher density of the replacement graph is at the origin of the hindering of the evolution of cooperation when moving from the graph approach (G taken as H) to the bigraph approach (the projection of B taken as H). Figures 3 and 4 show this effect for bigraphs built according to the graph approach from star and double-star graphs. In stars, and when the replacement graph is given by the projection of the interaction bigraph, leaves get connected to each other so that H is now a complete graph (see Panel D of Figure 3). This hinders the spreading of cooperative behavior from a Ccenter when defective leaves earn a higher payoff than cooperative leaves. In double-star graphs, leaves of the same star get interconnected and the center of one star gets connected to the leaves of the other star (see Panel D of Figure 4). This increased interconnection hinders cooperation by partially destroying both the positive feedback around C-centers and the negative feedback around D-centers on which inter-hub competition is based in the model of Ref. [38]. Note that, in all cases, the magnitude of these unfavorable effects depends on the weights attached to the links of the replacement graph. Indeed, different projection methods lead to different weight distributions, which in turn affect the topological importance of different nodes in the evolutionary process. Such topological importance can be captured by what we call in this paper replacement centrality, which we define as the expected number of times a given node/individual is selected for competition/imitation by its neighbors (see section 4 of Text S1). Other things being equal, nodes with a higher replacement centrality play a more influential role in the evolutionary dynamics. We find that the level of centralization of the replacement graph (defined as the degree to which a single node is more central than others in the network; see section 4 of Text S1) correlates with the amount of cooperation exhibited in these topologies, as measured by the inverse fixation time of a single Ccenter in a star graph (see Panel E of Figure 3) or by the cooperation level in BA scale-free networks (Panels B and C of Figure 2). Figure S2 shows that the relationship between projection method and centralization of the network found in star graphs (OG w NWP w WP w UP) is maintained in BA scale-free networks. As evidenced by Figure S2 and Panels B and C of Figure 2, weight distributions leading to more centralized replacement graphs are also responsible for higher cooperation levels.

Degree Distributions
In bigraphs constructed using the graph approach, group diversity (heterogeneity in the number of players per game) is inextricably intertwined with individual diversity (heterogeneity in the number of games per player). Indeed, the top degree distribution (determining group diversity) is exactly the same as the bottom degree distribution (determining individual diversity) in bigraphs built using the graph approach. In order to analyze group diversity and individual diversity independently of each other, we made use of random configuration model bigraphs (for which the degree sequences of top and bottom vertices can be specified independently of each other) as interaction bigraphs. We used two different degree sequences for top and bottom vertices: a constant sequence (all degrees are the same) and the degree sequence of a BA scale-free graph, which approximately follows a power-law. Combinations of these two degree sequences resulted in four bigraphs: config-reg-reg (with homogeneous top and bottom degree distributions), config-ba-reg (with heterogeneous bottom and homogeneous top degree distributions), config-reg-ba (with homogeneous bottom and heterogeneous top degree distributions), and config-ba-ba (with heterogeneous bottom and top degree distributions). The reason for using the degree sequence of a BA graph instead of determining the degree sequence by another method (for instance, by sampling the sequence from a random variable distributed according to a power-law distribution) is to be able to compare the results obtained for config-ba-ba with those obtained for fromgraph-ba in the subsec:replacement subsection. Indeed, config-ba-ba has the same top and bottom degree Figure 3. Evolutionary dynamics on stars. Consider a star graph G (Panel A) consisting of one C-center connected to Z{1 leaves (m of which are Cs) and the resulting interaction hypergraph B (Panel B) constructed from G using the graph approach. We assume that social interactions are modeled by the distributed NPD. In the graph approach, G is taken as the replacement graph (Panel C). In this case, competition/imitation occurs only between the center and the leaves. The C-center invades D-leaves for values of r above a critical value which reduces to a~2=(1{2=Z) if m~0 (the C-center is the only C). In the bigraph approach, the replacement graph is given by the projection of the interaction bigraph, so that leaves are now interconnected and the resulting topology is no longer a star but a complete graph (Panel D). The creation of these new links allows for interleaf competition/imitation, which is favorable to Ds if rv4. As a result, for avrv4, the time to fixation to the absorbing state where all individuals are Cs can become arbitrarily large depending on the weights attached to the links of the replacement graph, as it is shown in Panel E for Z~10, r~2:8 and for replacement graphs given by the OG, NWP, UWP and the UP methods. Panel F shows these replacement graphs together with the values of the weights of the links (w hl for the weight of the link between the center and a leaf; w ll for the weight of the link between two leaves) and their centralization indices (r X ). Note that more centralized graphs correspond to those more favorable to the spreading of cooperative behavior from the center. See section 2 of Text S1 for the derivation of the invasion conditions shown in Panels C and D, and sections 4 and 5 of Text S1 for the derivation of the centralization indices shown in Panel F. doi:10.1371/journal.pone.0044514.g003 sequences as fromgraph-ba, and can be effectively thought of as a randomization of such network. Figure 5 shows the results for the evolution of cooperation in the conventional and the distributed NPD for the four configuration model bigraphs and for the fromgraph-ba. Let us consider first the results for config-reg-reg, i.e. the homogeneous population structure lacking social diversity of any kind. As shown in the figure, this network is able to sustain cooperation for values of g above g c &0:7. Furthermore, cooperation is fully established for gwg d , with g d close to its value for infinite well-mixed populations (g d~1 ). For g c vgvg d , Cs and Ds co-exist in dynamical equilibrium. If group diversity is introduced (config-reg-ba), the co-existence zone grows so that g c &0:6, and g d w1:2. This shows that group diversity has mixed effects in the evolutionary dynamics, promoting cooperation (with respect to config-reg-reg) up to a critical value g Ã &0:85, and hindering cooperation above this value. If diversity is instead introduced at the individual level (config-ba-reg), cooperation is evolutionarily viable for gw0:55 in the conventional NPD and for gw0:45 in the distributed NPD. Note, nonetheless, that defective behavior is not completely eradicated, not even for gw1. From these results, it is evident that individual diversity leads to higher cooperation levels than group diversity (compare the curves for config-ba-reg with those for config-reg-ba) for all values of g. We also note that the levels of cooperation slightly improve when both kinds of social diversity are simultaneously present (compare config-ba-ba to config-ba-reg and config-reg-ba). Finally, the results obtained with config-ba-ba are almost the same as those obtained with fromgraph-ba, which suggests that the higher topological correlations present in fromgraph-ba and absent in config-ba-ba play a rather small role in the evolutionary dynamics.
The results for networks with homogeneous bottom degree distributions (config-reg-reg and config-reg-ba) and for networks with heterogeneous bottom degree distributions (config-ba-reg, config-ba-ba and fromgraph-ba-ba) differ not only quantitatively in their cooperation levels, but also qualitatively in their dynamics.  In the graph approach, G is taken as the replacement graph (Panel C). In this case, and for a wide range of values of r, spreading occurs preferentially from the centers (or hubs) to their respective leaves. Long-term evolution will ultimately depend on inter-hub competition, which is favorable to C-hubs due to the positive and negative feedback mechanisms brought about by the spreading from centers to (own) leaves. In the bigraph approach, the replacement graph is given by the projection of the interaction bigraph (Panel D), so that the center of one start gets connected to the leaves of the other star and leaves of the same star get connected with each other. This not only allows successful centers to breed copies of themselves in the leaves of the other star, but also makes inter-leaf competition possible, which is favorable to Ds if rv4. As a result, the feedback mechanisms on which the evolution of cooperation on heterogeneous graphs is based are diminished and the evolutionary outcome is more favorable to Ds. This is illustrated in Panels E and F, which show typical scenarios for the time evolution of the fraction of Cs under the distributed NPD (r~1:3) on the leaves of double-star graphs (Panel E: X~10, Y~20; Panel F: X~20, Y~10), for replacement graphs given by the OG, UP, UWP and NWP methods. In all cases we placed Cs on all nodes of the double-star, except for the left center, where we placed a D (see configurations a and e of Panel G). If the replacement graph is given by the original graph (OG), the dynamics are such that, typically, the D-center invades the leaves of his star (configuration b), then the C-center invades the D-center (c) and finally D-leaves on the left star are invaded by the C-center (d). When the replacement graph is given by the projection of the interaction graph, Ds can now easily spread from the initial center until they invade the whole population (e, f, g). Weights attached to the links of the projection play a key role in this case, with NWP still favoring Cs when the connectivity of the left star is small compared to that of the second star (e, f, h, i). See section 3 of Text S1 for the analytical derivation of the results shown in this figure. doi:10.1371/journal.pone.0044514.g004 almost entirely determined by the proportion of times the dynamics ended up in the full cooperation absorbing state. Figure 6 provides some insight on the different results obtained when diversity is introduced at the individual level (config-ba-reg) or at the group level (config-reg-ba). First, note that the degree distribution of the replacement graph for config-ba-reg is highly heterogeneous (see top panels of Figure 6). Indeed, it is well known that the degree distribution of the projection of a bigraph with a power law bottom degree distribution also follows a power law [64][65][66]. Contrastingly, the degree distribution of the replacement graph for config-ba-reg is less heterogeneous. A second important difference between config-ba-reg and config-reg-ba is the way received benefits are distributed on these networks. When the population consists of 50% Cs randomly placed on the bottom vertices of the bigraph, the distribution of received benefits closely follows a power-law in the case of config-ba-reg, but it approximately follows a normal distribution in the case of config-reg-ba (see middle panels of Figure 6). The reason behind these different distributions is that on config-ba-reg both the percapita per-game contribution and the number of games per individual are highly variable while on config-reg-ba they are constant. This leads to a highly heterogeneous distribution of received benefits on config-reg-ba and to a relatively homogeneous distribution on config-reg-ba. Finally, while there is a strong correlation between connectivity in the replacement graph and received benefit in config-ba-reg, such correlation is practically inexistent in config-reg-ba (see bottom panels of Figure 6). Indeed, for config-ba-reg hubs in the replacement graph are individuals participating in many games and hence accumulating large payoffs. Contrastingly, for config-reg-ba highly connected individuals in the replacement graphs are those participating in large groups, which have on average the same proportion of Cs and hence produce the same amount of public good than smaller groups. As a result, the evolutionary dynamics on config-ba-reg is dominated by a small number of very well connected and powerful individuals, while config-reg-ba is far more homogeneous, both concerning connectivity in the replacement graph and accumulated payoffs. These differences translate into two different modes of evolution.
In config-ba-reg (see Figure S3) the influence of hubs is decisive to the evolutionary outcome, so that a majority of C-hubs leads the whole population to the all-Cs absorbing state, while a majority of D-hubs leads the population to the all-Ds absorbing state. Additionally, the proportion of Cs is also higher in high-degree classes (very well connected individuals) than in low-degree classes (poorly connected individuals). Contrastingly, in config-reg-ba (see Figure S4) the evolutionary dynamics is largely independent of what happens with well-connected individuals, and evolution unfolds as a process of dynamical self-organization in which Cs tend to cluster in small groups which are more favorable to cooperation while Ds tend to do so in large groups which are more favorable to defection.

Bipartite Clustering Coefficient
The bipartite clustering coefficient captures the degree to which bottom vertices' neighborhoods overlap (see section 6 of Text S1 for details). As pointed out in the sec:introduction section, interaction bigraphs built using the graph approach lead, by construction, to relatively high bipartite clustering coefficients. In order to assess the real importance of clustering in the evolutionary dynamics, we considered four interaction bigraphs with the same top and bottom degree distributions (regular sequences in all cases) but different bipartite clustering coefficients: fromgraph-ring (constructed from a ring network of degree z~4), fromgraph-reg (constructed from a random regular network of degree z~4), fromgraph-vn (constructed from a square lattice with a von Neumann neighborhood), and config-reg-reg (random configuration model with regular top and bottom degree sequences). Figure 7 shows the cooperation levels under the conventional NPD and Figure 8 the bipartite clustering coefficient and the mean degree of the replacement graph for these different bigraphs. Interestingly, bigraphs with more bipartite clustering (and consequently lower mean degree in the replacement graph) lead in general to equal or higher cooperation levels for all the considered values of the normalized enhancement parameter g. These results make sense in the light of well established results on the effects of local interactions on the evolutionary dynamics of the pairwise and multiplayer versions of the NPD. It is well known that spatial structure enables Cs to form clusters within which they preferentially interact with other Cs, thus reducing the exploitation by surrounding Ds. Cluster formation is brought about by a feedback mechanism resulting from imitation/competition with direct neighbors that amplifies initial inhomogeneities in the distribution of strategies. As it is shown in Figure S5, large values of bipartite clustering coefficient favor cluster formation by allowing Cs to find each other more easily and to reduce the number of connections with surrounding Ds.

Discussion
Since the seminal works by Axelrod [67] and Nowak and May [9] on the evolution of cooperation on lattices, games on graphs have traditionally made use of unipartite graphs in order to model population structures. Despite its usefulness for exploring the effects of local interactions on the evolutionary dynamics of twoplayer games, the use of unipartite graphs as population structures entails a certain number of construction limitations when applied to general multiplayer games, leading not only to a lack of flexibility but also to unrealistic assumptions about the topological properties of networked populations. In this paper, we have shown how the use of bipartite graphs and of constructing procedures that fully take into account the bipartite nature of social and biological populations can circumvent the limitations of the standard graph approach, opening up new opportunities for studying the role of different properties of network topologies on the evolution of strategic interactions. In particular, it is important to emphasize the need of explicitly defining two graphs: the interaction bigraph, determining who plays with whom, and the replacement graph, determining who competes with whom. As demonstrated in this paper, different ways of constructing any of these two graphs or of deriving one from the other can have important consequences in the evolutionary dynamics of multiplayer games.
First, the implicit assumption that the replacement graph coincides with the original graph in the graph approach is crucial for the success of BA scale-free networks as cooperation-promoting topologies reported in Ref. [38]. When the replacement graph is derived in a more natural way, so that interaction and replacement neighborhoods perfectly overlap (the usual assumption in evolutionary two-person games on networks) cooperation is hindered in BA scale-free networks to a point that any advantage of social heterogeneity is effectively canceled by the resulting larger replacement neighborhoods (see Figure 2). The introduction of weights in the replacement graph somewhat alleviates this problem, as weighted links partly restore the high centralization characteristic of BA scale-free networks.
Second, while individual diversity (heterogeneous bottom degree distributions) systematically fosters cooperation, group diversity (heterogeneous top degree distributions) promotes cooperation up to a critical value of the enhancement factor, but hinders cooperation above such value (see Figure 5). We also showed that networks with both kinds of social diversity foster more cooperation than networks with only one kind of diversity, but that the difference between the cooperation levels of networks with both individual and group diversity and the cooperation levels of networks with only individual diversity are relatively small. Finally, intermediate cooperation levels in networks without individual diversity are mostly due to co-existence of Ds and Cs, while intermediate cooperation levels in networks with individual diversity are characterized by bi-stable evolutionary dynamics. In other words, the results for config-reg-reg and config-reg-ba shown in Figure 5 can be better understood as representing the final proportion of Cs in a population where both Cs and Ds are present. Contrastingly, the results for config-ba-ba, config-ba-reg and fromgraph-ba can be better interpreted as a probability of ending up in a fully cooperative state when starting from a condition where 50% Cs are randomly placed on the network.
Third, bipartite clustering, i.e. group overlap, plays an important role in the evolution of cooperation under the conventional NPD. We provided clear evidence of the beneficial role of bipartite clustering on cluster formation and, consequently, on the evolution of cooperation on regular structures. In this respect, our results mirror similar conclusions on the beneficial effects of unipartite clustering on the evolution of cooperation under the standard evolutionary two-player PD [68][69][70]. Figure 6. Statistics for config-ba-reg and config-reg-ba. The figure shows some statistics for config-ba-reg (left panels) and config-reg-ba (right panels). Top panels: degree distribution of the replacement graph. Middle panels: histograms for the received benefit. The received benefit is calculated as the payoff for Ds when approximately half of the population are Cs (randomly distributed) under the distributed NPD. Bottom panels: smooth scatter plots, regression lines and Pearson's correlation coefficients for the received benefit vs. degree in the replacement graph. Parameters: Z~512, m~n~5 and g~0:7. The figures show statistics for 100 randomly generated networks of each type. doi:10.1371/journal.pone.0044514.g006 Figure 7. Cooperation level for bigraphs with different bipartite clustering coefficients. The interaction bigraphs are constructed following the graph approach with a ring (fromgraphring), a square lattice with von Neumann neighborhoods (fromgraphvn), or a regular random network (fromgraph-reg) of degree z~4 as original graphs, or given by a configuration model bigraph with regular degree sequences for both top and bottom vertices (config-reg-reg). In all four cases the degree distributions of top and bottom vertices is a regular sequence with m~n~5, the replacement graph is given by the normalized weighted projection (NWP) of the interaction bigraph, and Z~1024. doi:10.1371/journal.pone.0044514.g007 Apart from the present paper and to the best of our knowledge, only two studies have made use of the bigraph approach for studying evolutionary multiplayer games: Ref. [60], where the use of bigraphs as population structures for evolutionary games on networks was first introduced, and Ref. [61], a subsequent study on the effects of social diversity on the evolution of cooperation under the NPD. In the first of these studies, the evolution of cooperation under the NPD on a real bipartite collaboration network is compared to the dynamics on its bottom projection. Higher cooperation levels are found for the bipartite network than for its projection. These results have been interpreted as hinting that ''the intrinsic group structure (described by means of the bipartite graph) promotes cooperation in PGGs, this being a new mechanism for this phenomenon beyond the scale-free character and other features of the one-mode (projected) complex network'' [60]. We would like to point out that a simpler explanation is that, by construction, the mean group size in the bigraph built from a projected network is always larger than the mean group size in the original bipartite network, and that larger group sizes hinder the evolution of cooperation under the NPD. In order to assess the influence of group structure and other mesoscopic properties on the evolutionary dynamics, a comparison of real bipartite networks with their ''randomized'' versions should be carried out, as it has been done for real unipartite networks and two-person games [71].
In the second study (Ref. [61]) the evolutionary dynamics of the conventional and distributed versions of the NPD were investigated on interaction bigraphs with tunable individual diversity but no group diversity at all. The main finding of this study is that bigraphs with low individual diversity (Poisson-like bottom degree distributions) can actually allow for more cooperation than bigraphs with high individual diversity (bottom degree distribu-tions following a power law) in the case of the conventional NPD. This result contrasts sharply with our own results, which suggest that individual diversity generally promotes cooperation. Note, however, that we used both a different network model (configuration random networks) and different degree distributions (with zero instead of moderate individual diversity). These different setups could account for the divergent results. We also note that Gómez-Gardeñ es et al. [61] suggest that the ability of BA scalefree networks to outperform homogeneous networks reported in Ref. [38] is ''intrinsically due to the entanglement of social and group heterogeneities''. Although our own results partially support this view, given the (moderate) synergy between individual and group diversity, we have provided evidence that the promotion of cooperation reported in Ref. [38] is mainly due to the implicit assumption that the replacement graph is equal to the original graph from which the interaction topology is constructed.
The choice of the NPD as case of study in this paper was based on the fact that most of the theoretical work on evolutionary multiplayer games has focused on this particular game. However, recent empirical [29] and theoretical [72] work testifies a growing discomfort with the NPD as model of realistic social dilemmas, in particular because of its linearity and because of the fact that cooperation is a strictly dominated strategy in this game. Several of the conclusions drawn in the present study will necessarily change if strategic interactions are modeled after PGGs different from the NPD. For instance, it has been recently shown that, even in the absence of a fixed topology, group diversity can importantly affect the evolutionary dynamics of non-linear PGGs [73]. In the light of these results, we would expect group diversity to play a more prominent role in the evolutionary dynamics of non-linear games played on bigraphs with highly heterogeneous top degree distributions. Also, bipartite clustering could be partially detrimental, instead of largely beneficial, for the evolution of cooperation if the social dilemma is modeled after a multiplayer game with a structure similar to the snowdrift game, as it is already the case for two-person games [62].

Population Structures
Population structures are modeled by means of two graphs: the interaction bigraph B~(T,\,E B ) and the replacement graph H~(\,E H ). The two sets of vertices of the interaction bigraph (T and \) represent, respectively, the set of groups/games and the set of individuals/players.
Graph approach. In what we call the graph approach [37,38] In what we call the bigraph approach [60], first the interaction bigraph B~(T,\,E B ) is defined, then the replacement graph H~(V ,E H ) is constructed by projecting the interaction bigraph into its set of bottom vertices. In addition, weights can be attached or not to the edges of H according to one of the following three methods: 1. Unweighted projection (UP). As done in [60], no weights are attached to the edges or, equivalently, the weights of all edges have a value of one. 2. Unnormalized weighted projection (UWP). The weight w ij of the link (i,j)[E H is given by the number of games i and j are connected to in the interaction bigraph [55]. From a social learning perspective, the reason behind this heuristic is that the more often i interacts with j, the better i is supposed to be acquainted with j and therefore the more often i should consider j as target for imitation. 3. Normalized weighted projection (NWP). The weight w ij is given by [56] w ij~P where d k i~1 if i participates in game k, d k i~0 otherwise, and n k is the number of players of game k. From a social learning perspective, the reason behind this heuristic is the assumption that individuals get acquainted with others more easily in smaller than in larger groups.
Bigraphs built from simple graphs using the graph approach. For fromgraph-X interaction bigraphs, we considered four different kinds of graphs: rings, scale-free networks, square lattices with von Neumann neighborhoods and regular random networks. Rings are one-dimensional lattices with degree z. Regular random networks (maximally random graphs where each node has the same degree z) were constructed using the igraph [74] implementation of the algorithm by Viger and Latapy [75]. Scale-free networks were obtained by means of the Barabási-Albert (BA) model [18], i.e. growing networks using a preferential attachment rule. In order to get graphs with average degrees exactly equal to SzT, we started the growing procedure from a fully connected graph of m 0~S zTz1 nodes, and added m~SzT=2 new edges per new node.
Configuration model bigraphs. Config-X-Y bigraphs were constructed using the configuration model [63,64,76] with a top degree distribution of type X and a bottom degree distribution of type Y. For the degree distributions, we used regular sequences (reg) and degree sequences from BA scale-free networks (ba), constructed following the procedure mentioned before.

Multiplayer Games
Each individual i participates in all games k such that (i,k)[E B . The social success of an individual is given by the sum of the payoffs obtained in all games it takes part in. We considered two versions of the NPD: the conventional NPD and the distributed NPD [38,39]. In the conventional NPD, the payoffs of a D and a C in a group k of size N k are respectively given by P D~r m k c=N k and P C~PD {c, where m k is the number of Cs in group k, c is the cost of cooperation and r is the enhancement factor. In the distributed NPD, each C of degree z i (i.e. taking part in z i games) contributes c=z i to each game, so that the overall contribution of any C is equal to c. In this case, the payoff of individual i with strategy s i (1 if C, 0 if D) is given by [38]

Evolutionary Dynamics
The success/fitness of each individual was calculated as the sum of the payoffs obtained in all the games it participates in. Strategies are updated synchronously using a finite population analogue of the replicator dynamics commonly used in the literature of games on networks [38,62]. When updating the strategy of individual i, a neighbor j of i in the replacement graph is randomly chosen with a probability p ij given by where w ij is the weight of the link (i,j)[E H . Denote by P i the accumulated payoff of individual i. Then, if P i §P j , i stays with its current strategy; otherwise it changes its strategy to j's with a probability given by (P j {P i )=M, where M is a normalization factor given by the highest possible difference between the accumulated payoffs of i and j.

Simulations
Simulations were started with 50% of Cs randomly placed on the graph. We measured the average fraction of Cs for 2000 additional generations after an initial transient of  Figure S1 Evolutionary dynamics on rings. In the inset, we plot a ring of degree z~4: the neighborhood of each node comprises the closest two nodes to the left and to the right. Following the graph approach, each node is the center of a game of size five so that each individual ends up interacting with the closest four neighbors to the left and the closest four neighbors to the right. We assume that the initial distribution of strategies is such that nodes iw0 are Cs and nodes iƒ0 are Ds. In the main panel, we plot the probabilities of switching strategies for the individuals at the boundary (nodes 0 and 1) when the replacement graph is given by the original graph (OG) and when it is given by the unweighted projection (UP) of the interaction bigraph. As shown, P(s 0 ?C)wP(s 1 ?D)urw5=7 for the graph approach, while P(s 0 ?C)wP(s 1 ?D)urw1=2 for the bigraph approach. See section 1 of Text S1 for the calculation of these probabilities. (TIFF) Figure S2 Centralization of the replacement graphs for interaction bigraphs built from Barabá si-Albert scalefree networks. Each boxplot shows the distribution of the centralization for a random sample of 10 4 replacement graphs given by the original graph (OG), the normalized weighted projection (NWP), the unnormalized weighted projection (UWP) and the unweighted projection (UP). In all cases, the original graph is a Barabási-Albert scale-free network of order Z~512 and mean degree SzT~4. The projections are taken from bipartite graphs constructed from the original graph using the graph approach. Notice that more centralized networks lead to higher cooperation levels in Panels B and C of Figure 2 in the main text. See section 4 of Text S1 for the definition of the centralization indices used in this figure.

Supporting Information
(TIFF) Figure S3 Time-dependence of the fraction of cooperators for different connectivity classes in the config-bareg network. The figure shows the fraction of Cs among lowdegree (z i vm), medium-degree (mƒz i vz max =3) and high-degree (z max =3ƒz i ƒz max ) individuals, for two different simulation runs. In Panel A, initially more than the 60% of the highly-connected individuals are Cs. C-hubs lead the evolutionary process and diffuse cooperative behavior among their less connected neigh-bors. In Panel B, initially less than 40% of the hubs are Cs. Less connected individuals quickly turn to defection, with mediumdegree and high-degree individuals eventually following the trend. Parameters: g~0:7, m~n~5 and Z~512. (TIFF) Figure S4 Time-dependence of the average experienced group size and of the fraction of cooperators in groups of different size for config-reg-ba. The figure shows the mean experienced group size for Cs and Ds (top panels) and the fraction of Cs in small (N i vn), medium-sized (nƒN i vN max =3) and large (N max =3ƒN i ƒN max ) groups (bottom panels) for g~0:7 (left panels) and g~1:2 (right panels). The evolutionary dynamics on this population structure is such that Cs preferentially cluster together in small groups and Ds cluster together in large groups. Parameters: m~n~5 and Z~512. Text S1 Supporting Text on analysis of the evolutionary dynamics on rings, stars and double-star graphs, and on definitions of replacement centrality, centralization, bipartite clustering coefficient and degree of assortment. (PDF)