Threshold for Cooperation on Irregular Spatial Networks

Cooperation is a defining attribute of life as we know it, from the delicate interactions of intracellular components to social behavior in groups. However, defection and exploitation are at least as ubiquitous. Evolutionary game theory is a successful tool for investigating how cooperation may be maintained despite large advantages for defection. The Prisoners Dilemma is one such game where spatial structure can maintain cooperation, but only if the benefit-to-cost ratio (b/c) is greater than some threshold, which appears to be the average number of neighbors (k). However, this inequality was tested only for regular spatial and irregular non-spatial networks. In this paper, we use networks in Cartesian space that are based on radii of interactions. We investigate whether the b/c > k threshold holds for these irregular spatial networks, and we use a much broader range of k than previously studied. We find that this rule, and other related inequalities, hold well for the larger radii even when there is noise in the expected neighborhood size. As the expected neighborhood size increases, so does the variation in the empirical edge distribution. However, the variation in the threshold for cooperation decreases. This paper is a first step in a broad investigation of how uncertainty affects the outcome of game theoretic simulations.


Introduction
Cooperation is a widely observed phenomenon in biology.Whether this involves complex societies in ants (Shik et al., 2012), group hunting in wolves (Escobedo et al., 2014), cooperation among cells in a multicellular individual (Aktipis et al., 2015), cooperation between organelles within a eukaryotic cell (Rand et al., 2004), or cooperation among bacterial cells (Brockhurst et al., 2010), a hallmark of cooperation is that individuals pay some form of cost by engaging in cooperative behavior.Defectors exploit this cooperation, gaining the benefits produced by the cooperators without paying the costs (Nowak, 2006).As such, the near ubiquity of cooperation within biological systems has long been a puzzle in evolutionary biology.
One of the most widely used tools to study the evolution of cooperation is game theory, a branch of mathematical science that addresses situations where the results of an action on the part of an agent depend on the actions of other agents in the system.A classic model in game theory is the Prisoners Dilemma (PD).In this simple game, an agent and its partner must each decide whether to cooperate or defect (Table 1).Over the course of many interactions, partnerships where both individuals cooperate can be highly advantageous; however, in any single instance, an individual will do best if it defects, regardless of whether or not its partner cooperates.How, then, can systems reach the beneficial point of mutual cooperation, when there is an immediate benefit to defection?
A variety of mechanisms have been proposed to explain cooperation.Some mechanisms, such as policing cheaters (Bausch, 2014), assessing reputation for cooperation (Cuesta et al., 2015;Romano et al., 2017), or even simple mutual choice to interact (Nguyen et al., 2016) require the ability to recognize specific individuals.Other mechanisms, such as kin selection (smith et al., 2016;Birch and Okasha, 2015) and spatial structure (reviewed in (Perc et al., 2013)) can operate even in the absence of advanced sensory or cognitive abilities on the part of the agents.As such, in this work we focus on the role of spatial structure in cooperation.Spatial structure is broadly observed in biological systems, even microbial ones (Green and Bohannan, 2006;Harcombe et al., 2014), and thus is widely applicable.
Previous work in PD has demonstrated that there exists a critical threshold ratio of benefits to costs in order for cooperation to be favored (Nowak, 2006;Szabó et al., 2005;Wahl and Nowak, 1999;Szolnoki and Perc, 2010).Within spatial networks, this threshold is proportional to k, the average degree of node of the network (i.e. the average number of neighbors with which an individual interacts).Researchers have derived this threshold for small values of k by using pair approximation (Matsuda et al., 1992), but these predictions assume k ≤ 10.While the number of individuals a given individual may interact with is hard to quantify in many biological systems, it is reasonable to assume that this will often exceed 10.In this work, we vary neighborhood size over a range of values not previously considered in analytical work for irregular spatial networks.We note that some excellent simulations have been done with a vari-Figure 1: Order parameter over time for a set of runs.We selected a radius that produces an expected 15 neighbor, and plot here 10 runs at range of u values that span from asymptotic behavior of full cooperation to full defection.We present the data for (A) the first 5000 generations and (B) the first 500 generations of these runs; note that the data in (B) is a subset of the data in (A).A more colorblind-friendly version of this figure is available in github.com/devosoft/SpatialPD/tree/master/ECAL2017ety of network structures, in which the authors have focused on other elements of the networks than we are examining here, such as on modularity and motifs (e.g.(Gianetto and Heydari, 2016)), or have shown that irregularity of networks may promote the evolution of cooperation without specifically considering its effect on the critical threshold of cooperation cost-to-benefits (e.g.(Pinheiro et al., 2012b;Santos et al., 2006;Pinheiro et al., 2012a)).In contrast, in our study, we examine the threshold criteria for evolution of cooperation across a variety of neighborhood sizes (k) and compare our empirical results to previous theoretical predictions.We also examine the effects of unpredictability in the number of neighbors an individual has on the threshold for cooperation.

