Perturbation theory for evolution of cooperation on networks

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 [12, 13, 14, 11, 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 [16, 17, 18, 5]. 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 $i$th entry encodes the type of the $i$th 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 [18, 13]. 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.

We assume that the graph $G$ is connected and undirected. We denote the set of nodes by $V=\{1,\ldots,N\}$, where $N$ is the number of nodes. For each pair of nodes $i,j\in V$, we denote the edge weight by $w_{ij}\geq 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}=\sum_{j=1}^{N}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\times 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=\mathrm{diag}(s_{1},\ldots,s_{N})$, i.e., the diagonal matrix whose diagonal entries are equal to $s_{1},s_{2},\ldots,s_{N}$. Let $\boldsymbol{\pi}=(\pi_{1},\ldots,\pi_{N})$ be the stationary probability vector of the random walk with transition probability matrix $P$, i.e., the solution of $\boldsymbol{\pi}P=\boldsymbol{\pi}$. 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 $\boldsymbol{x}=(x_{1},\ldots,x_{N})\in\{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 $\boldsymbol{x}$ is given by

where $\eta$ represents the strength of the selection. If $\eta=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 $\eta\to 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}(\boldsymbol{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 $\rho_{\mathrm{C}}$ the expectation that the cooperation fixates under this distribution of the initial state. If $\rho_{\mathrm{C}}>1/N$, natural selection favors cooperation [16, 13, 5, 15]. In Ref. [15], Allen et al. showed that

where

$p_{ij}^{(k)}$ is the $(i,j)$th entry of matrix $P^{k}$, which implies that $p_{ij}^{(1)}=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_{ij}^{(k)}$ is the $k$-step transition probability of the random walk from node $i$ to node $j$. Therefore, $\tau_{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., $\rho_{\mathrm{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\to\infty$ [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.

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 $\boldsymbol{1}=(1,\ldots,1)^{\top}$, where ^⊤ represents the transposition. Let $T=(t_{ij})$ be the $N\times N$ matrix of the mean coalescence time. Using these notations, we rewrite Eq. (6) as

where $k=1,2,3$, and $\circ$ represents the Hadamard product.

If one changes the weight of an edge $(i_{0},j_{0})$ by $\varepsilon$, where $|\varepsilon|\ll 1$, including the case in which we create a new edge with weight $\varepsilon$ $(>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 $\boldsymbol{\pi}$, $P$, and $T$. We denote these variables after the perturbation by $\boldsymbol{\pi}(\varepsilon)$, $P(\varepsilon)$, and $T(\varepsilon)$. To distinguish the quantities before and after the perturbation, we denote these variables before the perturbation by $\boldsymbol{\pi}(0)$, $P(0)$, and $T(0)$.

For writing down $\boldsymbol{\pi}(\varepsilon)$, 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 $\Delta\boldsymbol{\pi}=(\Delta\pi_{1},\ldots,\Delta\pi_{N})$. We obtain

where $\delta_{ij}$ is the Kronecker delta. We present the derivation of Eq. (16) in Appendix A.

To calculate $P(\varepsilon)$, we define a symmetric indicator function, denoted by $\chi_{i_{0}j_{0}}$, by

We obtain

where $\Theta^{(1)}=(\theta^{(1)}_{ij})$, $\Theta^{(2)}=(\theta^{(2)}_{ij})$, and $\Theta^{(3)}=(\theta^{(3)}_{ij})$ are $N\times N$ matrices whose entries are given by

and

We show the derivation of Eqs. (21), (22), and (23) in Appendix B.

We next calculate $T(\varepsilon)$. Matrix $T(0)=(t_{ij}(0))$ satisfies

which we obtain by applying $t_{ij}(0)=t_{ji}(0)$ to Eq. (7). Note that $\{p_{11}(0),p_{12}(0),\ldots,p_{NN}(0)\}$ are known from the network structure and that $\{t_{11}(0),t_{12}(0),\ldots,t_{NN}(0)\}$ are unknowns. We stack Eq. (24) for the different $i$ and $j$ values in lexicographical order of $(i,j)$ on the left-hand side. In other words, the first equation is $t_{11}(0)=0$, the second equation is $t_{12}(0)-\frac{1}{2}p_{11}(0)t_{12}(0)-\frac{1}{2}\sum_{k=3}^{N}p_{1k}(0)t_{2k}(0)-\frac{1}{2}\sum_{k=2}^{N}p_{2k}(0)t_{1k}(0)=1$, the third equation is $t_{13}(0)-\frac{1}{2}p_{11}(0)t_{13}(0)-\frac{1}{2}p_{12}(0)t_{23}(0)-\frac{1}{2}\sum_{k=4}^{N}p_{1k}(0)t_{3k}(0)-\frac{1}{2}\sum_{k=2}^{N}p_{3k}(0)t_{1k}(0)=1$, and so on. Denote by $\mathrm{vec}(T(0))$ the thus obtained vectorization of matrix $T(0)$, i.e.,

Equation (25) is a redundant expression because $T(0)$ is a symmetric matrix and its diagonal elements are equal to 0. However, we use Eq. (25) in the following text because it makes the theoretical derivations and computational implementation easier than the most compact vector form of $T(0)$, which would be $N(N-1)/2$-dimensional. Using Eq. (25), we rewrite Eq. (24) as

where $M(0)$ is the $N^{2}\times N^{2}$ matrix whose entries are determined by Eq. (24), and $\boldsymbol{d}$ is the $N^{2}$-dimensional column vector whose $((k-1)N+k)$th entry is equal to 0 for all $k\in\{1,\ldots N\}$, and all the other entries are equal to 1. Because it also holds true that $M(\varepsilon)\mathrm{vec}(T(\varepsilon))=\boldsymbol{d}$, the calculation of $T(\varepsilon)$ requires $M(\varepsilon)$, which is the matrix with perturbation, defined similarly to $M(0)$. We obtain the entries of $M(\varepsilon)$ by those of $M(0)$ with each $p_{ij}(0)$ (with $i,j\in\{1,\ldots,N\}$) being replaced by $p_{ij}(\varepsilon)$. We write the Taylor expansion of $M(\varepsilon)$ as

and calculate $\Delta M$ as follows.

We write $\Delta M$ as a block matrix

where each $\Delta_{ij}$ is an $N\times N$ matrix. We show the derivation of each $\Delta_{ij}$ in Appendix C. We point out that the number of nonzero rows of $\Delta M$ is equal to $2\times(N-2)+(N-1)\times 2=4N-6$, which is much smaller than $N^{2}$ for a large $N$ (see Appendix C).

To derive the first-order term of $T(\varepsilon)$ from $\Delta M$, we use Eq. (27) to obtain

Therefore, we obtain

where $\Delta T$ is the $N\times N$ matrix satisfying

Finally, using Eq. (13), we derive the perturbed $\tau_{k}(\varepsilon)$ as follows:

where

By substituting Eq. (32) in Eq. (8), we obtain

where

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\in\{1,\ldots,N\}$ and $i<j)$, which has a time complexity of $O(N^{6})$. With the Coppersmith-Winograd algorithm [21], the time complexity is reduced to $O(N^{4.75})$ [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 the Coppersmith-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.

The matrix $\Delta M$ defined by Eq. (28) is sparse and has a special pattern. If the $i$th row of $\Delta M$ is a zero row, then the $i$th element of vector $\Delta M\mathrm{vec}(T(0))$ is zero, and we do not need to calculate it. Therefore, to calculate $\Delta M\mathrm{vec}(T(0))$, we only need to focus on its $((i_{0}-1)N+k)$th entries, where $k\in\{1,\ldots,N\}\setminus\{i_{0}\}$, $((j_{0}-1)N+k)$th entries, where $k\in\{1,\ldots,N\}\setminus\{j_{0}\}$, and $((k-1)N+i_{0})$th and $((k-1)N+j_{0})$th entries, where $k\in\{1,\ldots,N\}\setminus\{i_{0},j_{0}\}$. All the other entries of $\Delta M\mathrm{vec}(T(0))$ are equal to $0$. We show a pseudo algorithm to calculate $\Delta T$ in Algorithm 1.

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 $\mathrm{vec}(\Delta 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.

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 $\langle k\rangle=6$.

With the Barabási–Albert (BA) model, we sequentially add new nodes each with $m=3$ edges that connect to existing nodes according to the linear preferential attachment rule [22]. We start the growth process from the star graph with four nodes. The degree distribution approximately obeys $p(k)\propto k^{-3}$, where $p(k)$ is the probability that a node has degree of $k$, and $\propto$ represents “proportional to”, in the limit of $N\to\infty$. We set $N=100$ and $m=3$, which yields $291$ edges, implying $\langle k\rangle=5.82\approx 6$.

The planted $\ell$-partition model, also called the random partition (RP) graph, partitions the set of $N$ nodes into $\ell$ groups, each of which has $N/\ell$ nodes [23]. Any pair of nodes in the same group is adjacent to each other with probability $p_{\mathrm{in}}$. Any pair of nodes belonging to different groups are adjacent to each other with probability $p_{\mathrm{out}}$. If $p_{\mathrm{in}}>p_{\mathrm{out}}$, the intra-cluster edge density exceeds the inter-cluster edge density such that the network has community structure. We set $N=100$, $\ell=2$, $p_{\mathrm{in}}=0.11$, and $p_{\mathrm{out}}=0.01$ such that the mean degree $\langle k\rangle=p_{\mathrm{in}}(N/\ell-1)+p_{\mathrm{out}}N(\ell-1)/\ell=5.89$ in theory. We use a network generated by this model having $\langle k\rangle=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 $\gamma$, and a power-law distribution of the size of the community with power-law exponent $\kappa$. The model also requires the maximal degree $k_{\mathrm{max}}$ and mean degree $\langle k\rangle$ as input. The mixing parameter $\overline{\mu}\in(0,1)$ specifies the fraction of edges that connect different communities. A small value of $\overline{\mu}$ leads to strong community structure. We set $N=100$, $\gamma=3$, $\kappa=2$, $\langle k\rangle=6$, $k_{\mathrm{max}}=100$, and $\overline{\mu}=0.1$. A network generated by this model that we use has $\langle k\rangle=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 golden-crowned 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.