## Abstract

Starting with Darwin, biologists have asked how populations evolve from a low fitness state that is evolutionarily stable to a high fitness state that is not. Specifically of interest is the emergence of cooperation and multicellularity where the fitness of individuals often appears in conflict with that of the population. Theories of social evolution and evolutionary game theory have produced a number of fruitful results employing two-state two-body frameworks. In this study, we depart from this tradition and instead consider a multi-player, multi-state evolutionary game, in which the fitness of an agent is determined by its relationship to an arbitrary number of other agents. We show that populations organize themselves in one of four distinct phases of interdependence depending on one parameter, selection strength. Some of these phases involve the formation of specialized large-scale structures. We then describe how the evolution of independence can be manipulated through various external perturbations.

## 1. Introduction

Cooperative behaviour, as exemplified by multicellular life, seems to have evolved at least 25 times independently—once for plants, once or twice for animals, once for brown algae, and possibly several times for fungi, slime moulds and red algae [1]. On shorter time scales, the social composition of eukaryotes such as *Saccharomyces cerevisiae*, and biofilm-forming bacteria such as *Pseudomonas aeruginosa* can dramatically change in a brief period [2–5]. In a related context, tumour formation is a rare example of the transition, taking place in the reverse direction, from a multicellular to an essentially unicellular lifestyle. Interestingly, cancer cells end up cooperating by collectively secreting angiogenic factors, and it seems possible, at least in principle, that there may even be cheaters (i.e. those who do not secrete the growth factors) among this collection of cooperating cheaters [6,7].

Evolutionary game theory provides excellent insight into how altruistic and cooperative behaviour can emerge to maximize the fitness of the group despite the apparent fitness advantage of cheating individuals [8,9]. In the context of evolution of cooperation, these models typically investigate the outcome of repeated runs of the prisoner's dilemma between pairs of agents that have two strategies, cheating and cooperating. Variants of the model include structured interactions, coupled populations, coevolution, stored reputation, punishment and preferential or random partner choosing [10–18].

However, real life is more complicated in a number of ways. First, many actual games are massively multi-player [19–22]. The fitness of an organism may depend on its simultaneous relationship with multiple players. Second, biology allows for a much larger variety of internal states beyond cooperating or defecting. For example, the genetic make-up of an organism may be suitable for cooperation with only an exclusive few, while some organisms may be incapable of defecting or cooperating all together. Third, real social evolution leads to highly organized dependence structures beyond the homogeneous mixtures or aggregates of cooperator–defector states that we often see in the majority of evolutionary game models. From biochemical to societal scales, life organizes itself in highly complex arrangements of cliques, communities, cycles and hierarchies.

Without compromising the simplicity and tractability offered by traditional evolutionary game theory, here we propose an evolutionary model in which the fitness of an agent is determined, not by the outcome of a two-player two-state game, but instead a multi-player multi-state one. Thus, our focus is not the cooperator–defector ratio in the population, but rather the large-scale structure of all exchanges; i.e. the *interdependence* between agents or groups of agents within a genetically heterogeneous population. We ask how independent agents become interdependent through the simple laws of evolution, whether positive selection is a necessary or a sufficient condition for the formation of interdependence, what kinds of interdependent structures are stable/unstable, and how these structures and processes depend on evolutionary parameters.

Accordingly, the present model offers a clear framework for classifying and categorizing different regimes of interdependence, as well as allowing for careful control of evolutionary parameters that may be influencing recent non-intuitive empirical outcomes [23]. We determine which kinds of external perturbations promote anti-sociality (e.g. in order to eradicate biofilms) and which other kinds can inhibit anti-sociality (e.g. as to suppress or reverse tumor growth) by simulating the introduction of selfish/altruistic strains into a population or the administration of anti-sociality/sociality promoting drugs. We evaluate the success rate of these evolutionary interventions as a function of the original population structure, drug dose, fraction of drug-resistant agents and reproduction speed of the target species.

## 2. Model

Our multi-state multi-player game can be best visualized as a network of *N* ≫ 1 agents, connected by directional edges. An edge from *A* to *B* indicates that *A* contributes to the fitness of *B* at the cost of its own. Unlike the typical evolutionary game theoretic models where the state of a player *i* is binary (cooperator/cheater), here the player states *ψ _{i}* are characterized by high-dimensional vectors, i.e.

*ψ*

_{i}= {

*x*

_{1},

*x*

_{2},…,

*x*

_{N}} with

*x*