Theoretical Predictions
Within spatial networks, cooperation will be favored if the ratio of benefits to costs exceeds a critical threshold.This threshold depends on the specific details of how the game is implemented.In Death-Birth methods, where an individual is chosen to die and then another individual is able to produce an offspring, Equation 1 holds for pair approximation (Matsuda et al., 1992).With Imitation methods, where an individual is chosen to update its strategy to that of one of its neighbors (Ohtsuki et al., 2006), Equation 2 accounts for the additional edges in the network.And when the ability of a small group of cooperators to jointly invade a population of defectors is considered (van Baalen and Rand, 1998), Equation 3 describes the estimated threshold.These theoretical expectations are derived explicitly for small values of k.We are interested in whether they apply to larger values of k, in irregular spatial networks.

Methods
To examine the transition between cooperation and defection, we built a computational model based on spatial networks (both the project code and analysis script are available at github.com/devosoft/SpatialPD/tree/master/ECAL2017).
Carole Knibbe et al, eds., Proceedings of the ECAL 2017, Lyon, France, 4-8 September 2017, (Cambridge, MA: The MIT Press, ©2017 Massachusetts Institute of Technology).This work is licensed to the public under a Creative Commons Attribution -NonCommercial -NoDerivatives 4.0 license (international): http://creativecommons.org/licenses/by-nc-nd/4.0/ In this model, we generate irregular spatial networks by randomly placing nodes in Cartesian space (with a uniform distribution) and connecting each pair of nodes that are within a specified distance, wrapping edges so as to avoid boundary conditions.We use Equation 4 to determine the radius for which we expect a particular number of neighbors, where k is the number of neighbors, N is the number of nodes on the entire network, and r is the radius of interaction.Using this method, we are able to vary the k parameter with fine granularity, and are no longer constrained by the geometric growth of neighborhoods on lattices.
For each treatment, we run 100 replicate populations for 5000 generations, which is long enough to reach their asymptotic limits of behavioral changes and so we do not have to rely on long-term projections from small initial differences in growth rate (see Figure 1).We seed populations in this system with 6400 nodes, each of which contains a single individual.Each individual plays one of two strategies cooperator or defector randomly assigned at the beginning of the experiment.We update these populations asynchronously.That is, during each generation, the system randomly chooses 6400 nodes with replacement.As each node is chosen, the strategy of the individual contained in the node is updated probabilistically, based on the total payoff received by its neighbors (the full rules are described below).This way, each individual is updated on average once per generation.The strategies in this simulation are fixed, and have no stochasticity or memory of past interactions.Thus, a cooperator is always a cooperator, and a defector is always a defector.
Individuals in this model compete directly with their immediate neighbors.As in other game-theoretic models, the outcomes of these competitions depend on the strategies played by each organism (Doebeli and Hauert, 2005).A cooperator provides a benefit to its opponent while incurring some cost.Defectors provide no benefit and incur no costs.In order to capture the effects of these interactions using one parameter, we define payoffs in terms of a cost-benefit ratio u = c/b as in (Langer et al., 2008).For the PD game, this setup results in the payoff matrix shown in Table 2, which indicates the payoff for the row player based on each combination of the row player and column players actions.As a cell is updated (and therefore replaced), that cell is compared against all of its neighbors.A single neighbor cell is chosen with probability proportional to the total payoffs it has accrued.The strategy of the winning neighbor is then copied into the focal cell, a process similar to previous studies (e.g.(Fu et al., 2010)).An individuals interactions are limited to its immediate neighbors in space; that is, just those within the interaction radius.A populations interactions are modeled using the irregular spatial networks described above, where each cell is represented as a vertex (Cartesian networks).
For this work, we vary the cost/benefit ratio (u) and the average neighborhood size (k) of the networks in order to determine how these factors affect the evolution of cooperation.We chose values for u to encompass the area of parameter space surrounding the area where the transition between cooperation and defection is predicted to occur based on the average per-node degree of the network as described by equation 5 in (Ohtsuki et al., 2006).
In order to characterize the asymptotic behavior of our simulations, we calculate an order parameter based on the number of cooperators and defectors present at the end of the run, as described in equation 5: N c is the number of cooperators in the population, and N d is the number of defectors, we effectively scale the extreme outcomes of all cooperators and all defectors to 1 and -1, respectively.Plotting the order parameter as it varies with u lets us visually inspect the critical u value for the evolution of cooperation as the point where the order parameter crosses zero (Iliopoulos et al., 2010).

