Perturbation theory for evolution of cooperation on networks

Network structure is a mechanism for promoting cooperation in social dilemma games. In the present study, we explore graph surgery, i.e., to slightly perturb the given network, towards a network that better fosters cooperation. To this end, we develop a perturbation theory to assess the change in the propensity of cooperation when we add or remove a single edge to/from the given network. Our perturbation theory is for a previously proposed random-walk-based theory that provides the threshold benefit-to-cost ratio, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(b/c)^*$$\end{document}(b/c)∗, which is the value of the benefit-to-cost ratio in the donation game above which the cooperator is more likely to fixate than in a control case, for any finite networks. We find that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(b/c)^*$$\end{document}(b/c)∗ decreases when we remove a single edge in a majority of cases and that our perturbation theory captures at a reasonable accuracy which edge removal makes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(b/c)^*$$\end{document}(b/c)∗ small to facilitate cooperation. In contrast, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(b/c)^*$$\end{document}(b/c)∗ tends to increase when we add an edge, and the perturbation theory is not good at predicting the edge addition that changes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(b/c)^*$$\end{document}(b/c)∗ by a large amount. Our perturbation theory significantly reduces the computational complexity for calculating the outcome of graph surgery.


Introduction
Since Darwin's time, explaining cooperative behavior in groups of self-interested individuals has been a challenge [1][2][3][4][5][6][7][8].Game theory including evolutionary game theory has shown that a population of self-interested individuals playing a social dilemma game of the prisoner's dilemma type does not sustain cooperation without an additional mechanism.To explain cooperation in social dilemma situations in nature including in biological populations and to promote cooperation in human society, there have been proposed various mathematical mechanisms to support cooperation.Population structure as represented by contact networks of individuals is one such mechanism.The structure of contact networks constrains who can interact with whom and promotes emergence and endurance of clusters of cooperative players in local regions in spatial lattices [2,[9][10][11] and adjacent pairs of nodes in general networks [11][12][13][14][15].
A major indicator of the success of a mutant trait in evolutionary dynamics is the fixation probability.It is defined as the probability that the mutant type will spread and eventually occupy the entire population as a result of evolutionary dynamics, given an initial distribution of mutants [5,[16][17][18].When each individual is in either of the two types (i.e., wild and mutant) at any given time and the population structure is described by a network on N nodes, the state of the network is specified by an N -dimensional binary vector of which the ith entry encodes the type of the ith node.In the absence of mutation, the fixation probability of the mutant starting from the state in which all the nodes are of the wild type is equal to 0. The fixation probability of the mutant is equal to 1 if all the nodes are initially mutant.For general initial conditions, the exact solution of the fixation probability requires solving a linear system of 2 N − 2 equations [13,18].Therefore, it is difficult to exactly compute the fixation probability except for small networks, highly symmetric networks, or networks with other mathematically convenient properties.
We focus on social dilemma situations, in particular the prisoner's dilemma game, in the present paper.In the prisoner's dilemma, the wild and mutant types correspond to cooperator and defector, respectively, or vice versa.The calculation of the fixation probability for the prisoner's dilemma game on networks, potentially with some additional assumptions, is usually more involved than the calculation in the case of the constant selection, in which the fitness of the wild and mutant types is fixed throughout the evolutionary dynamics.In games, the fitness of an individual generally depends on how other individuals behave, which makes setting up the linear system of 2 N − 2 equations and efficiently solving it, particularly the latter, a difficult task.Under this circumstance, weak selection is an assumption that often facilitates analytical evaluation of the fixation probability of the mutant type including in social dilemma games [16].Let us write down each individual's fitness as a sum of a constant term, called the baseline fitness, and the payoff that the individual receives by playing the game.By definition, weak selection means that the payoff is small compared to the baseline fitness.Under weak selection, Ohtsuki et al. developed a pair approximation theory that enables us to analytically derive the conditions under which cooperation fixates with a larger probability than a baseline on random regular graphs, i.e., random graphs in which all nodes have the same number of neighbors [13].Furthermore, Allen et al. extended this result to the case of arbitrary networks using coalescence times from random walk theory [15].With these methods, one can avoid dealing with a set of 2 N − 2 linear equations and calculate the leading term of the fixation probability in polynomial time in terms of N .
In Ref. [15], the authors derived a key indicator to quantify the ease of cooperation in networks, i.e., the threshold benefit-to-cost ratio above which selection favors cooperation, denoted by (b/c) * .In fact, substantial changes in (b/c) * may occur when one only slightly perturbs the network structure, which is an operation referred to as graph surgery [15].A carefully designed graph surgery may enhance cooperation by reducing (b/c) * by a larger amount than by a random graph surgery.For example, a small mean degree (i.e., the number of neighbors that a node has) of the network tends to induce cooperation [13,15].Therefore, decreasing the weight of an edge or removing an edge is expected to enhance cooperation.However, this may not be an optimal choice.Which particular edge should we perturb or remove to efficiently enhance cooperation?One can answer this question by removing just one edge from the original network, calculating (b/c) * for the perturbed network, and repeating the same procedure for each different perturbation of the original network.However, this procedure may be computationally costly.Note that the method to calculate the fixation probability for cooperation in arbitrary networks, developed in Ref. [15], is still computationally costly although its computational complexity is polynomial in N .
In the current study, we develop a perturbation theory with the aim of predicting the direction and amount of the change in (b/c) * when one slightly perturbs the weight of an arbitrary single edge.We find that, for most networks, the actual change in (b/c) * when we remove an edge and the change predicted by our perturbation theory are strongly correlated, which makes it possible to propose a single edge to be removed for efficiently enhancing cooperation.However, the correlation between the result of direct numerical simulations and the perturbation theory is considerably weaker when one adds a new edge to the existing network.Therefore, our perturbation theory is not practically useful when one adds new edges.Compared to the direct numerical simulations, our perturbation theory is much faster, which allows us to compute the fixation probability under graph surgery in larger networks.

Fixation of cooperation on networks under weak selection
We assume that the graph G is connected and undirected.We denote the set of nodes by V = {1, . . ., N }, where N is the number of nodes.For each pair of nodes i, j ∈ V , we denote the edge weight by w ij ≥ 0. If there is no edge between i and j, we set w ij = 0. We allow self-loops, i.e., positive values of w ii [15].The weighted degree of node i, denoted by s i = N j=1 w ij , also called the node strength, is the sum of the weight of the edges connected to the node.
A discrete-time random walker is said to be simple if the walker located at node i moves to one of its neighbors, denoted by j, in a single time step with probability proportional to w ij , i.e., with probability p ij = w ij /s i .Let W = (w ij ) be the N × N weighted adjacency matrix.The transition probability matrix P = (p ij ) of the simple random walk is given by P = D −1 W , where D = diag(s 1 , . . ., s N ), i.e., the diagonal matrix whose diagonal entries are equal to s 1 , s 2 , . . ., s N .Let π = (π 1 , . . ., π N ) be the stationary probability vector of the random walk with transition probability matrix P , i.e., the solution of πP = π.It holds true that [19,20] We use the donation game, which is a special case of the prisoner's dilemma game.In the donation game, which is a two-player game, one player, called the donor, decides whether or not to pay a cost c (> 0).If the donor pays c, which we refer to as cooperation, then the other player, called the recipient, receives benefit b (> c).If the donor does not pay c, which we refer to as defection, then the donor does not lose anything, and the recipient does not gain anything.Therefore, the payoff matrix of the donation game for a pair of players is given by where C and D represent cooperation and defection, respectively, and the payoff values represent those for the row player.We assume that each player on a node participates in the game as donor and recipient half of the times each.
We assign 0 and 1 to the defector and cooperator, respectively.Then, we can represent a state of the entire network by a binary vector x = (x 1 , . . ., x N ) ∈ {0, 1} N .With this notation, the payoff of node i averaged over all its neighbors is given by The reproductive rate of node i in state x is given by where η represents the strength of the selection.If η = 0, the reproductive rate does not depend on the payoff matrix or the action (i.e., cooperation or defection) of any node.This case is equivalent to the so-called voter model.If η → 0, the payoff weakly impacts the selection, and this limit is called the weak selection regime.The idea behind weak selection is that, in reality, many different factors may contribute to the overall fitness of an individual, and the game under consideration is just one such factor [13,15].We drive evolutionary dynamics by the death-birth process with selection on birth on an arbitrary network composed of cooperators and defectors [13,15].Specifically, we first select a node to be updated, denoted by i, uniformly at random.Second, we select one of the i's neighbors, denoted by j, for reproduction with the probability proportional to w ij R j (x).Third, the offspring, i, inherits the type of j.This completes a single round of the evolutionary dynamics, which we schematically show in Fig. 1.
The death-birth process in any finite population without mutation will eventually reach the state in which all individuals are cooperators or defectors and halt.In other words, the cooperation or defection fixates in finite time with probability 1. Suppose the initial condition in which one node is cooperator and the other N − 1 nodes are defectors.There are N such initial conditions depending on which node is the cooperator.We consider the initial probability distribution over all possible states that assigns probability 1/N to each of the states with exactly one cooperator and probability zero to all the other states.We denote by ρ C the expectation that the cooperation fixates under this distribution of the initial state.If ρ C > 1/N , natural selection favors cooperation [5,13,15,16].In Ref. [15], Allen et al. showed that where ij is the (i, j)th entry of matrix P k , which implies that p (1) ij = p ij , and Equation (7) implies that t ij = t ji is the mean coalescence time of two random walkers when one walker is initially located at node i and the other at node j.Note that p (k) ij is the k-step transition probability of the random walk from node i to node j.Therefore, τ k is the expected value of t ij when i and j are the two ends of a k-step random walk trajectory on G under the stationary distribution [15].Equation (5) implies that the threshold value of the benefit-to-cost ratio above which the natural selection favors cooperation (i.e., ρ C > 1/N ) is given by Natural selection favors cooperation if b/c > (b/c) * .For example, if the underlying network is regular with degree k, we have and such that as N → ∞ [15].Note that the right-hand side of Eq. ( 8) only depends on the adjacency matrix of the network.In other words, the structure of the contact network determines whether and how much natural selection favors cooperation.

Perturbation theory for graph surgery
In this section, we develop a perturbation theory to determine the change in (b/c) * when one perturbs the weight of a single edge.To this end, we start by rewriting Eq. ( 6) in terms of matrices and vectors.Let 1 = (1, . . ., 1) , where represents the transposition.Let T = (t ij ) be the N × N matrix of the mean coalescence time.Using these notations, we rewrite Eq. ( 6) as where k = 1, 2, 3, and • represents the Hadamard product.
If one changes the weight of an edge (i 0 , j 0 ) by ε, where |ε| 1, including the case in which we create a new edge with weight ε (> 0), the perturbed network remains connected and undirected.Therefore, one can still use Eq. ( 8) to compute (b/c) * .Equation (8) uses Eq. ( 6), which requires π, P , and T .We denote these variables after the perturbation by π(ε), P (ε), and T (ε).To distinguish the quantities before and after the perturbation, we denote these variables before the perturbation by π(0), P (0), and T (0).
For writing down π(ε), we denote by the sum of the weighted degree of over all the nodes.Under a small perturbation, we carry out Taylor expansion of Eq. ( 1) to obtain where ∆π = (∆π 1 , . . ., ∆π N ).We obtain where δ ij is the Kronecker delta.We present the derivation of Eq. ( 16) in Appendix A.
We write ∆M as a block matrix where each ∆ ij is an N ×N matrix.We show the derivation of each ∆ ij in Appendix C. We point out that the number of nonzero rows of ∆M is equal to 2 To derive the first-order term of T (ε) from ∆M , we use Eq. ( 27) to obtain Therefore, we obtain where ∆T is the Finally, using Eq. ( 13), we derive the perturbed τ k (ε) as follows: where By substituting Eq. ( 32) in Eq. ( 8), we obtain where 4 Time complexity To calculate (b/c) * for a network with N nodes, the original algorithm requires calculating the mean coalescence time by solving a linear system of N (N − 1)/2 variables, i.e., t ij (with i, j ∈ {1, . . ., N } and i < j), which has a time complexity of O(N 6 ).
We now discuss the computational complexity of our perturbation method.Because the inner product of N -dimensional vectors has a time complexity of O(N ), the first while loop in Algorithm 1 has a complexity of O(N 2 ).The second while loop computes vec(∆T ).Because the scalar multiplication of an N 2 -dimensional vector requires O(N 2 ) time, the entire while loop has a time complexity of O(N 3 ).Therefore, for a single perturbation experiment, one can carry out the entire algorithm in O(N 3 ) time to obtain the perturbed {t ij }, and hence (b/c) * .This is considerably smaller than O(N 4.75 ) and O(N 6 ) with the Coppersmith-Winograd algorithm and the standard algorithm, respectively.The entire procedure to determine the single edge to be removed to maximize cooperation with the perturbation theory requires O(N 3 |E|) time in general networks and O(N 4 ) time for sparse networks.

Data
We use the following four synthetic networks and seven empirical networks in our numerical analysis in section 6.We show the number of nodes and that of edges for each network in Table 1 and visualize them in Fig. 2. All the networks are connected networks without self-loops.
We use a network generated by the Erdős-Rényi (ER) random graph with N = 100 nodes.We connect 300 pairs of nodes out of the N (N − 1)/2 = 4950 pairs of nodes selected uniformly at random.The average degree k = 6.
The planted -partition model, also called the random partition (RP) graph, partitions the set of N nodes into groups, each of which has N/ nodes [23].Any pair of nodes in the same group is adjacent to each other with probability p in .Any pair of nodes belonging to different groups are adjacent to each other with probability p out .If p in > p out , the intra-cluster edge density exceeds the inter-cluster edge density such that the network has community structure.We set N = 100, = 2, p in = 0.11, and p out = 0.01 such that the mean degree k = p in (N/ − 1) + p out N ( − 1)/ = 5.89 in theory.We use a network generated by this model having k = 6.12.
The Lancichinetti-Fortunato-Radicchi (LFR) model generates networks with community structure [24].The model generates a power-law degree distribution with power-law exponent γ, and a power-law distribution of the size of the community with power-law exponent κ.The model also requires the maximal degree k max and mean degree k as input.The mixing parameter µ ∈ (0, 1) specifies the fraction of edges that connect different communities.A small value of µ leads to strong community structure.We set N = 100, γ = 3, κ = 2, k = 6, k max = 100, and µ = 0.1.A network generated by this model that we use has k = 6.08.
We consider the following seven empirical networks.The karate club network consists of 34 nodes and 78 edges [25].Each node represents a member of a karate club in a university in the United States, who were observed between 1970 and 1972.The edges represent interaction outside the activities of the club.
The weaver network has 42 nodes and 151 edges [26].Each node represents a sociable weaver (Philetairus socius) observed in Benfontein Game Farm, Kimberley, South Africa.The observation lasted for 10 months in total: September-December 2010 and 2011, and January-February 2013.Two nodes are adjacent to each other if the two weavers used the same nest chambers either for roosting or nest-building within a series of observations in the same year.
The sparrow network has 52 nodes and 516 edges [27].A node represents a goldencrowned sparrow (Zonotrichia atricapilla) observed at the University of California, Santa Cruz Arboretum.The data was recorded between January and March 2010 [27].Although the original network is weighted, we regard this network as an unweighted network.
The lizard network has 60 nodes and 318 edges [28].Each node represents a lizard (Tiliqua rugosa) observed in a chenopod shrubland near Bundey Bore Station in South Australia.Each lizard was attached to the dorsal surface of the tail a data logger unit, which recorded synchronized GPS locations every 10 minutes.Two lizards were regarded to be adjacent to each other if they were within 2 meters of each other in any GPS record.
The dolphin network has 62 nodes and 159 edges [29].Each node represents a bottlenose dolphin (Tursiops).An edge represents a frequent association between two dolphins.
The email network has 167 nodes and 3251 edges [30].Each node represents an employee of a mid-sized manufacturing company in Poland.An edge between two nodes (i.e., employees) indicates that there exists at least one email correspondence between the two individuals.We do not distinguish the senders and the recipients and treat the network as undirected network.
The bird network has 202 nodes and 11900 edges [31].In the experiment, they placed some nest boxes in Wytham Woods, Oxford, UK, for six days to record individuals that landed on the entrance hole while prospecting for breeding territories.Each node represents a wild bird, which is either great tit (Parus major ), blue tit (Cyanistes caeruleus), marsh tit (Poecile palustris), coal tit (Periparus ater ), or Eurasian nuthatch (Sitta europaea).An edge represents two birds that overlapped in nest-box exploration patterns on the same day.
6 Numerical results

Addition or removal of a single edge
We examine the accuracy at which our perturbation theory describes the change in (b/c) * when we add or remove an edge in the given unweighted network.Before the perturbation, w ij = w ji = 1 if there exists an edge between the ith and jth nodes, and w ij = w ji = 0 otherwise.In the case of edge addition, we add an edge with weight ε between a pair of nodes (i 0 , j 0 ) without an edge in the original network unless we state otherwise, where 0 < ε ≤ 1.Therefore, w i0j0 (= w j0i0 ) changes from 0 to ε, and all the other w ij ∈ {0, 1} values remain unchanged.The addition of an unweighted edge corresponds to ε = 1.In the case of edge removal, we reduce the weight of an edge (i 0 , j 0 ) in the original network by −ε, where −1 ≤ ε < 0. Therefore, w i0j0 (= w j0i0 ) changes from 1 to 1 + ε, and all the other w ij values remain unchanged.The complete removal of an unweighted edge corresponds to ε = −1.
We are interested in whether the linear approximation to (b/c) * (ε) given by Eq. ( 34 We show the relationship between (b/c) * (ε) − (b/c) * (0) and ε when we reduce the weight of a single edge in a BA network with N = 100 nodes in Fig. 3(a).Each line in the figure corresponds to an edge whose weight is gradually reduced.Note that ε = 0 corresponds to the original network.Figure 3(a) indicates that (b/c) * roughly monotonically decreases as we gradually decrease the edge weight (i.e., decrease ε from 0 to negative values) except near ε = 0.For this network, the removal of any single edge (i.e., ε = −1) leads to a decrease in (b/c) * , implying that the edge removal promotes cooperation.However, we note that a small decrease in the weight of an edge in the original network (e.g., ε = −0.3)increases (b/c) * for some edges, making cooperation more difficult than in the original network.Figure 3(a) implies that the perturbation theory is not accurate at describing the amount of the change in (b/c) * upon the edge removal because most of the curves shown in the figures, corresponding to the different edges in the original network, are far from being linear.However, we observe that the curves with the largest values of the slope of the curve at ε = 0 tend to yield the smallest values of (b/c) * at ε = −1.Therefore, the perturbation theory, which produces the slope value, is expected to be efficient at detecting the edges whose removal yields the largest decrease in (b/c) * .
We show in Fig. 3(b) the change in (b/c) * plotted against ε when we add a new edge with weight ε.Each line corresponds to a pair of nodes between which there is initially no edge.Note that ε = 1 corresponds to the addition of an unweighted edge.We find that the addition of any unweighted edge increases (b/c) * , making cooperation difficult.However, in contrast to the case of edge removal, the addition of an unweighted edge (i.e., with edge weight ε = 1) does not necessarily yield the largest change in (b/c) * among edges of different weights ε ∈ (0, 1].Specifically, for many node pairs that are initially not adjacent to each other, adding an edge with an intermediate edge weight (e.g., ε ≈ 0.7) maximizes the increase in (b/c) * (see Fig. 3(b)).Another observation is that the slope of the curve at ε = 0, corresponding to the perturbation theory, is apparently less predictive of the effect of adding an unweighted edge (i.e., ε = 1).Specifically, Fig. 3(b) indicates that, even if the slope at ε = 0 is large, (b/c) * at ε = 1 can be relatively small because (b/c) * decreases as ε increases when ε is close to 1. Furthermore, the curves with the largest slopes at ε = 0 do not yield the largest changes in the (b/c) * value at ε = 1, which implies that the perturbation theory is expected to be inefficient at predicting the edge addition that makes the cooperation most difficult.
We find similar results for the planted 2-partition model for the gradual removal of a single edge (see Fig. 3(c)).A notable difference from the case of the BA model is that there exists one edge whose complete removal increases (b/c) * , making the cooperation difficult.The two nodes forming this edge have degrees 2 and 9, which are not outstanding.Furthermore, we have confirmed by running a deterministic approximate modularity maximization algorithm [32], using function greedy modularity communities in NetworkX, that these two nodes belong to the same community among the four communities detected.Therefore, this particular edge looks like just a normal edge.
We show in Fig. 3(d) the dependence of (b/c) * on ε when we gradually increase the weight of an edge that is initially absent in the planted 2-partition network.The slope of the curve at (b/c) * at ε = 0 is apparently not strongly related to the change in (b/c) * at ε = 1.
We show the results of edge removal in the dolphin network in Fig. 3(e).There are two edges out of the 150 edges of which the removal (i.e., ε = −1) increases (b/c) * , making cooperation difficult.These two edges are formed by two nodes with degrees 2 and 5 and two other ones with degrees 2 and 7.These degree values are not outstanding in the entire network.The four nodes belong to the same community among the four communities detected by the same approximate modularity maximization algorithm [32].These results suggest that the two edges apparently look normal.The removal of any other edge decreases (b/c) * , enhancing cooperation.Similar to the BA model, the curves with the largest slopes at ε = 0 yield the largest decreases in (b/c) * at ε = −1.
We show in Fig. 3(f) the dependence of (b/c) * on ε when we gradually increase the weight of an edge that is initially absent in the dolphin network.The results are similar to those for the planted 2-partition model shown in Fig. 3(d).Many curves yield decrease in (b/c) * at ε = 1, implying that the edge addition can promote cooperation, whereas the converse is the case for many other curves.The slope of the curve of (b/c) * at ε = 0 is apparently not strongly related to the change in (b/c) * at ε = 1.
The nonlinearity in the curves shown in Fig. 3 indicates that our perturbation theory is not accurate at predicting the amount of change in (b/c) * when we completely remove or add an edge in most cases.Therefore, we turn to ask whether the slope obtained from the perturbation theory is useful at determining the edge whose Table 1 Pearson correlation coefficient, r, between the shift in (b/c) * obtained by direct numerical simulations and that predicted by the perturbation theory.We remind that N is the number of nodes and that |E| is the number of edges.A large positive value of r upon edge addition or enhancement implies that the perturbation theory is good at predicting the outcome of adding or enhancing an edge.A large negative value of r upon edge removal implies that the perturbation theory is good at predicting the outcome of removing an edge., suggesting that the perturbation theory is not good at predicting the outcome of adding an edge, whereas the correlation coefficient is significant due to a large sample size (r = −0.36,n = 4659, p < 0.01).Note that a large positive correlation coefficient when we add an edge would imply that the perturbation theory is good at predicting the outcome of adding an edge.We show in Figs.4(c) and 4(d) the results for the same correlation analysis for the planted 2-partition model network.When one removes an existing edge, the change in (b/c) * and slope ∆(b/c) * are strongly negatively correlated (r = −0.80,n = 306, p < 0.01; see Fig. 4(c)), which is similar to the result for the BA model shown in Fig. 4(a), suggesting that the perturbation theory is good at predicting the outcome of removing an edge.When one adds a new edge, the change in (b/c) * and slope ∆(b/c) * are weakly correlated for this network (r = −0.39,n = 4644, p < 0.01; see Fig. 4(d)), which is similar to the result for the BA model shown in Fig. 4(b).

Network
We show the corresponding results for the dolphin network in Figs.4(e) and 4(f).The change in (b/c) * and slope ∆(b/c) * are strongly negatively correlated when one removes an edge (r = −0.72,n = 150, p < 0.01; see Fig. 4(e)) and less strongly correlated when one adds a new edge (r = 0.56, n = 1732, p < 0.01; see Fig. 4(f)).A strongly negative correlation for the edge removal (i.e., r = −0.72) is similar to the result for the BA model.A positive correlation for the edge addition (i.e., r = 0.56) implies that the perturbation theory is to some extent good at predicting the outcome of adding an edge.
We show in Table 1 the same relationships for the other networks.For all synthetic and empirical networks, the slope ∆(b/c) * obtained from perturbation theory is strongly negatively correlated with the change in (b/c) * when we remove an existing edge (r ≤ −0.72).Therefore, the perturbation theory is effective at predicting the outcome of removing an edge across different networks.However, the correlation is strongly positive only for a small fraction of networks (i.e., r ≥ 0.5 for three out of the nine networks) when we add a new edge to the network.

Enhancement of the weight of an existing edge
In this section, we allow weighted networks and consider an increase or decrease in the weight of an existing edge of the network.Because we effectively analyzed the case of the decrease in the edge weight in section 6.1 (i.e., by setting −1 < ε < 0), here we only consider enhancement of the weight of an existing edge by ε.
We enhanced the weight of an existing edge by 0 < ε ≤ 1, making the edge weight 1 + ε, and numerically examined (b/c) * in the altered weighted networks.We plot in Figs.5(a), 5(c), and 5(e) the change in (b/c) * relative to the original network against ε for the three networks used in Figs. 3 and 4. For the BA model, increasing the weight of 74 out of the 291 existing edges from 1 to 2 (i.e., ε = 1) led to an increase in (b/c) * , making cooperation more difficult, whereas the opposite is the case when one enhances the weight of any other edge (see Fig. 5(a)).This result contrasts with the case of adding a new edge to the same network, which always increases (b/c) * (see Fig. 3(b)).In the planted 2-partition model (Fig. 5(c)) and the dolphin network (Fig. 5(e)), enhancing the edge weight of 7 out of the 306 edges and 23 out of the 159 edges, respectively, led to an increase in (b/c) * .Therefore, in a majority of cases, cooperation becomes easier by enhancing the weight of a single edge, which contrasts with the results for adding a new edge to these networks (see Figs. 3(d) and 3(f)).These results altogether suggest that adding new edges and enhancing the weight of existing edges often lead to different results.We have confirmed that a high accuracy also holds true for other networks (see the last column of Table 1).These high accuracy results are in stark contrast to the results in case of adding a new edge, with which the accuracy of the perturbation theory is low.

Sequential edge removal
The nonlinearity in the curves shown in Fig. 3, and the results shown in Fig. 4 and Table 1 indicate that our perturbation theory is not accurate at estimating the amount of change in (b/c) * upon an edge removal.Therefore, we turn to investigate whether our perturbation theory is good at finding edges to be sequentially removed to decrease (b/c) * by a large amount in larger networks.Denote by G 0 an original network.We remove the edge with the largest ∆(b/c) * , resulting in network G 1 .Then, we calculate ∆(b/c) * for each existing edge in G 1 and remove the edge with the largest ∆(b/c) * , resulting in network G 2 .We repeat this procedure another three times to eventually obtain network G 5 , which has five fewer edges than G 0 .
A simple rule of thumb to determine edges to be removed to enhance cooperation is to use the degree of nodes composing the edge.In particular, (b/c) * for the death-birth rule is small for random regular graphs with small degrees [13] and general networks with a small mean degree [15].Therefore, we test the performance of our perturbation theory against a degree-based heuristic to remove an edge for enhancing cooperation, which we define as follows.Denote by (i, j) the edge to be removed and by k i and k j the degree of the ith and jth nodes, respectively.Note that k i = N =1 w i (= N =1 w i ) for our networks, which are unweighted.For each network, we remove the edge whose k i + k j is largest.After removing an edge according to this criterion, we select the edge with the largest k i +k j in the reduced network and remove it.We repeat this procedure another three times to remove five edges in total.In our numerical experiments described below, we have verified that the selected edges are always the same if the score for the edge is defined by k i k j instead of k i + k j .
We carry out sequential edge removal experiments on three synthetic networks and three empirical networks.Note that the six networks are mostly larger than those used in the previous numerical simulations.For these networks, it is computationally difficult to exactly calculate (b/c) * for all possible networks with, for example, one edge being removed from the original network.
We show the change in (b/c) * relative to the original network as we sequentially remove five edges using our perturbation theory by the red lines in Fig. 6.As expected, (b/c) * decreases, corresponding to negative ∆(b/c) * values, as we remove edges one by one.We also show the result of the sequential edge removal based on the degree sum k i + k j by the blue lines in the same figure.For all networks, there are multiple edges that have the same value of k i + k j at least in one of the five steps to remove a single edge.In this case, we calculated ∆(b/c) * for all the possible scenarios of removing one of the edges that maximize k i + k j in each step of edge removal.This is why we have obtained multiple blue lines in the figure.In all cases, (b/c) * decreases as we sequentially remove edges with the largest k i + k j value.Figure 6 indicates that the edge removal based on our perturbation theory results in a larger decrease in (b/c) * than that based on k i + k j for all the networks.To be quantitative, we measured the decrease in (b/c) * after the removal of five edges compared to the original network with the perturbation theory and with the degree sum.The former was larger than the average of the latter (i.e., average of the blue lines in Fig. 6) by a factor of 1.02, 1.01, 1.02, 1.05, 1.02, and 1.02 for the ER random graph (Fig. 6

Conclusions
To determine (b/c) * for an arbitrary network, one needs to solve a system of N 2 linear equations such that the time complexity is O(N 6 ).With the Coppersmith-Winograd algorithm, the time complexity is reduced to O(N 4.75 ), but this is still large (see section 4).In particular, it is computationally costly to carry out graph surgery with various possible edges to be added or removed to compare the results in terms of (b/c) * .Therefore, we have developed a perturbation theory for the graph surgery with which we can evaluate the perturbed (b/c) * in O(N 3 ) time.We have verified that the first-order term ∆(b/c) * obtained from our perturbation theory predicts the rank of the change in (b/c) * when one removes an edge from the network with a high accuracy.Specifically, we have numerically shown that the edge with the largest ∆(b/c) * value is the one whose actual removal decreases (b/c) * by the largest amount in two out of the three networks (see Fig. 4(a), 4(c), and 4(e)).Therefore, we conclude that our perturbation theory is useful for finding the edge whose removal efficiently enhances cooperation in the given network with a reduced computational cost.
We focused on the death-birth process because it tends to foster cooperation compared to other rules of strategy updating [11,13].However, it is straightforward to formulate similar perturbation methods in the case of other updating rules such as the birth-death process [13,18] and the pairwise comparison rule [16,[33][34][35] as well as in the case of other payoff matrices.In particular, our theory should be applicable to the case of constant selection [18,36], with which the payoff matrix is independent of the opponent's action.The perturbation theory may be more accurate for other update rules or games than the combination of the death-birth rule and the prisoner's dilemma game examined in the present study.Exploitation of our perturbation approach in these directions is left for future work.
Another direction of future work is interaction between the selection strength and network perturbation.In the present work, we have assumed the weak selection limit.However, one can retain a selection strength parameter (which is η in this article) to be finite and write down a formal solution.Then, it may be interesting to consider the simultaneous limit of weak selection η → 0 and weak network perturbation ε → 0 in a way η and ε are interrelated.Apart from this research direction, assessing the validity of the present perturbation theory under strong selection is left for future work.To this end, we first need to understand the accuracy of the original theory of fixation of cooperation in networks [15], which our theory is based on, under strong selection.
We do not know why the perturbation theory is more accurate when one removes an edge than when one adds an edge.Furthermore, we have found that the perturbation theory is fairly accurate at predicting the result for adding a parallel edge where an edge already exists, whereas it is not accurate when adding an edge where an edge does not exist in the original network.In a related vein, we observed nonmonotonic behavior in the cooperativity in terms of (b/c) * especially when we gradually added a weighted edge (Figs.3(b) and 3(f)).These results lead us to hypothesize that we can engineer networks that promote cooperation better by considering weighted networks than unweighted networks.These topics also warrant future work.

Fig. 2
Fig. 2 Visualization of the networks used in the numerical analysis.(a) ER random graph.(b) BA model.(c) Planted 2-partition model.(d) LFR model.(e) Karate club.(f) Weaver.(g) Sparrow.(h) Lizard.(i) Dolphin.(j) Email.(k) Bird.All the networks are undirected.The linear size of the node is proportional to its degree.We have ignored the weight of the edge in this figure and our analysis.
removal or addition changes (b/c) * by a large amount, representing strong promotion or suppression of cooperation in networks.We show in Fig.4(a) the relationship between the change in (b/c) * when we remove an edge from the BA network and the slope ∆(b/c) * obtained from Eq.(34).The two quantities are strongly negatively correlated (Pearson correlation coefficient r = −0.86,sample size n = 291, p < 0.01).This result indicates that the perturbation theory, which is theoretically accurate only in the vicinity of ε = 0, is good at predicting the outcome of removing an edge.We show in Fig.4(b) the change in (b/c) * when we add a new edge to the same BA network as a function of the slope, ∆(b/c) * .The change in (b/c) * is not strongly positively correlated with ∆(b/c) *

Fig. 3
Fig. 3 Change in (b/c) * as a function of the change in the edge weight, ε.(a) BA model, removal of an existing edge.(b) BA model, addition of a new edge.(c) Planted 2-partition model, removal of an existing edge.(d) Planted 2-partition model, addition of a new edge.(e) Dolphin network, removal of an existing edge.(f) Dolphin network, addition of a new edge.In (a), (c), and (e), each line represents an edge in the original network.In (b), (d), and (f), each line represents a pair of nodes that is not adjacent to each other in the original network.The line color is only as a guide to the eyes.

Fig. 4
Fig. 4 Change in (b/c) * when we remove or add an unweighted edge as a function of the slope ∆(b/c) * of the curves shown in Fig. 3 at ε = 0. (a) BA model, removal of an existing edge.(b) BA model, addition of a new edge.(c) Planted 2-partition model, removal of an existing edge.(d) Planted 2-partition model, addition of a new edge.(e) Dolphin network, removal of an existing edge.(f) Dolphin network, addition of a new edge.Each circle in (a), (c), and (e) represents an edge in the original network.Each circle in (b), (d), and (f) represents a pair of nodes that is not adjacent to each other in the original network.

Figures 5 (
a), 5(c), and 5(e) also indicate that the change in (b/c) * is close to linear as a function of ε.Therefore, our perturbation theory should be accurate at estimating the change in (b/c) * with ε = 1.To verify this prediction, we show in Figs.5(b), 5(d), and 5(f) the relationship between the change in (b/c) * in response to changing the weight of a single edge from 1 to 2 and ∆(b/c) * obtained by the perturbation theory for the three networks.As expected, the accuracy of the perturbation theory is high.

Fig. 5
Fig. 5 Change in (b/c) * when one enhances the weight of an existing edge.Panels (a), (c), (e): Change in (b/c) * as a function of the increase in the weight of an existing edge, ε.Each line represents an edge in the original network.Panels (b), (d), (f): Change in (b/c) * when we enhance an edge weight by ε = 1, plotted against the slope ∆(b/c) * of the curves shown in panels (a), (c), and (e) at ε = 0.Each circle represents an edge in the original network.(a) and (b): BA model.(c) and (d): Planted 2-partition model.(e) and (f): Dolphin network.

Fig. 6
Fig. 6 Changes in (b/c) * upon sequential removal of five edges.(a) ER random graph with 300 nodes and 900 edges.(b) BA model network with 300 nodes and 891 edges.(c) Planted 2-partition network with 300 nodes and 939 edges.(d) Lizard network with 60 nodes and 318 edges.(e) Email network with 167 nodes and 3251 edges.(f) Bird network with 202 nodes and 11900 edges.The red lines represent the edge removal according to the perturbation theory.The blue lines represent the edge removal according to the rank of the degree sum.
[15])[15]. To determine the single edge whose removal decreases (b/c) * by the largest amount, for example, one needs to repeat this procedure for each edge.Therefore, the entire procedure with an ordinary algorithm and theCoppersmith-Winograd algorithm requires O(N 6 |E|) and O(N 4.75 |E|) time, respectively, where |E| is the number of edges.For a sparse network, for which |E| = O(N ), the time complexity is O(N 7 ) and O(N 5.75 ), respectively.