How mutation alters the evolutionary dynamics of cooperation on network

Cooperation is ubiquitous in every level of living organisms. It is known that spatial (network) structure is a viable mechanism for cooperation to evolve. Until recently, it has been difficult to predict whether cooperation can evolve at a network (population) level. To address this problem, a previous study proposed a numerical metric, called Average Gradient of Selection (AGoS). AGoS can characterize and forecast the evolutionary fate of cooperation at a population level. However, stochastic mutation of strategies was not considered in the analysis of AGoS. Here we extend AGoS so that it can analyze the evolution of cooperation where mutation may occur to strategies of individuals on networks. We show that our extended AGoS correctly predicts the final states of cooperation with mutation in the individual-based simulations. Our analyses revealed that mutation always has a negative effect on the evolution of cooperation regardless of the payoff functions, fraction of cooperators, and network structures. Moreover, we found that scale-free networks are most vulnerable to mutation and thus the dynamics of cooperation is altered from bistability to coexistence on those networks, undergoing an imperfect pitchfork bifurcation.


Introduction
Cooperation is ubiquitous in every level of living organisms and has played an important role in the major evolutionary transitions [1][2][3]. In principle, cooperators benefit others by incurring some costs to themselves, while defectors do not pay any costs. Therefore, cooperation cannot be an evolutionarily stable strategy for a noniterative game in a well-mixed population [4][5][6][7][8].
In such a situation, spatial (network) structure is a viable mechanism for cooperation to evolve [9]. In particular, cooperation is significantly enhanced on scale-free networks with accumulated payoffs [10]. In this case, cooperative strategies can spread into the whole population because cooperative hubs surrounded by other cooperators obtain much higher payoffs compared to other nodes [10]. However, until recently, there has been no macroscopic metric that can be used to predict whether cooperation can evolve at a network (population) level, because of the complex interactions between evolution of strategies and topologies of networks (except for some limited cases, such as weak selection [11,12]). To address this problem, Pinheiro et al. [13] proposed a numerical metric, called Average Gradient of Selection (AGoS), to characterize and forecast the evolutionary fate of cooperation at a population level. AGoS can analyze the dynamics of the evolution of cooperation in structured populations even when nontrivial selection pressure is introduced [14], when nonlinear imitation probability is used [15], and also when structures and states of networks change over time (adaptive social networks) [16].
In these earlier studies, however, stochastic mutation of strategies was not considered (although Allen et al. mathematically analyzed the effect of mutation for cooperation in regular graphs under weak selection [17]). It is important to incorporate such mutation because they frequently occur in real societies and also because results obtained with stochastic fluctuations of strategies would provide more robust observations and conclusions. Here we analyze the evolution of cooperation using AGoS where mutation may occur to strategies of individuals in networks.

Models
We developed an agent-based model for the analysis of AGoS. Individuals are placed on the nodes in a network and they interact only with their neighbors. The networks used in this paper are described in the next subsection. Each individual can take one of two strategies: Cooperation (C) or Defection (D).
Once the composition of individuals in a network is given, we can calculate the probability of increasing or decreasing the number of cooperators by one, called Gradient of Selection (GoS) at time t [13]. Simultaneously, we can also update the strategies of individuals based on the framework of evolutionary games at time t. Right after the strategy updating, mutation (flipping strategy from C to D or vice versa) occurs with probability m. Therefore, in evolutionary simulations, the calculations of GoS and the strategy updating with mutation take place alternately in one time step and these processes are repeated. The calculation of GoS and the strategy updating with mutation are repeated for Λ time steps. Simulations are repeated for Ω times in total.

Network structure
Pinheiro et al. [13] revealed that, by using AGoS, cooperation was sustained in a network level and network structures led to the different evolutionary results of cooperation. Following the existing work [13,14], we focus on two classes of network structures: homogeneous and heterogeneous. We use homogeneous in the sense that every individual has the same degree. In the case of heterogeneous, individuals can have different degrees.
For homogeneous networks, we use two types: homogeneous random networks (HR) [18] and square lattice (Lattice). HR is created by randomizing links from homogeneous regular ring networks. For heterogeneous networks, we use Barabási-Albert scalfe-free networks (BA) [19]. In BA, a small number of nodes called hubs connect with a substantial number of links while most other nodes connect with a few nodes.
In the initial setting, the number of cooperators in a population, j ∈ [0, N ], is given to each simulation. The locations of j cooperators and N − j defectors in a network are random.