Statistical Methods
For all of our networks, we considered only those nodes with at least one neighbor when we calculated average k or average predicted u in that network, as nodes with no neighbors can never change state in our model.We calculated non-parametric 95% Confidence Intervals around theoretical prediction thresholds by taking the average k parameter for each network involved in the analysis, calculating the threshold for that particular average k, and retaining the central 95% of all such calculated thresholds.We performed all data analysis in R 3.2.13(R Core Team, 2015).

Results and Discussion
Harnessing the ability to continuously vary average neighborhood size, we investigated various related threshold criteria for the evolution of cooperation on spatial networks over a much broader range than previous studies.The simulations run long enough to see asymptotic behavior between strategies, instead of simply relying on initial growth rates or strategy expansion when rare.Looking at asymptotic  In order to gain a better intuition about the dynamics within a run, we built a visualization tool for our program, accessible at http://www.cse.msu.edu/∼ofria/SpatialPD/ Figure 2 shows a time series of images from a single evolutionary run on one such network, with parameter settings such that we expect cooperators to outcompete defectors.At the start of any individual run, cooperators and defectors are well mixed throughout the world.Even a single generation into the run, clusters begin to form; groups of high cooperator density are able to gain substantial benefits from mutual while cooperators at lower density are vulnerable to invasion by defectors that exploit the benefits of cooperation without paying the costs.
Since we re-parameterize the payoff matrices based on the u parameter, we can also express the cooperation critical point estimates in terms of this re-parameterization, from Equations 1-3 we have: where u = c/b.As k increases, these estimates converge.Figure 3 shows these three critical u value estimates as a function of k.Interestingly, they all converge quickly for Carole Knibbe et al, eds., Proceedings of the ECAL 2017, Lyon, France, 4-8 September 2017, (Cambridge, MA: The MIT Press, ©2017 Massachusetts Institute of Technology).This work is licensed to the public under a Creative Commons Attribution -NonCommercial -NoDerivatives 4.0 license (international): http://creativecommons.org/licenses/by-nc-nd/4.0/ even moderate values of k; that is, they estimate essentially the same critical u value as each other by the time k reaches 15.This result suggests that the critical u value is quite robust despite different update rules, fixation measures, and other assumptions for reasonable neighborhood sizes.
Figure 3: Predicted critical u thresholds as a function of neighborhood size (k).This visually represents equations 6-8, and their convergence of estimates This robustness is also evident in Figure 4, where the final order parameter values are reported over the experimentally manipulated u values.Where these curves cross zero represents the empirical thresholds for cooperation.The three vertical lines are drawn where the estimated threshold value fall.In each case, the largest estimate is for Equation 6(black); the middle estimate is for Equation 8 (red), and the smallest estimate is for Equation 7 (light blue).
As we expect, there is more divergence between estimates at low values of k, as is evident in Figure 4 (a) vs (d), where with k = 8 the three estimates are distinct, and with k = 24 they become strikingly similar.Thus, all three derivations of the critical u value perform well, as they typically flank the empirical estimates for the thresholds of cooperation.Since each estimate is derived under slightly different conditions, and potentially is used to explain different definitions of invasibility, we do not attempt to argue which performs best.Rather, we show how all predicted thresholds fit reasonably well with the asymptotic critical value of u we measure empirically, and that these values converge as the k parameter value increases.
When initially approaching this work, we expected that increased neighborhood sizes (k) would lead to an increased variation in the estimated cooperation threshold.Previous studies using this same method of generating irregular spa-Figure 4: Asymptotic behavior for different neighborhood sizes across configured values of u.The horizontal line is drawn through y = 0 to aid visual interpretation.Vertical lines are the estimated thresholds for cooperation drawn with the same legend as Figure 3.Note that for each panel, we set the r parameter to produce the expected neighborhood size according to Equation 4 and label the panels based on these expectations; we shade the non-parametric 95% CI for each threshold based on the average k observed in each network.All three estimates predict the empirical threshold of cooperation well, and clearly converge as the neighborhood size increases.tial networks reported an increase in variation of per-node edges as k increased (Connelly et al., 2010).Figure 5 (a) shows how the standard deviation of the per-node edge distribution increases as a function of k.We term this sort of variation social noise.There may be high amounts of social noise when players can potentially interact with a wide distribution of neighbors at any given time; conversely, there is no social noise in traditional lattice interaction networks.We expected that this increase in noise around the value of k would lead to more noise in the empirical transition between cooperation and defection.However, Figure 5 (b) shows a decrease in variation of the per-node critical u values despite the increased amount of variability in k.This reduction in variation is due to the critical u value being estimated as 1/k, and thus the increasing k value drives the estimated critical u lower, and any noise in k is overwhelmed by the exponential decrease in the u estimate.Additionally, while the magnitude of variation in k increases with increasing k, it does so at a sublinear rate.Figure 6