_{j}∈ {0, 1} indicating whether

*i*provides a fitness benefit to

*j*. The evolutionary dynamics is governed by the following assumptions (figure 1):

(1) The fitness

*ω*(*x*) of a node_{i}*i*is assumed to be a monotonically increasing function of received net benefit*x*_{i}=*bn*_{i,in}−*cn*_{i,out}, where*n*_{i}_{,in}and*n*_{i}_{,out}are the number of in and out edges for node*i*. The parameter*β*=*b*/*c*quantifies the benefit of an edge (to the receiver) relative to its cost (to the provider).(2) Every generation, the

*r*most fit nodes produce offspring that replace the*r*least-fit nodes. Reproduction preserves all edge relationships of the parent, i.e. parents and offspring connect to the same agents.(3) There is a small mutation probability

*p*per generation with which edges are added/removed randomly.

In other words, we have a fitness-based selection rule keeping the number of agents *N* constant. Our simulation code is available as an electronic supplement.

Our model has four parameters, kept constant throughout the course of evolution: population size *N*, mutation probability *p*, number selected for replacement *r*, and the relative benefit *β*. For every run, we keep track of the total number of edges *E*(*t*) as a function of generation number *t. E*(*t*) is a measure of the interdependence of the population as well as the average fitness (the latter follows from which can be positive or negative depending on the value of *β*). In addition, we study the community structures and genetic composition within the population, which are defined in terms of the connectivity matrix *C* of the network. We use the convention that *C _{ij}* = 1 if

*j*depends on

*i*, and 0 otherwise, and represent these by black and white pixels in array plots. Our simulations were run for

*N*= 200, partly due to computational constraints. Although this might appear small, we note that the relevant degree of freedom here is the number of edge slots

*N*

^{2}= 4 × 10

^{4}, ensuring that the evolutionary transitions we report are not accidental fluctuations.

## 3. Results

The evolutionary dynamics and final interdependence states depend on the values of *β* and relative selection pressure *r*/*m*, where *m* = *N*^{2}*p* is the expected number of mutations per generation (which is equal to the number of mutants if *p* ≪ 1). We can divide the parameter space into four regimes: neutral constructive (*r* ≪ *m*, *β* > 1), selective constructive (), neutral destructive (*r* ≪ *m*, *β* < 1) and selective destructive (). In all cases, the formation of an edge is beneficial for one of the nodes and deleterious for the other. However, the value of *β* determines the change in average fitness per edge, d〈*ω*〉/d*E* = *c*(*β* − 1)/*N*, and therefore one intuitively expects *E*(*t*) to decrease for *β* < 1 and increase for *β* > 1 as long as selection strength is finite. We will see that this expectation will not always be satisfied, particularly when *β* > 1.

Let us start with the straightforward case of *β* < 1. Here the formation of an edge is more deleterious to its originator than beneficial to its target, and the fitness of the population changes by *b* − *c* < 0 per edge. Thus, in the long run, only if the selection is weak (*r* ≪ *m*) can such deleterious edges accumulate, and we get a random interdependence network with *E* = *N*^{2}/2. As expected, increasing *r*/*m* causes the network to become sparse and fragmented, and all structure vanishes as .

We now overview the constructive regime *β* > 1, which produces distinct phases of complexity (figure 2*a–f*). The long-term behaviour of *E*(*t*), which can be viewed as a proxy for average fitness as well as interdependence and complexity, depends non-monotonically on selective pressure *r*/*m*: for small values of *r*/*m*, the asymptotic value *E*(*t* → ∞) decreases with *r*/*m* (figure 2*a–c*). However, if *r*/*m* exceeds a critical point we see sudden transitions between well-defined discrete levels (figure 2*d–f*). As selective pressure is increased further, we see increasingly larger fluctuations around these levels.

We describe these asymptotic states in more detail by connectivity matrices (figure 3*a*) and phylogenetic trees (figure 3*b*) for varying levels of *r*/*m*. The phylogenetic trees are obtained by quantifying the similarity distance between all pairs of nodes *i* and *j*. In other words, if *i* and *j* receive from and provide to the same nodes, they are considered to be genetically related, consistent with our reproduction rule (cf. figure 1).

While the destructive *β* < 1 regime produces either random or sparse networks, the constructive regime *β* > 1 can be summarized in terms of a sequence of complex phases governed by the value of *r*/*m* (figure 4): a transition from cooperation to competition between individuals (figure 2*a–c*) is followed by unstable interactions between individuals and ‘bunches’ (figure 2*d*), followed by a transition from cooperation to competition between ‘bunches’ (figure 2*e,f*). We define a *bunch* to be the opposite of a graph theoretical community; a group of nodes that form denser connections towards other groups, than they do within (cf. last two panels in figure 3). Dense outwards connections and sparse intra-connections are the key qualitative characteristic of a highly specialized system. For example, nearly all energy spent by a heart muscle cell is directed at serving other tissues. The same holds true in a specialized society, e.g. a lawyer dedicates most of her effort defending non-lawyers. The interdependence structures we report in the last two panels of figure 3, conform to these biological and social examples of specialization.