AGoS calculation
Once each simulation starts, GoS of the given population at time t is calculated. If there were no population (network) structure, the dynamics of cooperation could be calculated analytically using the replicator equation [4,8]. However, in structured populations, the calculation is difficult due to the complex topologies of networks (unless the network structure is regular or the selection is weak [11]).
The GoS gives the numerical characterization of the evolution of cooperation even in such a situation. Let φ i (t) be the probability that i's strategy changes to the other different strategy (from C to D or from D to C) and n i is a set of i's neighbors that have the strategy different from i's. Then, φ i (t) can be the product of two terms: the probability of selecting a neighbor with the different strategy, |ni| ki , and the average probability that i imitates the different strategy, 1 |ni| l∈ni p(i, l, t), where k i is the number of i's neighbors and p(i, l, t) is the probability that individual i imitates l's strategy at time t, which is defined in the next subsection. Thus, φ i (t) can be written as follows: From the definition, for a given simulation r, the probability to increase the number of C at j (the number of Cs in the network) is φ + r (j, t) = 1 N i∈sD φ i (t) while the probability for the decrease of the number of Cs is where N is the population size, s D is a set of defectors (|s D | = N − j), and s C is a set of cooperators (|s C | = j). Then, GoS at time t in a given simulation r is defined as the difference between them, as Finally, we obtain the Average Gradient of Selection (AGoS) [13], which averages the GoS by the number of time steps Λ and the number of simulations Ω, as The above derivation is the original AGoS without mutation [13]. Here we need to consider the effect of mutation on cooperation. To incorporate mutation, Eq. (1) needs to be revised so that the changed strategy by selection remains the same state with probability 1 − m while the unchanged strategy by selection changes to the opposite strategy with mutation probability m. Thus, the probability that i's strategy changes to the other different strategy is altered to Note that Eq. (4) returns to Eq. (1) when m = 0 (without mutation). Instead of Eq. (1), we use Eq. (4) for calculation of GoS (Eq. (2)) and AGoS (Eq. (3)). We use these extended versions for our analyses.

Strategy updating with mutation
To check the accuracy of the AGoS, we simultaneously conducted individual-based simulations. After the GoS calculation, the strategy updating takes place as follows. In each time step, we randomly choose one individual i from the population. The individual plays the Prisoner's Dilemma (PD) game with its k i neighbors and obtains the payoffs resulting from the games. In each game, both individuals obtain payoff R for mutual cooperation while P for mutual defection. If one selects cooperation while the other does defection, the former obtains the sucker's payoff S while the latter obtains the highest payoff T , the temptation to defect. The relationship of the four payoffs is usually T > R > P > S in PD games. AGoS is used for the other types of collective games including Stag Hunt [20] and Snowdrift [21] games.
Following the parameter settings used in the model by [13,14], we used P = 0, R = 1, and S = 1 − T , while T > 1. All the neighbors of i also play the PD game with their neighbors and obtains the payoffs. Let π i and π l be the payoffs of individual i and l (a randomly selected neighbor of i), respectively. Here we consider the two types of payoff functions: accumulated payoff and average payoff, because results may be quite different between those two conditions. In the accumulated payoff, cooperation is greatly promoted due to the hub effect [10]. In the average payoff [22][23][24][25][26](or in the case that having links takes some costs [27]), however, cooperation is inhibited because the hub effect is canceled. In the former case, the payoffs obtained from each game is accumulated over all the neighbors. In contrast, in the latter case, the accumulated payoffs for i at the end of the games is divided by his degree k i . In either case, the payoff is reset to zero at every time step.
Based on the framework of evolutionary games, higher fitness will be imitated more. In order to realize this, we use the pairwise comparison rule [28][29][30]. Individual i imitates l's strategy with the probability where β ≥ 0 controls the intensity of selection. Eq. (5) is also used in Eq. (1) for the calculation of AGoS. For β = 0, there is no selection pressure, meaning that evolutionary dynamics proceeds by random drift. As β becomes larger, the tendency that strategies with higher payoffs will be imitated increases.
In the previous studies [13,14], the next time step immediately follows after the strategy updating. In contrast, we incorporate mutation in this paper. Specifically, the strategy of individual i changes to the different strategy with probability m after the strategy updating. We focus on how mutation alters fitness of cooperation in networks.

