Shape matters: lifecycle of cooperative patches promotes cooperation in bulky populations.

Natural cooperative systems take many forms, ranging from one-dimensional cyanobacteria arrays to fractal-like biofilms. We use in silico experimental systems to study a previously overlooked factor in the evolution of cooperation, physical shape of the population. We compare the emergence and maintenance of cooperation in populations of digital organisms that inhabit bulky (100 × 100 cells) or slender (4 × 2500) toroidal grids. Although more isolated subpopulations of secretors in a slender population could be expected to favor cooperation, we find the opposite: secretion evolves to higher levels in bulky populations. We identify the mechanistic explanation for the shape effect by analyzing the lifecycle and dynamics of cooperator patches, from their emergence and growth, to invasion by noncooperators and extinction. Because they are constrained by the population shape, the cooperator patches expand less in slender than in bulky populations, leading to fewer cooperators, less public good secretion, and generally lower cooperation. The patch dynamics and mechanisms of shape effect are robust across several digital cooperation systems and independent of the underlying basis for cooperation (public good secretion or a cooperation game). Our results urge for a greater consideration of population shape in the study of the evolution of cooperation across experimental and modeling systems.

Evolution is often perceived as being “red in tooth and claw” and its processes interpreted through some variation of the “survival of the fittest” principle (Russe 1999). And yet, nature teems with instances of seemingly altruistic acts: organisms help one another, paying a direct energetic or reproductive cost without receiving an obvious or immediate benefit (Axelrod 1984; Dugatkin 1997; Rankin et al. 2007). Given that each vampire bat that regurgitates blood to a nest mate effectively challenges one of the basic tenants of modern evolutionary theory, it is no wonder that the evolution of cooperation remains a popular and contentious topic (Wingreen and Levin 2006). Additionally, instances and theories of cooperation in nature can easily inspire reflections concerning human society, maintaining both scientific and public interest (West et al. 2011).

Currently, a large fraction of empirical research on cooperation is conducted using microbial systems (Crespi 2001; Velicer 2003; West et al. 2006; Chuang et al. 2009; Dimitriu et al. 2014). Microbes exhibit a stunning variety of complex social interactions, such as punishment and spite (Kiers et al. 2003; Inglis et al. 2009), but remain relatively easy to perpetuate and manipulate. Although the evolution of cooperation never fully disappeared from the research focus, it is currently going through a renaissance because of bacteria such as Pseudomonas (Rainey and Rainey 2003; Griffin et al. 2004; Kümmerli et al. 2009; Popat et al. 2012). Although prisoner's dilemma tournaments motivated much of the research in 1980s and early 90s (Axelrod 1984; Nowak and May 1992), modern cooperation work is largely driven by the observations and experiments on microbial communities. At the core of our understanding of the evolution of cooperation is its first formal explanation, the succinct but powerful Hamilton's rule: rb > c (Hamilton 1963, 1964). Simply, cooperators will prosper if the relatedness between the interacting individuals (r), multiplied by the benefit experienced by the cooperators (b), is greater than the cost they pay for performing the cooperative action (c). Hamilton's rule continues to motivate research into quantifying the costs and benefits of cooperation (Brown and Taddei 2007), ways to measure relatedness (Pepper 2000; Damore and Gore 2012), and specific conditions that may hinder or aid its applicability (Lehmann and Keller 2006; Chuang et al. 2010). In our study we rely on mechanistic explanations, but do measure relatedness and show it is consistent with the phenomena we observed.

In parallel with, and as a result of theoretical and experimental advancements, many properties and processes favoring or disfavoring the evolution of cooperation have been identified. For example, we know that population viscosity generally favors cooperation (Kümmerli et al. 2008, 2009), that diffusion and degradation of the public good nonlinearly affect public good production (Misevic et al. 2012), and that costs and benefits of cooperation are affected by public good durability (Kümmerli and Brown 2010). Sometimes helping, at other times hindering the growth of the field, the passionate debate on the relative importance of group versus kin selection (Nowak et al. 2010; Abbot et al. 2011) frames all the work on the evolution of cooperation, but remains unresolved. Here we examine a typically overlooked parameter, the physical shape of a population, which reveals unexpected evolutionary dynamics specifically when coupled with cooperation.

Much has been made out of the importance of population structure for the evolution of cooperation, both in theory and experiments (Nowak and May 1992; Kümmerli et al. 2008, 2009; Perc et al. 2013). The discourse typically focuses on increased assortment (relatedness), the elevated probability that cooperators will encounter each other rather than noncooperators (Dimitriu et al. 2014). Significantly less attention is given to the exact location of the interacting individuals in the physical space and in the population they constitute. By default, starting with the early evolution of cooperation literature, everything was always happening on a square lattice grid (Nowak and May 1992). There are some exceptions, including, for example, the work on honeycomb (hexagonal) grids, which were shown to promote cooperative prisoner's dilemma strategies more often than square ones (Hauert and Szabó 2003). However, such setups were treated primarily as interesting but nonessential extensions of the basic, square grid models and simulations. The game theory studies of cooperation have also considered one-dimensional and three-dimensional populations (Nowak et al. 1994a), but few general results exist. For example, in a very specific case of prisoner's dilemma with only two strategies (all defect and tit for tat), Nakamaru and colleagues claim that “two-dimensional lattice model behaves somewhat in between the one-dimensional model and the complete mixing model” but also show that dimensionality does not matter when using the mean field approximation (Nakamaru et al. 1997). More recently, increased computational power allowed some simulations to consider continuous, rather than discrete space, further adding to the biological realism of the in silico cooperation system (Nadell et al. 2010). Perhaps the most interesting extension is based on graph theory, where each individual is represented by a vertex and the potential interactions between individuals by edges connecting those vertices. Although not a new concept (Myerson 1977), it has the potential to connect and extend spatial models of different dimensionality, by adding more connections between individuals (Vukov et al. 2008). As real networks of interactions are starting to be fully mapped in nature (Salathé et al. 2010), graph theory is quickly becoming a promising methodology for analysis of cooperation dynamics (Ohtsuki et al. 2006; Szabo and Fath 2007). In this work, we consider differently shaped two-dimensional lattices in three distinct computational models of varying complexity. In spite of the extensive literature on cooperation in space, to our knowledge, there have been no studies directly addressing the effect of the population shape on the dynamics and outcome of the evolution of cooperation. Given a variety of sizes and shapes occupied by natural cooperative system, this study has potential to inform and affect a broad range of the theoretical and experimental studies of cooperation.