Conclusions
We have explored the interplay between the cost-to-benefit ratio and the number of interactions an individual has on the evolution of cooperation.This project builds upon previous works that were limited to small neighborhood sizes and regular spatial or non-spatial networks.Specifically, we ex- Though the absolute level of variation in k increases with k, it does so slower than the rate of increase in k, and thus has a relatively smaller impact.plored threshold estimates using an irregular spatial network based on radii of interactions in Cartesian space, as well as over larger neighborhood sizes than originally investigated.We found that all three estimates perform well at predicting the critical cost-to-benefit ratio (u) where the populations transition between all cooperators and all defectors.Despite an increasing amount of variation in k as neighborhood size increased, the noise associated with the critical u value decreased due to their functional relationship.
In order to experimentally test if the noise in neighborhood size actually has a significant effect on the evolution of cooperation, we would ideally manipulate the standard deviation in the distribution of neighborhood sizes independently of the expected value k.We have begun to explore this avenue, but only preliminary data have been collected.In the future, we plan to explore how this and other kinds of noise in spatial game theory may interact to promote or inhibit the evolution of cooperation. b for the Prisoner's Dilemma.C represents cooperation; D represents defection.b is the benefit of cooperation; c is the cost of cooperation.Payoffs are shown for the row player.
Table 2: Payoff matrix for the Prisoner's Dilemma.C represents cooperation; D represents defection.u = c/b where b is the benefit of cooperation; c is the cost of cooperation.Payoffs are shown for the row player.

Figure 2 :
Figure 2: Visualization of one example run, with settings r = 0.02, u = 0.1, N = 6400.Each panel shows a different snapshot of one network.Blue circles represent cooperators; red circles represent defectors.As time progresses, clusters form of each type.Different panels show the population after different numbers of generations.Panel A) Initial population.Panel B) 1 generation.Panel C) 50 generations.Panel D) 100 generations.Panel E) 500 generations.Panel F) 2000 generations.
shows how the coefficient of variation in k as a function of increasing k; the amount of variation in the k parameter relative to its average value Carole Knibbe et al, eds., Proceedings of the ECAL 2017, Lyon, France, 4-8 September 2017, (Cambridge, MA: The MIT Press, ©2017 Massachusetts Institute of Technology).This work is licensed to the public under a Creative Commons Attribution -NonCommercial -NoDerivatives 4.0 license (international): http://creativecommons.org/licenses/by-nc-nd/4.0/

Figure 5 :
Figure 5: Relationship between variation in neighborhood size and the empirical threshold for cooperation.The values here are the standard deviation in the per-node empirical neighborhood size (A) and an estimate of the critical u parameter (B) for each neighborhood size.Despite the increasing variation in neighborhood size, the variation in the threshold for cooperation decreases.

Figure 6 :
Figure 6: Coefficient of variation in k as a function of k.Though the absolute level of variation in k increases with k, it does so slower than the rate of increase in k, and thus has a relatively smaller impact.