Results
In each simulation, we calculated the GoS and conducted the strategy updating with mutation every time step, which was iterated Λ = 10 5 time steps. The population was composed of N = 1000 (except for Lattice for which N = 31 2 = 961) individuals and an average degree was k = 4. For each network type, we constructed 30 networks. For each network, we selected the number of initial cooperators from j = 0 to j = N randomly and this was repeated 1000 times. Thus, we ran Ω = 30 × 1000 = 30000 simulations in total for each network type and each mutation probability. We varied mutation probability m as the main experimental parameter. In the AGoS results, we obtained the locations of equilibrium points by numerically finding where the AGoS curve becomes zero. In doing so, the AGoS curves were smoothed using Gaussian kernels to reduce the effects of stochasticity in the results.

Case of accumulated payoffs
We first see the results of accumulated payoffs. Figure 1A shows the AGoS (G A (j)) on HR networks where the mutation probabilities m are varied. When m = 0, the result perfectly matches the corresponding case of Ref. [13]. It has two stable equilibrium points and two unstable equilibrium points. One of the unstable points, x L , exists at j/N = 0.043. One of the stable points, x R , exists at j/N = 0.578. The other unstable and stable points exist at j/N = 0 and j/N = 1, respectively. As the black arrows suggest, as long as j/N > x L , the population composition converges to the stable point x R . Thus, unlike the case of well-mixed populations where defectors are dominant, cooperation and defection can co-exist in the network. In other words, from a global, population-level perspective, HR networks can sustain cooperation even though all individuals play PD.
As m becomes larger, G A (j) goes down overall. Thus, the stable coexistence point x R becomes lower (see the red arrow in Fig. 1A). This means that mutation is harmful for the evolution of cooperation. Mutation itself is neutral because both strategy changes (cooperation to defection and vice versa) take place without any bias. Nevertheless, mutation negatively affects cooperation because it easily destroys the clusters of cooperators for cooperation to evolve.
At j/N = 0, G A (j) naturally rises when mutation is considered. At this point, G A (j) is the same as m, where there are only defectors and therefore any change from defector to cooperator must occur simply by mutation.
To verify the accuracy of our defined AGoS, we conducted the simulations of the strategy evolution ( Fig. 2A). Results of five hundred simulations are averaged for each. One million time steps are repeated for each simulation. We used two different initial fractions of cooperators, j/N (0) = 0.2 and j/N (0) = 0.8, to check those converging points. Without mutation (m = 0; left), the fraction of cooperators converges to 0.58 irrespective of the starting points. With mutation (m = 0.01; right), it converges to 0.41 irrespective of the starting points. We find that both converging points are nearly equal to the stable equilibrium points obtained through the analysis of AGoS. Thus, AGoS is appropriate to predict the fate of cooperation in HR networks.
We move to the case of lattice networks. Figure 1B shows the AGoS on lattice networks where the mutation probabilities m are varied. The tendency of this result is the same with the HR networks because lattice is classified into homogeneous networks (i.e., each node has the same degree). When m = 0, it has two stable equilibrium points and two unstable equilibrium points. One of the unstable points, x L , exists at j/N = 0.083. One of the stable points, x R , exists at j/N = 0.538. Hence, the density of cooperators in the coexistence situations is lower than HR networks (0.538 < 0.578). The other stable and unstable points exist at j/N = 0 and j/N = 1, respectively. If j/N > x L , cooperation and defection can coexist in the network. As m becomes larger, G A (j) goes down overall. Thus, the stable coexistence point x R becomes lower. Mutation negatively effects the evolution of cooperation in lattice networks, too. Figure 2B shows the simulations of the strategy evolution. The setting is the same with HR networks. We used two different initial fractions of cooperators, j/N (0) ≈ 0.2 and j/N (0) ≈ 0.8 (these are approximate values because N = 961 for Lattice), to check those converging point, respectively. Without mutation (m = 0; left), the fraction of cooperators converges to 0.55 irrespective of the starting points. With mutation (m = 0.01; right), it converges to 0.28 irrespective of the stating points. We find that both converging points are nearly equal to the stable equilibrium points in the analysis of AGoS. Thus, AGoS is appropriate to predict the fate of cooperation in lattice networks. Compared to HR networks every stable equilibrium point is lower in lattice networks. In lattice networks, cooperative clusters tend to be localized and separated from each other. It may prevent cooperators from expanding their regions in networks.
We next see the results of BA networks. Figure 1C shows the AGoS on BA networks where the mutation probabilities m are varied. When m = 0, there are one unstable equilibrium point x L = 0.402 and two stable equilibrium points, j/N = 0 and j/N = 1. This is the same with the one observed in the previous study [13]. This means that cooperation becomes dominant once j/N > x L . BA networks with accumulated payoffs changed the game from defector dominance to a bistable one. As m increases, the number of cooperators needed for sustaining the dominance of cooperation becomes larger. Interestingly, we found that, when m = 0.01, the dynamics completely change from bistability to coexistence. There exists only one stable equilibrium point at j/N = 0.142. In other words, a relatively small population of cooperators will coexist with the majority of defectors in this case. Figure 2C shows sample simulation results of the strategy evolution. The setting is the same with HR and Lattice networks. Without mutation (m = 0; left), the fraction of cooperators converges to 0.186 when j/N (0) = 0.2 while 0.974 when j/N (0) = 0.8. Thus, there are some deviations from the corresponding analytical values, 0 and 1, which are predicted by AGoS. In particular, the difference between 0.186 (see the red circle in the left panel of Fig. 2C) in the simulation and the value predicted by the AGoS (0). Similarly, with mutation (m = 0.01; right), both lines (j/N (0) = 0.2 and j/N (0) = 0.8) converge to around 0.475 in the simulations but the stable point predicted by AGoS is 0.142. Thus, there is also quite a difference between simulation and AGoS.
To clarify the reason, we examined individual simulation results. Figure 3 shows all 500 runs without mutation (m = 0) when j/N (0) = 0.2. From the results, we found that the fraction of cooperators did not converge to 0.186. Instead, four hundreds and seven simulations converged to 0 while ninety three simulations converged to 1 in 500 total runs. Thus, AGoS correctly predicted the final fate of cooperation with probability 0.814 but failed with probability 0.186. In BA networks, degree-strategy correlations significantly affect the increase and decrease of cooperation. However, the information cannot be captured by AGoS. This is why AGoS sometimes fails to predict the final fate of cooperation in BA networks.