We now move towards an understanding of the *control and manipulation* of the evolution of interdependence, which is now experimentally possible (albeit with mixed success) in biomedical and ecological settings. For example, the sociality of *P. aeruginosa* can be manipulated by drugs that suppress the microbe's production of a common good (iron scavenging siderophores). As the microbes that are resistant to the drug will altruistically continue to produce the expensive siderophores, they are taken over by their selfish counterparts affected by the drug [5,24,25]. As a result, the iron-deficient population can be easily annihilated by the host's immune system [26]. Note that the evolutionary fate of the drug-resistant group would have been the opposite, had the drug been an antibiotic instead of a quorum blocker. On the other hand, there have also been experiments yielding the *exact opposite* outcome, where the drug aggravates the infection instead of impairing it, presumably by levelling the relative advantage of cheaters [23]. We will use our model to quantify these mixed outcomes. Social evolution is complex, and its manipulation and control requires a detailed quantitative understanding of the evolutionary outcomes of varying initial states and system parameters.

To manipulate the sociality of a highly interdependent, *β* > 1 population, we start with an initial network that has a given community structure, select a fraction *η* of the population and block a fraction *γ* of their outgoing connections of those that are selected. Following this perturbation, we track the evolution of the network and check if our perturbation causes the entire population to lose all connections (which, for *β* > 1, amounts to minimal fitness). If *E*(*t*) drops to and remains at zero we count it as a success and we determine the fraction of successes for every parameter value. If *η* ≪ 1, the perturbation can be interpreted as an external introduction of a new strain/species, or a novel mutation which introduces a very small number of selfish individuals in the population. If *η* ≃ 1, the perturbation can be interpreted as a drug or event that inflicts nearly everyone, such as the quorum blocker discussed earlier. Accordingly, the quantity *γ* can be interpreted as the *dose* of the drug, or the degree of ‘selfishness’ of the newly introduced species/strain.

Figure 5 displays the dependence of success as a function of *η* (empty versus closed plot markers correspond to *η* = 2% and 98%), initial population structure *k* (quantifying the number of bunches) and *γ*. We consider fast reproducing and slow reproducing species separately, shown in figure 5 (*a*) and (*b*), respectively.

It is important to distinguish between two very different mechanisms that bring the population back to its pre-perturbed state. The first is determined by the time required for interdependence to evolve anew from *E* = 0. The original factors causing the establishment of cooperation in the first place is present regardless of our perturbation, and the effect of even the strongest drug (*η* = 1, *γ* = 1) is to simply reset the evolutionary clock. The second mechanism is evolution through the repopulation of the drug-resistant fraction, which happens much faster, on reproductive time scales. To clearly distinguish between these two mechanisms, we set *p* = 0 in figure 5; using a finite *p* scales down all the success rates but does not otherwise change the qualitative dependence on *η*, *γ* and *k*.

We observe a number of interesting features in the response of the population to external perturbations. For *η* ≃ 1, we find a non-monotonic dependence of success rate to *β* for populations with few bunches: for slow reproducing populations a moderate dose works as well as, or better than, a strong one. For larger numbers of bunches, and faster reproduction rates the non-monotonicity vanishes: the stronger the dose, the better the outcome. A second remarkable outcome is the degree to which a few individuals can make a difference: targeting *η* = 2% of the population is as effective as targeting *η* = 98% of the population provided the drug has a high enough dose. This is because few selfish individuals, as is the case in tumours or invasive species, can devastate an entire population. Finally, we observe a very strong dependence of the success rate on the initial community structure. With increasing *k* and *r* this difference vanishes.

## 4. Discussion