Material and Methods

Our approach to studying the evolution of cooperation rests on the foundations of both social microbiology and biological modeling. We use three different simulations systems: Aevol (public good-based cooperation, complex genetic), Aevol-lite (public good-based cooperation, single binary locus genome), and CAevol (simple probabilistic cellular-automata model of cooperation based on prisoner's dilemma). We fully describe the relevant properties all three models here, and provide additional details in the Supporting Information.

AEVOL EXPERIMENTAL SYSTEM

Aevol is an open-source simulation platform that uses the genetic algorithm approach to individual-based modeling of evolutionary processes. Digital organisms that cooperate in the Aevol in silico system have double-stranded genomes and nontrivial genotype–phenotype maps (Misevic et al. 2012; Frénoy et al. 2013). Although certainly less complex than even the simplest bacteria, they go well beyond cooperate/defect strategies in capturing the natural evolutionary processes. Cooperation in Aevol is based on a costly public good that diffuses and degrades in a spatially structured environment. The interactions among organisms occur over both ecological and evolutionary time scales, and our experiments span tens of thousands of generations of evolution via natural selection. Inspired by experimental microbial evolution research, Aevol is capable of investigating scenarios and regions of parameter space that may be important in nature but are unattainable in the laboratory (Hindré et al. 2012). Additional information and the program files can be found at http://www.aevol.fr as well as in the previous publications that have used this platform (Knibbe et al. 2006, 2007; Misevic et al. 2012; Frénoy et al. 2013). Aevol is an abstract instantiation of evolution and specific features of its in silico biology are not meant to directly model any particular organism but are inspired and motivated by known properties of microbial genetics. We treat it as an experimental system in its own right and analyze all the data accordingly (Peck 2004). Aevol is well suited for this study because it captures not only the dynamics of cooperation phenotypes, but also the genetic-level interaction between cooperation and metabolic traits. These interactions have recently been shown to strongly influence the evolution and maintenance of cooperation (Morgan et al. 2012; Frénoy et al. 2013)

Genotype–phenotype fitness mapping

Each digital organism has a genotype consisting of a double-stranded binary number string, typically 5–10,000 bases long. What sets Aevol apart from other frequently used in silico systems, such as Avida (Adami 2006), is the bioinspired genotype–phenotype fitness map. The phenotype of an organism is a result of transcription and translation of its genes, based on locating within the genome, the predetermined promoter sequence, followed by the consensus Shine–Dalgarno sequence, the gene itself, and the palindromic terminator sequences. More broadly speaking, the phenotype of a digital organism is a collection of traits and their values. Traits are arranged on a trait axis and associated with real numbers between 0 and 1. We thus represent a protein as a triangle resting on the trait axis, and the organism's phenotype as the sum of all proteins encoded in its digital genome (Fig.1). Each protein targets a particular trait, but also affects the neighboring traits, to a lesser degree. Genetic code, including start, stop, and six other triplet codons, specifies the (mean, height, width) triplet that describes the size and location of a protein triangle (see Fig. S1). Mean determines the trait that is most affected by a protein, height the value of that trait, and width the most distant trait on the axis that this protein contributes to. Traits that are close in the phenotypic space are not necessarily close in the genotypic space and there are many possible ways to encode same or similar proteins. However, each gene has only one possible translation into a protein.

Figure 1

Schematic representation of phenotype and fitness in Aevol. (A) All proteins of an organism, each represented by a triangle. A protein is numerically defined by three values, m (location of the triangle on the phenotypic axis, the primary trait this protein affects), h (triangle height, level at which it affects the trait), and w (width of the triangle, representing the pleiotropy with neighboring traits). Depending on its location on the phenotypic trait axis, a triangle either contributes to organism's fitness directly (blue) or modifies the amount of the public good secreted (orange). (B) Gaussians used to describe the best possible phenotype, the target phenotype, colored by the location on the phenotypic axis. (C) Organism's phenotype, the jagged black line, which is the sum of triangles from (A), and the target phenotype, which is the smooth curve calculated as a sum of Gaussians from (B). The difference between metabolic proteins and the target protein levels (area shaded blue) is a direct component of individual fitness, whereas the same difference for secretion proteins (area shaded orange) determines the amount of public good individual secretes (see methods and Supporting Information for the numerical formula and more details on genotype to phenotype to fitness mapping.)

The full phenotype is composed of two types of traits, relating either to metabolism or secretion. Metabolic traits affect fitness directly, thus changing the probability of the organism's offspring populating the next generation. Secretion traits determine the amount of the public good molecule that is secreted by an organism. Intuitively, the direct fitness contribution of the metabolism is inversely proportional to the distance between the organism's phenotype and the target phenotype. The target phenotype is an arbitrarily chosen collection of traits and their values, that is, the optimal phenotype for the environment (Fig.1B, C). It is specified by adding the functions of the form y = H exp(−(x − M)2/(2W2)), where (H,M,W) = {(0.2, 0.3, 0.04), (0.3, 0.2, 0.02), (0.2, 0.7, 0.02), (0.3, 0.8, 0.04)} and was constant and identical across all the experiments (see Supporting Information for discussion of the parameter choice). The amount of secreted public good is inversely proportional to the distance between the target secretion phenotype and the phenotype curve formed by all secretion traits encoded in the organism's genome. The metabolic and the public good portion of the fitness are combined multiplicatively. Precisely, the fitness of a digital organism is calculated as W × (1 + PG – C × S), where W is organism's metabolic fitness, PG is the amount of the public good molecule present in the organism's grid cell, C is the cost of secreting a unit of public good, and S is the amount of public good molecule that the organism secretes. Additionally, W and S = e, where Gm is the area of the gap between the target phenotype and organism's phenotype for metabolism (Fig.1C, shaded blue), Gs is the gap between the target phenotype and organism's phenotype for secretion (Fig.1C, shaded orange), As the area under the secretion part of the target phenotype, and a the selection pressure constant set to 0.7 for all the experiments (see Supporting Information for additional discussion of fitness calculation and ranges for parameters and variables). The fitness, amount of secretion as well as diffusion, degradation, and other parameters are reported in dimensionless, arbitrary units. In principle, using various constants we can adjust the units and report fitness as energy, in joules, or secretion as the number of public good molecules, but this could be interpreted as a superficial and ultimately unnecessary modification.

Cooperation in Aevol

The cooperative trait in Aevol, secretion of public good molecules, is continuous, unlike in classical game-theoretical models of cooperation (Axelrod 1984). We thus consider the average amount of secretion, rather than a proportion of cooperators in the population, as the measure of the level of cooperation. Public good diffuses at the rate of 0.05, meaning that at each generation 5% of the public good molecules present in a cell diffuses into each of the neighboring cells. At the same time, public good degrades at a rate of 0.1, meaning that in every generation, 10% of the public good present in each population location degrades. After both processes, 54% of the original amount of public good remains at that location for the next generation. The cost of secretion was 0.03, meaning that for secreting S units of public good that may increase the fitness of its offspring and their neighbors in the next generation, an organisms pays immediate fitness cost equal to 0.03S, as detailed in the fitness formula in the section above. Previous work with Aevol demonstrates that this setup does create the classical evolutionary dilemma of cooperation (Misevic et al. 2012; Frénoy et al. 2013).

Population structure

All populations contained 10,000 individuals positioned in a quadrilateral grid. As there is no explicit death of individuals during a generation, all positions in a gird are occupied, only to be effectively all emptied during the reproduction step in the lifecycle. Each digital organism has exactly eight neighbors (the classical 3 × 3 Moore neighborhood, Fig. S2). Other plane tessellations, based on anything from regular hexagons to Escher's lizards, are of course possible, but the grid of squares persists as a default setup due to its intuitive simplicity. To avoid edge effects, our digital world is folded onto itself, becoming a torus. Depending on the dimensions of the starting grid, the resulting toroidal world may be akin to a bulky doughnut or a slender bicycle tire.

Aevol lifecycle

Generations are nonoverlapping in Aevol and each consists of seven steps performed in order: fitness evaluation, secretion, diffusion, degradation, selection, and reproduction with mutations. A generation starts with evaluating the fitness of all organisms, as described above. Only afterwards do individuals deposit in their grid-cell location an amount of public good determined by their genotype and, in turn, their phenotype. The public good then degrades and diffuses as already described. Due to this order of steps, an organism cannot directly benefit from its own secretion. Selection starts by ranking the organisms in each 3 × 3 Moore neighborhood (Fig. S2) based on their fitness, where the fittest organism receives the rank 9 and the least fit the rank 1. An individual with the rank R has (a – 1) × a9 − / (a9 – 1) probability of reproducing into the central position of the neighborhood, where a is the selection pressure constant and in these experiments is set to 0.7. Put differently, the individual populating the location i in generation n + 1 is an offspring of one of the nine individuals in the neighborhood of i in generation n, the parent having been selected based on the probability of reproduction. The benefit of a rank-based regime is that even small fitness differences result in significant selection pressures, extending the period during which a population is continuously adapting. The strength of selection decreases with higher values of a. During reproduction, the binary sequence of a new individual changes randomly, with a fixed probability. Specifically, the mutation rates were 10−5 per base pair for base pair substitutions and small insertions/deletions (less than 6 base pairs), and 10−6 for inversions, translocations, and large insertions/deletions (see Supporting Information for discussion of the parameter choice). Reproduction with mutations completes the cycle of a single generation in Aevol.

There are no cooperators at the beginning of the simulation. For each replicate simulation, we repeatedly generate 5000 long random strings of 0's and 1's, until one of those strings contains a functional metabolic gene (see Genetic Encoding in Aevol section in Supporting Information). The initial population consists of 10,000 clonal replicates of this randomly generated organism. Finally, in our main experiments, the lifecycle is repeated 50,000 times, which directly corresponds to 50,000 generations of evolution.

AEVOL-LITE EXPERIMENTAL SYSTEM

Aevol-lite is a simplification of the Aevol system, implemented in Matlab. It preserves the full setup and parameterization of Aevol, but eliminates the underlying complex genetics. Instead of a double-stranded binary string, which is transcribed and translated, a genome of an Aevol-lite individual contains a single number, 0 or 1, which determines whether they are a nonsecretor or secretor. There are no metabolism traits, and unlike with Aevol, the cooperation trait is not continuous, but binary. The lifecycle of fitness evaluation, secretion, diffusion, degradation, selection, and reproduction with mutations is the same as before. The population size is 10,000 individuals, all of which are nonsecretors at the beginning of the simulation. Mutation rate between the two types of individuals is set to 0.005 per generation, meaning there is a 0.5% chance that a cooperator will mutate into noncooperator, and vice versa. Diffusion rate was 0.05 and degradation rate was 0.1, the population shapes were 100 × 100 and 4 × 2500, as before. Cooperators secreted 0.4 units of public good and their fitness was calculated as 1 + PG − C, where PG is the amount of the public good molecule present in the organism's grid cell, and C is the cost of secreting, here set to 0.1. For the sake of simplicity, the fittest organism from each 3 × 3 neighborhood was selected to reproduce. The broad parameter range we have tested displayed qualitatively the same evolutionary dynamics and the parameter values used here are a representative example.

CAEVOL EXPERIMENTAL SYSTEM

CAevol model is a further simplification of the Aevol system, based on Cellular Automata dynamics (thus the name) and is implemented in Matlab. The cooperation is encoded by 0 or 1, marking an individual as a noncooperator (defector) or cooperator, respectively. However, there is no public good, and cooperation is instead based on the classic prisoner's dilemma game. During a generation, each individual interacts with all of its eight neighbors (found in the 3 × 3 Moore neighborhood, as before), and collects payoffs according to the following simplified payoff matrix: (CC, CD, DC, DD) = (1, 0, T, 0), where the letter pair signifies the organisms interacting (C for cooperator and D for defector) and the number determines the payoff to the first individual in the pair. T, typically called temptation, is the payoff for the defector when interacting with the cheater, and was set to 1.03 our experiments. The fitness of an individual is simply the sum of its payoffs. In prisoner's dilemma scenarios, selecting the fittest individual to leave an offspring in the next generation is known to produce spatial artifacts (Nowak et al. 1994b; Oliphant 1994) and thus we used stochastic, fitness-proportionate selection, where the probability to reproduce is proportional to individual's fitness. The fitness evaluation selection reproduction lifecycle is unchanged from Aevol-lite. There are no cooperators at the start of the experiments, and the mutation rate is 0.005 per generation, as before.

EXPERIMENTAL DESIGN

Until now we have focused on the general properties of different systems. Here we present on the specific protocols for the different experiments.

Population shape experiments. In all experiments, using Aevol, Aevol-lite, or CAevol, the population shapes are derived from and denoted by the dimensions of their unfolded grids, namely 100 × 100 (bulky) and 4 × 2500 (slender). Fifty independent populations evolved in Aevol for 50,000 generations at each shape treatment (see Supporting Information for the discussion of experiment duration). The number of individuals was always constant and equal to 10,000, independently of the world shape (there were no empty positions in the grid). The shape of an evolving population is thus a direct consequence of the toroidal grid that it inhabits. This setup was reused in Aevol-lite and CAevol, but due to quicker evolution and less variation between replicates, it was sufficient to evolve 30 populations for 5000 generations in each Aevol-lite and CAevol shape treatment. We use the average amount of public good secreted by an individual in an Aevol population and the proportion of cooperators in Aevol-lite and CAevol to quantify evolutionary dynamics and outcomes.

No secretion experiments. To investigate the connection between shape effect and cooperation, we evolved 50 Aevol populations in which individuals could not secrete or cooperate for 50,000 generations. The phenotype target values typically associated with secretion (right side of the Fig.1 panels) remained numerically unchanged, but rather than determining the amount of public good secreted, they were considered as metabolic traits and contributed directly to fitness.

Reshaping experiments. In these experiments, we test the robustness of different cooperation levels to a change in population shape. The Aevol populations were reshaped by first unfolding a torus into a flat grid, and then, in the case of transforming a bulky into a slender population, cutting away 4 × 100 strips of the grid and putting them one next to the other thus forming a 4 × 2500 grid to be refolded into a torus. Similarly, for slender populations, we cut 4 × 100 strips that we then stacked on top of each other to form a 100 × 100 grid (see Supporting Information for mathematical formulation of the transformation, and Fig.2B for an example 4 × 2500 population similarly reshaped for ease of display). Both the organisms themselves and the information about the amount of the secreted molecule present in their respective grid cells were conserved in the reshaped populations. We compared the final amount of secretion among 50 populations in each of the four treatments characterized by the amount of time spent in each population shape: (1) slender for 70,000, (2) bulky for 70,000, (3) slender for 50,000 then bulky for 20,000, and (4) bulky for 50,000 then slender for 20,000 generations.

Figure 2

Within-population diversity in (A) a bulky and (B) slender Aevol population. Each square represents an individual and individuals that secrete the same amount share the same color. Black lines mark the edges of the population that connect to form a torus. The slender population is separated into 4 × 100 strips that are stacked onto each other to display it at the same scale as the bulky population. Based on the amount secreted, there are 579 different organism types in the bulky and 1210 in the slender population. The colors are assigned at random to better differentiate types with numerically similar secretion.

LIFECYCLE AND DYNAMIC OF COOPERATION PATCHES

To quantify the mechanism underlying the difference in cooperation between bulky and slender populations, we analyzed the dynamics of a cooperator patch in Aevol-lite. Specifically, we first ran controlled single patch experiments by placing a single cooperator in the center of the population of noncooperators and deposited six units of secretion in the same location, enough for the expanding patch to be formed. We propagate the population for 100 generations, branching out separate experiments during it. At each of those generations, we consider all possible mutations from single cooperator to a noncooperator. After we introduce one of the mutations, we propagate the population for an additional 400 generations, without any further mutations. We used the symmetrical nature of the expanding patches to examine only unique scenarios, significantly decreasing the computational load, but still analyzing thousands of expanding populations. For each of them, we counted the number of cooperators present thought the length of the experiment and weighted it by the probability of the scenario occurring. The scenario probability is simply the product of the probability of the mutation occurring in that generation and the probability there have been no previous mutations, both a function of the mutation rate. We examined 99 mutation rates evenly spaced between 0.001 and 0.1 in addition to those two.

Because these experiments examine only the lifecycle and dynamics of single, isolated patches, we wanted to compare them with regular, unconstrained Aevol-lite populations. We thus evolved 30 populations for 5000 generations in Aevol-lite, using the standard setup from the initial population shape experiments, at seven similar mutations rates, namely 0.001, 0003, 0.005, 0.01, 0.03, 0.05, and 0.1. In this case, we calculate and report the ratio between the percentage of cheaters and cooperators as a measure of the difference between slender and bulky populations. The standard deviation of the ratio for each mutation rate was calculated as the square root of the variance, approximated as , where x and y are variables for the percentage of cooperators in bulky and slender populations, respectively.

RELATEDNESS METRICS

In a digital system such as Aevol, it would be in principle possible to record the full pedigree of an organism and extract the exact degree of relatedness or similarity between any two organisms. However, given the population size, length of evolution, and number of replicates, neither the storage nor analysis of such exhaustive data is technically possible. Instead, based on previous work (Taylor and Frank 1996; Pepper 2000; Damore and Gore 2012), we compute the coefficient of correlation as a simple and intuitive relatedness metrics using data on the relevant trait, the public good secreted by each organism. For each population, we compute the correlation between the amount of public good secreted by an individual and the average amount secreted by individuals in its 3 × 3 Moore neighborhood. The higher values of the correlation coefficient indicate higher level of relatedness among organisms in a population.

POPULATION CLUSTERING

Although just two organism types exist in Aevol-lite and CAevol, the number of cooperator types is much larger in Aevol populations, due to the continuous nature of the public good secretion. To examine the population composition and patchiness, we identify the clusters of similar individual using the Markov Cluster Algorithm, MCL (Van Dongen 2008). MCL considers an Aevol population as a graph in which individuals represent the vertices and edges are drawn between each two neighboring organisms. The weight of an edge is , where s1 and s2 are the amount of public good secretion for two organisms connected by the edge and max(s) is the maximal amount of secretion by any organism in the population. MLC has a single tunable parameter, inflation (R), which determines the rescaling of the stochastic matrix of the Markov process used to identify the clusters in the graph. It primarily affects the granularity of the clusters, with higher values of R producing smaller clusters than lower values. In preliminary tests, we varied the inflation parameter, before settling on 1.5 based on comparison of clustering of populations that evolved for 50,000 generations and their randomized versions. We apply a random permutation to the organism location indices to derive a randomized population from the regular one. For each population we recorded the number and size (number of vertices) of clusters calculated by MCL. For R = 1.5, the number of clusters was smaller and they were larger in regular than in randomized populations, which was an internal control for our methodology: when the number of clusters is smaller in regular than in shuffled populations, the reported clusters represent a meaningful grouping of individuals.

STATISTICS

The statistical analysis was done using the Matlab software package. For each comparison between different treatments, we first test the null hypothesis that the data are coming from a normally distributed population, using Lilliefors test. The details of the test and P values are reported in the Supporting Information. We used Wilcoxon rank sum test when the hypothesis of normality was rejected and two-sample Student's t-test for data when it could not be rejected. The exception was the comparison of clustering results between regular and shuffled populations where Wilcoxon sign-rank test for related samples was used, due to the shared ancestry of the populations being compared.

Results

POPULATION SHAPE EXPERIMENTS

We start by quantifying the central result of this study: there is more cooperation in bulky than in slender populations. The first data we look at are from Aevol experiments with populations of two different shapes, namely 100 × 100 (bulky) and 4 × 2500 (slender). We find a significant difference in the evolved secretion amount between the two shapes (Wilcoxon rank sum, P < 10−12, Fig.3A). Similar results are obtained for Aevol-lite (Wilcoxon rank sum, P < 10−12, Fig.3A) and CAevol experiments (Fig.3C, Wilcoxon rank sum, P < 10−17). All the comparisons were done using data from the final time point of the experiments. However, the effect of shape was present throughout the experiments and was robust to many parameter combinations (see Figs. S2, S3, S6–S8). Over the entire course of the experiment, slender populations in Aevol secrete on average 10.1% less than bulky ones, which contributes to them also being 12.5% less fit, indicating a biological and not just statistical significance (see Supporting Information for further discussion).

Figure 3

The difference in cooperation between populations that evolved in bulky (100 × 100) versus slender (4 × 2500) populations in (A) Aevol, (B) Aevol-lite, and (C) CAevol. Cooperation was quantified at the final generation of the experiment as the average amount of public good secreted in Aevol and the percentage of cooperators present in the population in Aevol-lite and CAevol experiments, which are represented here by boxplots. The horizontal line is the median, the box edges are the 25th (q1) and the 75th (q3) percentile, the whiskers mark the most extreme data points still smaller than q3 + w(q3 − q1) and larger than q1 − w(q3 − q3), where w = 1.5. Differences between shapes in each panel are significant (see main text).

POPULATION DYNAMICS OF PATCHES

To identify the explanation for the effect of population shape on secretion, we examined the dynamics of cooperating populations. It was immediately evident that the population averages from Figure3 do not capture the existing within-population diversity (Fig.2). In Aevol there are hundreds of different cooperation phenotypes, and even in Aevol-lite and CAevol, populations are a constantly changing mixture of cooperators and cheaters (see Video S1). The individual types are in continued competition, which is reflected by the perpetually expanding and shrinking patches they form. Although patch dynamics is difficult to follow in Aevol due to large number of organism types, it is visually striking in the two other systems. In Aevol-lite we can clearly identify the following stages of a cooperative patch lifecycle (Fig.4): (1) emergence, persistent group of cooperators appears; (2) expansion, cooperator offspring are invading neighboring locations; (3) invasion, at least one noncooperator mutant appears within the cooperator patch; (4) decline, an expanding noncooperator patch forms; and (5) extinction, noncooperator patch completely takes over, no individuals from the original cooperator patch remain.

Figure 4

Example of a full lifecycle of a cooperator patch in (A) a bulky and (B) slender Aevol-lite population. Each square represents an individual, with cooperators drawn as dark and noncooperators as white cells. Each square in panel (A) and rectangle in panel (B) represents a snapshot of a portion of the population, with the red number in the top left corner indicating the generation. In both panels a cooperator patch can be tracked from its appearance, through initial expansion, the appearance of noncooperator, cheater patches within it, and until it disappears by being completely taken over by the cheater patches.

There are several scenarios for the emergence of cooperator patches, including multiple cooperators appearing by mutation in a neighborhood in the same or consecutive generations or even a single cooperator reproducing by chance, in spite of having lower fitness than its noncooperating neighbors. These complex, stochastic scenarios make analytical, equation-based quantification of the patch dynamics difficult, which is why we simulated the isolated patches in Aevol-lite. Following initial appearance, a cooperator patch in a bulky population can expand on all four sides, potentially adding dozens of new cooperators in a single generation (e.g., Fig.4A, from generation 343 to 344). On the other hand, after reaching the size of about 16 individuals (4 × 4 square), the corresponding patch in a slender population is constrained by the population shape and cannot grow by more than eight cooperators per generations, four on each side (e.g., Fig.4B, from generation 619 to 620). Due to faster patch growth, we expect the noncooperators to appear within a patch in bulky sooner than in slim populations. The invasion of noncooperators cannot happen from the outside, but requires a mutation of an individual in the growing patch, because of the accumulation of the public good and the corresponding benefit it confers. This will happen at different times in populations of different shape. For example, at mutation rate of 0.0005, there is more than a 99.9% probability of invasion starting within the first 68 generations in slim and 32 generations in bulky populations. The nonsecretors likely start their takeover earlier in bulky than in slim populations. However, cooperative patches in bulky populations still grow faster and to a larger size, because they are not constrained by population shape. The difference in the dynamics and lifecycle of cooperative patches directly leads to a greater number of cooperators being present in bulky than in slender populations–this is a simple, effective, mechanistic explanation for the effect of population shape on cooperation.

We directly test the patch dynamics hypothesis by tracking expanding patches from their emergence to extinction and summing the total number of cooperators. We do so by isolating a single cooperative patch and only one noncooperative mutant in a reduced Aevol-lite setup. Independently of the mutation rate, the total number of cooperators during the patch lifetime is greater in bulky than in slender populations (Fig.5A), directly supporting our analysis. Further confirmation comes from comparing the ratio of cooperators in bulky and slender populations at different mutation rates in the full Aevol-lite simulations, with many patches simultaneously going through their lifecycles (Fig.5B). As seen previously, there are more cooperators in bulky than in slender populations. Moreover, the ratio of cooperators to cheaters was generally greater at smaller mutation rates. The change in the ratio can be interpreted as being a direct consequence of the greater difference in patch size between slender and bulky populations at lower mutation rates, and supports our explanation of the effect of population shape on cooperation.

Figure 5

Quantification of patch lifecycle and dynamics. (A) Total number of cooperators present during the lifetime of a patch in a bulky (blue, triangles) and slender (green, circles) population as a function of mutation rate. Lifecycle of a single, cooperative patch was examined in isolation, using a reduced version of Aevol-lite. Mutation rate corresponds to the probability of cooperator mutating into a cheater. All possible scenarios for appearance of a single cheater mutant are considered and the number of cooperators summed over time. (B) Ratio of the mean percentage of cooperators in bulky versus slender populations for different mutation rates after 5000 generations of evolution in Aevol-lite. Multiple patches of cheaters and cooperators exist and the mutation rate corresponds to the probability of an offspring having a different phenotype than the parent. The error bars represent one standard deviation.

PATCHINESS IN AEVOL

Based on the amount of secretion, we evaluate the patchiness of Aevol populations. To do so, we have clustered individuals in the bulky and slender populations at the end of our experiments using the MCL algorithm described in the methods. Comparing the number and size of detected clusters from the actual populations with the ones of the shuffled populations confirmed that there is some significant clustering of individuals. Populations were shuffled by randomizing the positions of all individuals. There are significantly fewer and significantly larger clusters in the actual than in the shuffled populations (P < 10−21 for all comparisons, paired two-sample t-test, Fig.6). Finally, we confirm that the bulky populations have both fewer (P < 10−58, two-sample t-test) and larger clusters (P < 10−41, two-sample t-test, Fig.6) consistent with the mechanisms of patch spread we have described in Aevol-lite.

Figure 6

Clustering of Aevol populations. We used the MCL algorithm to identify (A) the number and (B) the average size of clusters in evolved populations. Individuals were clustered based on the amount of the public good they secrete. We evaluated the significance of our clustering method by comparing the clustering of regular evolved populations (blue, first and third boxplot from the left, in each panel) and populations constructed from them by randomly permuting the location of all individuals (green, second and fourth boxplot from the left). Stars mark statistically highly significant differences between clustering results (see main text for statistics). The boxplots have the same properties as in Figure3, and the points outside of the whiskers are considered outliers and marked with plus symbols (+).

RELATEDNESS IN AEVOL

We measured relatedness by calculating the coefficient of correlation between the neighboring individuals’ secretion phenotype in the evolved populations. Higher coefficient of correlation observed in bulky populations (Wilcoxon rank sum P < 10−17, Fig.7A) directly points toward higher relatedness in bulky than in slender populations. The relatedness difference is consistent with the corresponding smaller, more numerous clusters in slender populations, arising from lower similarity among neighboring individuals. To investigate the potential causal connection between patch dynamics, relatedness, and cooperation, we repeated the analysis using data from Aevol experiments in which secretion was not possible, calculating relatedness based on the metabolism phenotype. The difference in relatedness between the shapes we observed in cooperating populations was completely absent when individuals had no means of cooperating with each other (Wilcoxon rank sum P = 0.754, Fig.7B).

Figure 7

Relatedness in Aevol populations. Relatedness is calculated as the coefficient of correlation between secretion of an individual and average secretion in its neighborhood. We measure relatedness in bulky and slender populations in which secretion, and thus also cooperation, are (A) possible or (B) not possible. The boxplots have the same properties as in Figure3.

RESHAPING EXPERIMENTS IN AEVOL

So far we have established that the population shape consistently affects cooperation. Here we investigate the potential genetic constraints the evolution in different shapes may convey. Namely, if the shape of the population changes, would the level of cooperation change as well, or is it genetically “hard-coded” and difficult to change. Such a test is meaningful only in the Aevol system, where we know that organisms evolve genetic architecture that may affect cooperation evolvability (Frénoy et al. 2013), but not in the context of the Aevol-lite and CAevol simple genetic encoding.

We compared populations that spent 70,000 generations in the same shape against ones that switched the shape after 50,000 generations (Fig.8). The reshaped populations evolved significantly different levels of secretion, compared to populations that remained in their respective ancestral shapes (P < 10−5, in both comparisons, Wilcoxon rank sum). After the reshaping, the population secretion evolved toward the levels observed in populations that spent their entire evolutionary history in that shape. Specifically, at generation 70,000, there was no significant difference between either the two groups of 50 slender populations (P = 0.171, two-sample t-test), or the groups of bulky populations (P = 0.053, Wilcoxon rank sum). The difference between slender and bulky populations that were not reshaped remained highly significant (P < 10−11, Wilcoxon rank sum).

Figure 8

The average amount of public good secreted after 70,000 generations of evolution in either constant or changing population shapes. Fifty replicate populations with 100 × 100 (bulky) and 4 × 2500 (slender) shapes from the generation 50,000 of the initial Aevol experiments were propagated for additional 20,000 generations either in the same (blue, first and third boxplot from the left) or different (green, second and fourth boxplot from the left) shape. The different shape was slender for bulky populations and vice versa. Stars mark statistically highly significant differences between clustering results and “n.s.,” the nonsignificant ones (see main text for statistics). The boxplots have the same properties as in Figure3.

Discussion

Our experiments have a clear and easily stated outcome: organisms in bulky populations evolve to cooperate more than organisms living in slender populations. Simply put, shape matters. This outcome went against our intuition. At the core of the current understanding of cooperation is the idea that it will be selected for when cooperators preferentially interact with each other—for example, when cooperators form some type of refuges, regions that are free of noncooperators. Cooperator-rich regions, abundant in the secreted public good, effectively create a different niche, one that favors the maintenance of cooperation. Such islands of cooperation, already observed in some of the first game-theoretical simulations (Nowak and May 1992), would have fewer cells in contact with noncooperators in slender than in bulky populations just due to the constraints of grid geometry: in slender populations the borders between two types are more likely to exist only on two and not on all four sides of the neighborhood. However, because this prediction was incorrect, our experiments serve a double purpose: to provide an explanation of the observed effect and to revise our intuition about the evolution of cooperation via public good secretion.

To find a satisfying, mechanism-based answer we switched from Aevol to a simpler system, Aevol-lite. Because there are only two types of Aevol-lite individuals, secretors and nonsecretors, patterns are much easier to represent graphically, identify, and qualify than in Aevol (e.g., Fig.4, Video S1). We started from the observation that each population is a mix of secretors and nonsecretors, who are grouped in ever-expanding and shrinking patches. The dynamic of these patches holds the key to the population shape effect: in slender populations, the cooperator patches do not have as much space to expand as in the bulky populations, resulting in smaller number of cooperators over time. We quantified the difference and found that during the total patch lifetime there is greater number of cooperators present in bulky than in slender populations (Fig.5A). Considering one patch over time is generally equivalent to considering the entire population at one time point, as the population is composed of many patches in different stages of their lifecycle. Thus, we can conclude that the population shape effect is caused by differences in the area available to expanding patches of cooperators between slender and bulky populations.

Until now we have examined cooperation based on public good secretion. To test whether secretion is a necessary condition for the population shape effect we performed simulations with CAevol, where cooperation has significantly different properties. For example, in CAevol there is no memory of past population states and the fitness depends only on properties of the current population. Additionally, the speed of patch expansion is always the same, unlike in Aevol and Aevol-lite, where the expansion of secretor patches is delayed because it depends on accumulation and in turn the fitness benefit of the secreted public good (see, e.g., generations 336–338 in Fig.4A). However, the patch lifecycle and dynamics in CAevol actually causes greater differences in cooperation between slender and bulky populations, as observed in Figure3, because noncooperative patches cannot ever “catch-up” and overtake cooperative ones. Overall, the simulations with CAevol show that the shape effect is robust across evolving systems and not contingent on a particular cooperation mechanism.

Having identified the cause of the different level of cooperation in Aevol-lite and CAevol, we finally return to Aevol itself. Reconstructing the full lineage of all individuals and thus delineating the patches and their size remains computationally intractable. Given the continuous levels of secretion, and thus cooperation, distinct classes of cooperators and noncooperators cannot be identified to clearly characterize their interactions, as before. The same Aevol individual may be considered a “cheater” when compared to individuals secreting more, or a cooperator when interacting with individuals secreting less. However, the same mechanism governing the two-type cooperation dynamics in Aevol-lite should hold: constrained patch size of the expanding (more) cooperative strains leads to less overall cooperation. We identify patches of organisms with the MCL clustering algorithm and confirm that Aevol bulky populations indeed have a smaller number of patches than slender populations, and that those patches are larger, as predicted by Aevol-lite experiments (Fig.6). Furthermore, although we do not record “digital pedigrees” and are unable to find Haldane's proverbial two brothers or eighth cousins for Aevol organisms to cooperate with, we have measured relatedness and found higher relatedness in bulky than in slender populations, as would be expected based on the clustering data and Hamilton's rule.

Given the observation from simpler systems, we must conclude that the Aevol genotype–phenotype mapping is not critical for the shape effect, and in turn question the necessity of using Aevol in this study. However, while CAevol and Aevol-lite eliminate the possibility that complexity single handedly drives the effect, they cannot assure us that it would not hinder it. Genetic hitchhiking, epistasis between secretion and nonsecretion genes, evolving switching rates between cooperators and noncooperators are just some of the high-level processes that may shape the evolution of cooperation in bulky and slim populations. They all are dependent on the genetic architecture and thus can play a role in Aevol but not Aevol-lite and CAevol. The reshaping experiments specifically show that shape does not irreversibly modify cooperation. Even though Aevol individuals can exhibit constrained, hard to modify genetic architecture (Frénoy et al. 2013), average cooperation level does change after the change in population shape. Organisms adapt to the present shape, as they would to any other environmental constraint, such as public good diffusion or degradation. The same way a given set of public good physical properties determines the “optimal,” stable amount of cooperation (Brown and Taddei 2007; Kümmerli and Brown 2010; Wintermute and Silver 2010; Misevic et al. 2012), so does the population shape. Overall, the presence of the shape effect in Aevol populations increases the likelihood of also observing it in even more complex, biological systems.

It would have been convenient to explain the shape effect via relatedness alone, the effect somehow being directly created by the population shape. However, the analysis of populations that do not cooperate and do not show a difference in relatedness directly contradicts arguments purely based on Hamilton's rule (Fig.7B). Rather than “shape causes relatedness that causes cooperation,” we must consider a closer, more integrated relationship between these population properties: relatedness is not the result of the population shape alone, but of the dynamics of cooperation within a shape. Without secretion, the patches of different organisms are not “chasing each other” around the population, one outcompeting the other only for itself to be outcompeted by the same type later. Instead, noncooperative populations are characterized by only occasional sweeps of beneficial mutations, making the issues of patch lifecycle and constrained patch size irrelevant. In other words, the effects of shape are contingent in the presence of cooperation and the within-population dynamics the cooperation creates.

What do our results tell us about the cooperation in the “real world”? After all, organisms in nature do not live on slender, bicycle tire shaped worlds, and although bulky doughnuts may host some microbes, they are not the most typical of habitats. However, even a brief review of microbial cooperative systems reveals much habitat diversity and highlights some strange and unexpected shapes. For example, filamentous Anabaena cyanobacteria inhabit a roughly linear, 1D world and exhibit cell differentiation with a stable ratio of vegetative (carbon fixing) and heterocystous (nitrogen fixing) cells (Kumar et al. 2010; Kim et al. 2011). The vegetative cells are effectively a germline and the nonreproducing heterocysts are much like soma, their differentiation caused by a molecular signal diffusing though neighboring cells (Yoon and Golden 1998). Although not a classical cooperation system, the division of labor in Anabaena is also evolutionarily unstable, where populations could be invaded and overtaken by “cheater” cells that forgo becoming heterocysts (Rossetti et al. 2010; Rossetti and Bagheri 2012). 2D and 3D bacterial worlds are certainly not difficult to find in nature and have routinely been studied in the laboratory. For example, one of the workhorses of microbial sociobiology, Pseudomonas aeruginosa, causes opportunistic infections via biofilm formation in the lungs of immunocompromised patients (Griffin et al. 2004; Köhler et al. 2009). The already complex structure of biofilms (Xavier et al. 2009), combined with potentially fractal shape of the lung itself (Nelson and Manchester 1988) may lead to bacterial populations with noninteger, fractal dimensions greater than 2 for their surface and even greater than 3 for the volume. Indeed, bacterial colonies have been found to grow in complex, fractal shapes, driven by both genetics and environment (Obert et al. 1990; Ben-Jacob et al. 1998). If some of those shapes select for or against cooperation, as our work suggests, then population shape is another avenue for evolution to fine-tune natural systems and for us to disrupt them in medically relevant situations (Brown et al. 2009).

In conclusion, cooperation is present in a wide variety of natural systems, ranging from microbes to humans, and it is inevitably instantiated in populations of individuals occupying a certain shape. In our work with Aevol and related digital platforms, we show that the population shape profoundly affects the evolution of cooperation and in some cases may differentiate between stable cooperation and no cooperation at all. To fully describe and predict the evolution of cooperation in all its diversity, together with classic factors such as population structure, future studies should also consider the effect of population shape.