Case of average payoffs
We now focus on how cooperation evolves with the average payoff compared to the accumulated payoff. It is known that, in networked games, cooperation is inhibited in the case of average payoffs because the hub effect is dismissed [22][23][24][25][26]. Figure 4 show the case of average payoffs. In HR networks (Fig. 4A), the internal stable equilibrium point is 0.598 without mutation (m = 0) and is 0.437 with mutation (m = 0.01). The corresponding simulations shown in Fig. 5A are 0.6 without mutation (m = 0; left) and are 0.44 (m = 0.01; right), which means that the prediction by AGoS is quite accurate. In Lattice networks (Fig. 4B), the internal stable equilibrium point is 0.560 without mutation (m = 0) and is 0.303 with mutation (m = 0.01). The corresponding simulations shown in Fig. 5B are 0.55 without mutation (m = 0; left) and are 0.28 (m = 0.01; right). The AGoS predicts the results of simulations well in Lattice networks as well as the HR networks.
For those two networks, if we compare the results of AGoS between accumulated ( Fig. 1A and B) and average ( Fig. 4A and B) payoffs, both equilibrium points in the average payoff are slightly higher than those in the accumulated payoff (see also Fig. 7). There is not so much difference between the average and accumulate payoffs because each node has the same degree in those networks. However, if the payoffs are averaged, the payoff difference between cooperators and defectors shrinks, hence the chance that cooperators can survive slightly increases in the average payoffs.
We move to the case of BA networks. Figure 4C shows the AGoS on BA networks where the mutation Compared to the case of accumulated payoff (Fig. 1C), the critical point x L greatly moves to the rightward, which means that the critical fraction of cooperators needed for cooperation to become dominant is quite large. This is because the hub effect is dismissed in the case of average payoffs. The corresponding simulations (Fig. 5C) show that AGoS predicts well when j/N (0) = 0.2 where the simulations converge to 0 but largely fails when j/N (0) = 0.8 where the simulations converge to around 0.2. The reason is that degree-strategy correlations are missed in AGoS in BA networks.
When m = 0.01, the dynamics changes from bistability to coexistence in the same way as with the accumulated payoffs. The stable equilibrium point is 0.372 (Fig. 4C) and the fraction of cooperators converge to 0.39 in the corresponding simulations (right in Fig. 5C). We consider the possible reason for the change of this dynamics in the next subsection.