The phase diagram for the evolution of interdependence is shown in figure 4 and compactly summarizes our results. Figure 4*a* shows the number of dependences, while Figure 4*b* quantifies their structure through ‘bunch modularity’. We define the latter by exchanging 1 ↔ 0 in the connectivity matrix and determining the community modularity. We will label regimes of interdependence with A–D and discuss them in the following subsections. Regimes A/B refer to cooperation/competition between individuals, while regimes C/D refer to cooperation/competition between bunches.

### 4.1. Cooperation between individuals

In the neutral regime (*r* ≪ *m*), additions and deletions of edges are equally likely. Thus, if the population starts fully independent, *E*(*t*) increases until the network is fully randomized with *N*^{2}/2 edges. This increase is statistically irreversible and is the analogue of the scenario described in [27–29]. In the neutral regime, individuals have high fitness due to the benefits of indirect reciprocity; interdependence emerges not due to higher fitness but simply due to higher likelihood. Figure 1*a* shows the dynamics and final outcome of nearly neutral evolution. Despite following a similar trajectory to fully neutral evolution ( indicated by the dashed curve in figure 2*a–f*), the connectivity matrix shows the onset of community formation; the interdependence structure and genetic composition of the population is far from random (cf. figure 3*a*,*b*).

As the selection strength is increased (*r* < *m* but not *r* ≪ *m*), the fluctuations in *E*(*t*) are amplified. This is caused by the random formation of nodes for which the number of in-edges are different than out-edges. However, the system is self-stabilizing (figure 2*a,b*); e.g. when the fit defectors reproduce, they typically replace their unfit providers, which in turn reduces their own fitness. Consequently, they are taken over by the fair and fit nodes that dominate the *r* ≪ *m* population. In figure 3, third panel, two large reciprocating groups can be distinguished. They are taken advantage by smaller scale opportunistic sub-populations. It is also possible to see a smaller sub-cooperative group sustaining itself within a larger cooperative group.

### 4.2. Competition between individuals

As *r* approaches *r*_{c} ∼ *m*/2 from below, the fluctuations in *E*(*t*) become comparable to *E*(*t*) itself. Here the selective competition is just high enough to allow for small cooperative communities to form and grow at a rate much higher than random chance, but also high enough for cheaters to spread over their providers in one step, beyond recovery (figure 1*c*). Although regime A and B have similar destabilizing factors, their re-stabilization is very different. The drops in *E*(*t*) in regime A can recover through *re-population*, over time scales approximately 1/*r*. By contrast, regime B exhibits system-size losses from which the only way to recover is *re-mutation*, over longer time scales determined by approximately 1/*m*^{2}, as the smallest cooperative group requires two mutations.

Comparing A to B reveals that higher selection strength in this case leads to lower fitness; had one mixed the stronger selected population B with the weaker selected one A, the former would be driven to extinction. The behaviour of B is similar to that expected in a classical prisoner's dilemma, which emerges from our model as a special case—survival of the fittest produces the globally least-fit outcome.

### 4.3. Formation of specialized bunches

As *r* is increased above the critical point *r*_{c} ∼ *m*/2 higher structures start to form. While the connectivity matrix *C* is sparse and random for *r* just below *r*_{c}, we start seeing metastable bunches at *r* > *r*_{c}, the number and stability of which increases with *r*. The sudden jump in the edge number we see in figure 2*d*,*e* is analogous to that found in [30].

For a very large window of selective strength we see that the system can only maintain certain discrete values of *E*. These are the stable configurations corresponding to an integer number of equal-sized bunches (*k*) given by the relation *E*_{k} = *N*^{2}(1 − 1/*k*), *k* = 0, 1, 2 , … , *k*_{max} (figure 2*d*). The maximum number of bunches *k*_{max} is determined by the mutation rate, *k*_{max} ∼ *N*/*m* (i.e. so that in steady state there is one mutation per bunch per time step). However, we have observed *k* transiently increasing to 50% higher than this value. Note that the degree of interdependence (and hence the average fitness) in the strong selection limit well exceeds that in the weak selection limit.

It is interesting that in the limit *r*, *m* ∼ 1, we see structures more complex than bunches. These include hierarchies (smaller bunches within a bunch), cycles (3 or more groups providing to one other), and hierarchies of cycles (cycles within a cycle). In this limit, the dynamics of *E*(*t*) still exhibits discrete steps similar to figure 2*e*, however with more possible metastable plateaus corresponding to unequal-sized matrix blocks.

### 4.4. Competition between bunches

With increasing *r* the fluctuation in the number of edges around the stable *k* starts increasing, and we see destructive competition similar to that near the phase boundary of B; however, now the competition is between the bunches rather than the individuals, which creates significant size differences between them. These fluctuations can lead to one bunch replacing another, causing *E _{k}* to make large transitions between different values of

*k*. Despite the apparent noise (figure 2

*f*) the dependence structure remains in a highly ordered state with high reciprocity (figure 3

*a*,

*b*, fifth panel). As

*r*is increased further we see that competition between bunches cause fluctuations comparable to the size of bunches, i.e. a small bunch can increase in size by spreading over others until another metastable structure re-evolves.

## 5. Conclusion

We have constructed a simple model that allows us to study the coevolution of self-replicating interdependent structures, and reported multiple evolutionary transitions as *β* and *r* are varied. This model is quite general and has very few assumptions—the fitness function is only assumed to be an arbitrary increasing function of *x* and there are only two relevant parameters governing the dynamics (selection strength *r*/*m* and relative benefit *b*/*c*) as the population size *N* does not make a qualitative difference as long as both *Np* and *r* are much smaller than *N*. Furthermore, the value *β* does not make a qualitative difference apart from whether it is larger or smaller than unity, and no *quantitative* difference if |1−*β*| is smaller than 1/*N*. Unlike the typical simplified models of evolutionary game theory, we do not assume that an individual's behaviour is the same towards all others (although some individuals can end up in a state where they give to all and receive from all). In this respect, the states allowed in this work are a generalization of the two-state models common in the literature. Thus, we hope that our model can serve as a guiding framework for understanding the emergence of sociality.

Even in this simple case, we observe a number of surprising and important phenomena. First, we report that even the weakest selection strengths (*m* ≫ *r*) can produce interdependence structures that are far from random. Thus, assumptions regarding ‘random interdependence’ invoked by neutral evolutionary arguments may be too strong [27–29]. Our second observation is the natural emergence of specialized bunches and multi-scale structures from the simple laws of evolutionary dynamics. As we probe the response of the system to various selection strengths, we see regimes of random interdependence, competition between nodes, cooperation between nodes (bunches) and competition between bunches.

Thirdly, we report that the regime *β* > 1, *r* > 0 does not ensure complex interdependence. There exists a ‘dead zone’ within the constructive regime (figure 2*c*) due to the competition between agents. This non-monotonic dependence of *E* on selection strength can have important implications in medicine. For example, biofilm populations may be induced into a less virulent non-cooperative state by decreasing the selective pressure, so that a cooperative film behaving as figure 2*e* evolves into an intermediate non-cooperative state behaving as figure 2*c*. This may be experimentally verified in *P. aeruginosa* by increasing the available iron while keeping their population constant by limiting their carbon source.

Another remarkable result is the non-monotonic connection between anti-social drug dose and the successful annihilation of cooperativeness. Indeed, the model exhibits a ‘contagion’ effect which allows the manipulation of a few individuals to have population-wide effects. It has been noted that introducing several selfish mutants (or using an anti-social drug effective on a few individuals) may be far more effective than manipulating an entire population [23] and is consistent with experimental observations [5,24,25].

Finally, when mutation rate is set to zero we observe that the behaviour of *E*(*t*) resembles that of classical population dynamics. The dynamics between providers and receivers becomes qualitatively similar to that between the predators and prey of a Lotka–Volterra-type system. Starting from a randomized connectivity matrix and setting *p* = 0, it is common for *E*(*t*) to reach a fixed value and oscillate around it.

Owing to its generality and applicability, our model has room for many natural extensions. For example, the distribution of parameters and fitness functions in a more realistic model could include spatial, temporal and individual heterogeneity. Further, the quantities *p*, *b* and *c* can be dynamic as they are themselves, to an extent, subject to evolutionary forces. This can lead to a very interesting set of potential future studies exploring connections between interdependence and evolvability/efficiency. Another factor not taken into account here is the possibility of the change in population size due to statistical fluctuations (e.g. due to a time-dependent energy input, or infection/predation). Such extensions would be appropriate to address systems in ecology, structured biological population, and provide insight into complicated social trends.

## Competing Interests

We declare we have no competing interests.

## Funding

This work was funded by the Wyss Institute for Biologically Inspired Engineering, the Harvard Kavli Institute for Bio-nano Science and Technology, the MacArthur Foundation and government support (A.I.) under FA9550-11-C-0028 awarded by the Department of Defense, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate Fellowship, 32CFR168a.

- Received January 17, 2015.
- Accepted May 7, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.