Mutation effect on BA networks
We saw that the dynamics completely changed on BA networks in both accumulated and average payoffs. To consider which factor contribute to the change of dynamics, we show a case in which only mutation occurs with probability 0.01 (m = 0.01) without strategy updating (red line in Fig. 6). If there is no cooperator (j/N = 0), the expected probability that cooperator increases by one is 0.01 (G A (j) = 0.01). In contrast, if there is no defector (j/N = 1), the probability that cooperator decreases by one is 0.01 (G A (j) = −0.01). If cooperators and defectors are fifty-fifty, that probability becomes 0. Note that, if there is no strategy updating, these results always hold irrespective of network structures because the locations of individuals and their payoffs do not matter.
We first see the comparison of the accumulated payoffs with the case of only mutation. In Fig. 6, it seems that the dynamics are heavily influenced by mutation because the solid orange line curves in accordance with the red line. In other words, the dynamics changed from bistability to coexistence due to the effect of mutation. In BA networks, contrary to HR or Lattice networks, once a cooperative hub is destroyed by mutation, there is almost no chance for cooperation to evolve. This is why BA networks are most vulnerable to mutation compared to HR or Lattice networks.
As in the case of accumulated payoffs, the change of dynamics in the average payoffs arises from the effect of mutation. As seen in Fig. 6, the AGoS in the average payoff (dashed line) is closer to the case where only mutation is considered because there is no hub effect in this case. Thus, the dynamics in the average payoffs is highly influenced by mutation compared to the case of the accumulated payoffs.

Comparison of equilibrium points
Finally, we take a close look at the location of stable equilibrium points against mutation probability on the three networks. Here, m is varied from 0 to 0.1. Figure 7 shows the results in (A) HR, (B) Lattice, and (C) BA networks. In each network, the solid lines indicate the stable equilibrium points and the dashed lines indicate the unstable equilibrium points. Blue (Red) is the case of the accumulated (average) payoffs.
By comparing (A) with (B), we find that the stable equilibrium points in HR networks are always higher than those in Lattice networks. This is because cooperators are often locally clustered in Lattice networks compared to HR networks as described before. In those two networks, as m becomes larger, the stable equilibrium points go down unless m is too large (e.g. m = 0.1). The reason is that cooperative clusters are destroyed by mutation. In contrast, if m is too large, cooperators are continuously produced by mutation irrespective of the strategy updating. It contributes to the rise of equilibrium points.
In BA networks, as we described before, the dynamics changes from bistability to coexistence. Moreover, by changing m in more detail, we find the imperfect pitchfork bifurcations in the accumulated payoffs (Fig. 7). The main reason of the change to coexistence was that scale-free networks are vulnerable to mutation. Once cooperative hubs (only key for cooperation to evolve in scale-free networks) are destroyed by mutation, the dynamics is mostly governed by mutation, which lead to the coexistence. Because there is no hub effect in the average payoffs, the coexistence takes place at smaller m values.

Discussion and conclusions
In this paper, we analyzed the evolution of cooperation at a population level in the presence of mutation by the AGoS. We used two classes of networks: homogeneous (HR and Lattice) and heterogeneous (BA). We also examined the two different payoff mechanisms: accumulated and average payoffs. Our analyses revealed that mutation has a negative effect on the evolution of cooperation regardless of the payoff functions, fraction of cooperators, and network structures, because clusters of cooperators can easily be destroyed by mutation. Moreover, we found that scale-free networks are most vulnerable to mutation because a few cooperative hubs are the only key to the success for cooperation despite that cooperation is greatly enhanced with accumulated payoffs without mutation in those networks. If a few hubs are destroyed by mutation, cooperation easily collapses. This may be relevant to the robustness of scale-free networks. In general, scale-free networks are quite weak when hubs are destroyed by targeted attacks [31]. Moreover, cooperation on scale-free networks easily collapses when hubs are removed [32,33]. Our paper shows that, even if hubs are not directly removed, scale-free networks are vulnerable to mutation because hubs are easily invaded by defectors through mutation. Therefore, to prevent cooperative societies from collapsing, it may be needed to develop some mechanisms how hubs are protected from mutation.
Our results indicate the importance of considering random noise (mutation), which was largely overlooked in the literature, in studying the evolution of cooperative behavior in social networks. Mutation can be considered genetic changes if we assume biological systems. If we assume cultural systems, mutation can be considered a stochastic behavior to explore new behaviors. Although we showed such a random exploration was harmful to keep cooperative societies, if we consider different forms of exploration, those explorations may work beneficially for societies, such as collaborative problem solving [34]. We would like to consider those cases in the future work.
Finally, we remark the limitation of this paper. Our extended AGoS considering mutation sometimes fails to predict the final states of cooperation on BA networks although it perfectly works for HR and Lattice networks. This was because AGoS aggregate every state and topology of networks into one value. To consider those complicated network interactions, we may need to take into account the influence of nodes separated by the degree. If the new version will be developed in the future, it may offer better predictions even if networks are heterogenetic. In this paper, we analyzed the evolution of cooperation at a population level in the presence of mutation by the AGoS. We used two classes of networks: homogeneous (HR and Lattice) and heterogeneous (BA). We also examined the two different payoff mechanisms: accumulated and average payoffs. Our analyses revealed that mutation has a negative effect on the evolution of cooperation regardless of the payoff functions, fraction of cooperators, and network structures, because clusters of cooperators can easily be destroyed by mutation. Moreover, we found that scale-free networks are most vulnerable to mutation because a few cooperative hubs are the only key to the success for cooperation despite that cooperation is greatly enhanced with accumulated payoffs without mutation in those networks. If a few hubs are destroyed by mutation, cooperation easily collapses. This may be relevant to the robustness of scale-free networks. In general, scale-free networks are quite weak when hubs are destroyed by targeted attacks [31]. Moreover, cooperation on scale-free networks easily collapses when hubs are removed [32,33]. Our paper shows that, even if hubs are not  directly removed, scale-free networks are vulnerable to mutation because hubs are easily invaded by defectors through mutation. Therefore, to prevent cooperative societies from collapsing, it may be needed to develop some mechanisms how hubs are protected from mutation.
Our results indicate the importance of considering random noise (mutation), which was largely overlooked in the literature, in studying the evolution of cooperative behavior in social networks. Mutation can be considered genetic changes if we assume biological systems. If we assume cultural systems, mutation can be considered a stochastic behavior to explore new behaviors. Although we showed such a random exploration was harmful to keep cooperative societies, if we consider different forms of exploration, those explorations may work beneficially for societies, such as collaborative problem solving [34]. We would like to consider those cases in the future work.
Finally, we remark the limitation of this paper. Our extended AGoS considering mutation sometimes fails to predict the final states of cooperation on BA networks although it perfectly works for HR and Lattice networks. This was because AGoS aggregate every state and topology of networks into one value. To consider those complicated network interactions, we may need to take into account the influence of nodes separated by the degree. If the new version will be developed in the future, it may offer better predictions even if networks are heterogenetic.

Competing interests
We have no competing interests.
Author's contributions GI and HS designed the research. GI and YS developed the model. GI and YS performed the simulation. GI and HS discussed and analyzed the results. GI and HS wrote the main manuscript text. All authors gave final approval for publication.