Reciprocity in spatial evolutionary public goods game on double-layered network

Spatial evolutionary games have mainly been studied on a single, isolated network. However, in real world systems, many interaction topologies are not isolated but many different types of networks are inter-connected to each other. In this study, we investigate the spatial evolutionary public goods game (SEPGG) on double-layered random networks (DRN). Based on the mean-field type arguments and numerical simulations, we find that SEPGG on DRN shows very rich interesting phenomena, especially, depending on the size of each layer, intra-connectivity, and inter-connected couplings, the network reciprocity of SEPGG on DRN can be drastically enhanced through the inter-connected coupling. Furthermore, SEPGG on DRN can provide a more general framework which includes the evolutionary dynamics on multiplex networks and inter-connected networks at the same time.


Previous Study
To understand "Tragedy of the commons" 48 problem with large participants has been studied through the SEPGG on the complete graph (CG) and dense random networks 17 . Depending on the multiplication factor r and the size of graph N, either Loner-only state (L-state), which is the anomalous state with no active participants, or Defector-only state (D-state), which means the state of "tragedy of the commons", has been shown to appear on CG 17 . Furthermore, we have shown the following crossover behaviors as the mean-degree 〈 k〉 of underlying random networks changes 17 . For small r, the L-state crosses over to the D-state and the D-state successively crosses over to the Cooperator-only state (C-state) as 〈 k〉 decreases. For large r, the direct crossover from the D-state to the C-state occurs as 〈 k〉 decreases. We have been found that cooperation gradually increases as the number of participants or 〈 k〉 decreases, which is the origin of these crossovers. Hence, the crossovers describe how the enhanced cooperation on sparse networks overcomes "tragedy of the commons" on dense networks.

Results
SEPGG on double-layered random network. Now we want to explain how the double-layered random networks (DRNs) are composed. The first random network (layer) with the size N 1 and the average intradegree 〈 k intra 〉 1 and the second random layer with N 2 (≥ N 1 ) and 〈 k intra 〉 2 are separately constructed. Then, to make DRN with n inter (≤ N 1 ) interlinks, n inter different nodes both on the first layer and the second layer are chosen randomly. Each chosen node on the layer 1 is made to be randomly linked to a chosen node of the second layer without making multiple interlinks to a certain node. (See Fig. 1 ). We call a DRN with N 1 = N 2 (= N) and 〈 k intra 〉 1 = 〈 k intra 〉 2 a symmetric DRN and a DRN with N 1 ≤ N 2 or 〈 k intra 〉 1 ≠ 〈 k intra 〉 2 an asymmetric DRN.
SEPGG model on a constructed DRN is defined as what follows. Each agent is assigned to a node on DRN. The strategy s i of the agent on a node i can be Cooperator (C), Defector (D) or Loner (L). In each update an agent i is randomly chosen. First, we calculate the payoff P i of i using the following rule. Let n i,C be the number of agents with C, n i,D be that with D and n i,L be that with L among the k i,intra + k i,inter + 1 agents. Here k i,intera (k i,inter ) is the intradegree (interdegree) of node i and n i,C + n i,D + n i,L = k i,intra + k i,inter + 1. P i is given by, Here, c is the cost contributed to the common pool by a C, r(> 1) is the multiplication factor and σ is the fixed payoff of an L. We imposed the condition 0 < σ < c(r − 1) 35 . Then, i changes its strategy through a biased imitation process as what follows 24 . With the probability p, the interlinked neighbor j on the opposite layer is selected. With the probability 1 − p, an intralinked neighbor j on the same layer is randomly selected. If i has no interlink, we choose a neighbor j from the same layer regardless of p. The strategy of i is changed into the strategy of j with the transition probability f ij , where Here, β(≥ 0) controls the amount of noise. When β → 0 i randomly adopts the strategy of j. However, for β > 0 we have shown that there exist distinctive states on random network 17 and 〈 k intra 〉 17 . Therefore, in this paper, we mainly use β = c = σ = 1 for the numerical analyses without loss of generality.
Results on the DRN with N 1 = N 2 = N and k inter = 1. Let's first study the SEPGG model on the DRN with N 1 = N 2 = N on which any agent i has one interlink or k i,inter = 1. We focus the steady-state densities . Here, ρ α,C (t) means the density of C on the layer α(= 1, 2) at the time t and etc. All quantities are averaged over 2000 realization of networks. For each network realization strategies are randomly assigned to each agent with given ρ ρ ρ , . When p = 0, there exists no coupling between the layers and the steady-states of each layers are the same as those on a single random network, which we have already studied in ref. 17. The followings are brief summary of previous results. Depending on r and 〈 k〉 α (= 〈 k intra 〉 α + 〈 k inter 〉 α ), the steady-state on each layer for p = 0 becomes one of the following 5 states. For r 0 = 0.3(< 1) For β = c = σ = 1 and σ Next we consider the case of p = 1(complete coupling). In complete coupling, any pair of interlinked nodes with identical strategies, i.e., a C-C pair, a D-D pair or an L-L pair, cannot be changed. Any pair with different strategies should change into a pair with identical strategies by the first successful transition. Therefore, the final steady-state is the absorbing state in which any interlinked pair has a common strategy and ρ ρ ρ ρ ρ ρ = = = , , 2, , there cannot be anomalous effects to break the density symmetry between layers during evolution. So, we expect ρ 1,C (t) = ρ 2,C (t), ρ 1,D (t) = ρ 2,D (t), ρ 1,L (t) = ρ 2,L (t) for p = 1, which is confirmed by simulations on various DRNs. To understand the steady-state behavior for p = 1 in a mean-field level, we study the PGG model on the symmetric DRN of two complete graphs (CGs). The payoff of a node on CG simply depends on ρ C , ρ D and ρ L 17 . When  N 1, the payoff of a node with D on both CG 1 and 2 is written as = Therefore, the transition probability f (1,D)(2,C) (t) that a D-node on the CG 1 accepts the C-strategy of the interlinked node on the CG 2 is equal to f (2,D)(1,C) and The first random layer with the size N 1 and the average intradegree 〈 k intra 〉 1 and the second random layer with N 2 (≥ N 1 ) and 〈 k intra 〉 2 are separately constructed.
To construct random layers we use the Erdös-Rényi (ER) network model 57 whose degree distribution is known to satisfy the Poisson distribution. Then, for DRN with n inter (≤ N 1 ) interlinks, n inter (≤ N 1 ) different nodes are randomly chosen on both layers. Then, n inter links are made, so that one-to-one correspondence between the chosen nodes on the first layer and those on the second layer occurs.
[ρ 1, Since the first successful update changes an interlinked pair of nodes with different strategies (active pair) into the pair with the same strategies (dead pair), the final absorbing state appears very rapidly. If the active pair is changed into the dead pair in accordance with the initial transition probabilities, then the steady-state for p = 1 on the DRN of two CGs is calculated by  As shown in Fig. 2(a), the mean-field equation (3) explains the simulation results on the DRN of two CGs for any r very well. When  r 1, , we get ρ = .  , gets smaller as r decreases for r < 20. Next, we study how the steady-state on the DRN with general 〈 k intra 〉 α 's for p = 1 varies from the mean-field result. In Fig. 2(b,c), the simulation results on symmetric DRNs are shown. The result for large r (or r = 30) in Fig. 2(b) deviates from the mean-field result. The deviation becomes much more enhanced as 〈 k intra 〉 α gets smaller.
also gets larger (smaller) compared to the mean-field result. ρ α L S , is nearly the same as the mean-field expectation for the very small 〈 k intra 〉 α . The result for small r (or r = 2) in Fig. 2(c) is nearly the same as the mean-field result for small 〈 k intra 〉 α . For small 〈 k intra 〉 α , ρ α C S , is slightly larger than the mean-field result and both ρ α D S , and ρ α L S , are slightly smaller than the mean-field results. In Fig. 3, the simulation results for p = 1 on asymmetric DRNs are shown. When p = 1 the interlinked pairs of the same strategies make dead pairs as in the symmetric DRN, the ρ α S 's on both layers are the same when N 1 = N 2 . As shown in Fig. 3, for a given Thus, the steady-state ρ α S 's on the asymmetric DRNs depends mainly on 〈 k intra 〉 1 and independent of 〈 k intra 〉 2 except for the symmetric range 〈 k intra 〉 2 ∼ 〈 k intra 〉 1 . For large r (r = 30) both ρ α C S , and ρ α L S , on the asymmetric DRN with the small 〈 k intra 〉 1 (Fig. 3(a)) is larger than that on the DRN with the large 〈 k intra 〉 1 (Fig. 3(b)), while ρ α D S , for the small 〈 k intra 〉 1 is smaller. For small r (r = 2) as shown in Fig. 3(c,d), the dependences of ρ α C S , and ρ α D S , on the asymmetric DRN on 〈 k intra 〉 1 are the same as those for large r. In contrast, ρ α L S , on the asymmetric DRN is nearly constant of the mean-field values regardless of 〈 k intra 〉 1 . When r gets larger in the interval r ≳ 20, ρ α C S , (ρ α D S , ) on the asymmetric DRN becomes slightly larger (smaller) than that for mean-field expectation, but ρ α L S , hardly varies. Now we want to explain the theoretical origins of the results for p = 1 in Figs 2 and 3. As explained when deriving the mean-field result, we use the fact that the final absorbing state for p = 1 appears very rapidly. Thus, the initial configuration is very important to decide the final absorbing state. Topologically localized cluster of a certain strategy on CG is impossible to form if there exist Cs, Ds and Ls simultaneously. In contrast, on the layer with relatively small 〈 k intra 〉 α , it is quite easy to form localized clusters of the same strategies due to fluctuation of the distribution of Cs, Ds and Ls. Such localized clusters reinforces the network reciprocity [16][17][18][19][20][21][22][23]25 . As 〈 k intra 〉 α decreases, it becomes much easier to form localized clusters. If a node i of a C-cluster on one layer is interlinked to a node j of a D-cluster on the opposite layer, then the node j is changed into a C-node when P i > P j or r > (n i,C + 1)(n j,D + 1)/[n i,C (n j,D + 1) − (n i,C + 1)] for p = 1. Since in RN n i,C and n j,D is roughly comparable with 〈 k intra 〉 α , small 〈 k intra 〉 α enhances the fluctuation. Thus, large r and small 〈 k intra 〉 α make more D-nodes in D-cluster be changed into C-nodes compared to the mean-field expectation. This enhances the network reciprocity as shown in Fig. 2(b,c). If a node of a C-cluster on one layer is interlinked to a node of an L-cluster on the opposite layer, the C-node wins over the L-node for large r, but the L-node wins over the C-node for small r. Therefore, the fluctuation effect makes more C and suppresses D compared to the mean-field result and this effect becomes more enhanced as 〈 k intra 〉 α decreases. In addition, it enhances L for small r and suppresses L for large r. Thus, on the symmetric DRN, ρ α C S , increases and ρ α D S , decreases as 〈 k intra 〉 α for large r as in Fig. 2(b). For small r, D-node is relatively hard to change its strategy into C. As a result the fluctuation effect is suppressed a little bit compared with large r and appears only on the DRN with small 〈 k intra 〉 α as in Fig. 2(c). On the asymmetric DRN, the C-cluster is more easily formed on the layer with smaller 〈 k intra 〉 α (or the layer 1), which makes the D-nodes in the opposite layer change their strategy into C. Thus the steady-state is decided by the layer 1. This effect becomes more enhanced as 〈 k intra 〉 1 decreases. This fluctuation effect on the asymmetric DRN explains the behaviors of ρ α C S , and ρ α D S , in Fig. 3. The fluctuation effect on ρ α L S , is relatively weaker as in Fig. 3. We now explain the results for the steady-state on the asymmetric DRN for 0 < p < 1, which show novel network reciprocity. On the asymmetric DRN with 〈 k intra 〉 1 and 〈 k intra 〉 2 (> 〈 k intra 〉 1 ), the steady-state should mainly depend on 〈 k intra 〉 1 as expected from the result for p = 1. To see this effect on the asymmetric DRN we first set 〈 k intra 〉 1 = 10, which is small enough that the steady-state on the layer 1 is the C-state when p = 0 (i.e., the state I and IV). If 〈 k intra 〉 2 is also small enough, the state is the trivial C-state on both layers. If r 0 < 1 and 〈 k intra 〉 2 is in the moderately large, then the D-state (the state II) appears on the layer 2 for p = 0.
We first carry out the simulations on the asymmetric DRN with 〈 k intra 〉 2 = 40 which is moderately large intradegree for r 0 = 0.3(< 1). As shown in Fig. 4(a) . The C-dominant state on the layer 1 exists for any p(< 1) as expected from the C-state(state I) for p = 0. The steady-state on the layer 2 for p > 0 is rather surprising, because ρ C S 2, becomes very large for even very small p. This indicates that the inter-layer coupling drastically enhances the network reciprocity. To see the origin of this result, we study ρ α,C (t), ρ α,D (t), ρ α,L (t). For small p(= 0.1), ρ 1,C (t) (ρ 1,D (t)) increases (decreases) rapidly at early time t as shown in Fig. 4(b). Due to the weak inter-layer coupling, when 20 ≲ t as the case p = 0. Since the layer 2 would be in D-state when p = 0 as shown in Fig. 4(a), ρ 2,D (t) increases and is slightly larger than ρ 2,C (t) at small t. Through the inter-layer coupling, D's on layer 2 can easily change their strategy into C due to the larger payoff of C on the layer 1 when ρ 1,C (t) becomes large enough. As a result ρ 2,C (t) increases to ρ .
when 20 ≲ t ≲ 30 both ρ 2,C (t) and ρ 2,D (t) follow ρ 1,C (t) and ρ 1,D (t). Thus, for small p, the C-dominance on the layer 1, which is fully developed due to weak inter-layer coupling, induces the C-dominance on the layer 2 at later time. For moderate p(= 0.6), the dependences of ρ 1,C (t) and ρ 1,D (t) on t are much more similar to those of ρ 2,C (t) and ρ 2,D (t) as shown in Fig. 4(c). For moderate p, the effects from the inter-layer coupling are nearly equal to the effect of the intra-layer interactions. Due to the increased inter-layer coupling, ρ 2,C (t) increases nearly synchronously with ρ 1,C (t). Such rapid increase of ρ 2,C (t) makes remnant D's on the layer 2. Those remnant D's on the layer 2 get relatively high payoff from the intra-layer interactions with relatively dense C's induced by inter-layer coupling. Then, through the inter-layer coupling, ρ C S 1, slightly decreases by remnant D's on the layer 2 for 0.2 ≲ p ≲ 0.6. For large p(> 0.6), ρ C S 1, increases as p increases. In this case, ρ α (t)'s first reach ρ α S for p = 1 very rapidly as shown in Fig. 4(d). Then, through the intra-layer interactions,  Fig. 4(a-d). ρ α L S , is nearly equal to 0, which can be understood from the result p = 0. These mechanisms explain the dependence of steady-state on p in Fig. 4(a) rather well. Next, we study the inter-layer coupling effects on the asymmetric DRN with 〈 k intra 〉 1 = 10 and 〈 k intra 〉 2 = 500 for r 0 = 0.3(< 1). For p = 0, the steady-state on the layer 1 is the C-state (the state I) and the steady-state on the layer 2 is the L-state (the state III) as shown in Fig. 5(a). As p increase, ρ = 1 for 0 < p < 1. When p ≲ 0.34, ρ α (t)'s show nearly the same behavior as those in Fig. 4(b). For larger 〈 k intra 〉 2 it is more difficult to make C's through intra-layer interaction and more inter-layer coupling or large p is needed to increase ρ C S 2, . The time-dependences as in Fig. 4(b) do not occur for the larger p(> 0.34) on the DRN with 〈 k intra 〉 1 = 10 and 〈 k intra 〉 2 = 500. As p increases further (p > 0.34), ρ α (t)'s behave nearly the same as in Fig. 4(c), which decreases (increases) ρ C S 's (ρ D S 's) on both layers. Due to large 〈 k intra 〉 2 (=500), ρ α (t)'s as in Fig. 4(d) do not happen. Instead, ρ α S 's for larger p approach to S 's smoothly. These mechanisms explain the dependence of steady-state on p in Fig. 5(a) rather well. We also study the SEPGG for r 0 = 10.0(> 1) on the asymmetric DRN with 〈 k intra 〉 1 = 100 and 〈 k intra 〉 2 = 2000. As shown in Fig. 5(b), the steady-state for p = 0 on the layer 1 is the C-state (state IV) and that on the layer 2 is the D-state (state V). For p > 0, ρ except for very small p. Thus, on this DRN, ρ α (t)'s show nearly the same behavior as those in Fig. 4(b) for p ≠ 0 due to the very large r. Furthermore, the payoff of D on layer 2 is less than that of C's in layer 1 in general. Thus, the D's in the layer 2 easily change into C's through the inter-layer coupling as p increases, which makes ρ p ( ) approaches to zero for p > 0.1. In Fig. 6 we display the behavior of ρ α p ( ) S 's for r 0 = 0.3(< 1) on the asymmetric DRN with 〈 k intra 〉 1 = 40 and 〈 k intra 〉 2 = 500. As shown in Fig. 6, the steady-state for p = 0 on the layer 1 is the D-state (state II) and that on the layer 2 is the L-state (state III). Except for . . In contrast, the delicate anomalous behavior around .  p 0 5 is confirmed to be originated from the balance between the inter-layer processes and the intra-layer processes, because for .  p 0 5, the intra-layer processes happen at nearly the same rate as the inter-layer processes. We confirm the following cyclic process I)-II)-III)-IV)-I) by simulations for .  p 0 5. I) The intra-layer interactions on the layer 1 make C-dense regions. II) This C-dense regions induce C's on the layer 2 through the inter-layer coupling. III) The C's on the layer 2 make D's by the intra-layer interactions. IV) The D's on the layer 2 shrink the C-dense regions through the inter-layer coupling. This cyclic process makes non-zero ρ α C S , and ρ α D S , around .  p 0 5. This anomalous effect also makes to appear the network reciprocity that ρ > α 0 C S , from the inter-layer coupling, without which there cannot exist any C on both layers.
We also study the model for 0 < p < 1 on the symmetric DRN. We check for various p and r and find that the steady-state densities on the symmetric DRN show exactly the same behavior as those on the single random network. We found that the steady-state on one layer inevitably is the same as that on the other layer. From the comparison of ρ α S 's for various p and r to confirm that the steady-state on the symmetric DRN for 0 < p < 1 is exactly the same as on the corresponding single network, and find that p only makes the time-delay (see Supplementary  Information).
Results on the DRN with 〈k intra 〉 < 1 or N 1 ≠ N 2 . We also study the model on the DRN with N 1 = N 2 = N and 〈 k intra 〉 < 1. In Fig. 7(a), the results on the asymmetric DRN with 〈 k intra 〉 1 = 10, 〈 k intra 〉 2 = 40 and 〈 k intra 〉 = 0.5 for r 0 = 0.3(< 1) are shown. Comparing Fig. 7(a) to Fig. 4(a), the dependences of ρ α S 's on p(< 1) for 〈 k intra 〉 = 0.5 are nearly the same as those up to p = 0.5 for k inter = 1. In contrast ρ α S 's at p = 1 for 〈 k intra 〉 = 0.5 show the nontrivial behaviors, because N〈 k intra 〉 = 0.5N interlinked pairs of nodes become dead pairs in the steady-state. On the DRN with 〈 k intra 〉 < 1 the dependences of ρ α S 's on p(< 1) are generally confirmed to be nearly the same as those up to p = 〈 k intra 〉 except at p = 1.
To know the effects of the difference between two layer sizes, we study the model on the DRN with N 1 ≠ N 2 . In Fig. 7(b), the results on the DRN with N 1 = 16000, 〈 k intra 〉 1 = 10 and N 2 = 32000, 〈 k intra 〉 2 = 57 for r 0 = 0.3(< 1) are shown. At p = 0, the steady-state of the layer 1 is the C-state (State I) and the steady-state of the layer 2 is the D-state (State II). As p increase, ρ = 1 . Because N 1 < N 2 , the intra-layer interactions on the layer 2 are stronger than those for N 1 = N 2 . The difference between the results in Fig. 7(b) and those in Fig. 5(a) are originated from the enhanced intra-layer interactions on the layer 2. Due to the enhanced intra-layer interactions in layer 2, the there are more D's than to the case of N 1 = N 2 and ρ D S 's (ρ C S 's) on both layer increase (decrease) compared to those in Fig. 5(a). We also study the model on the DRN with N 1 = 16000, 〈 k intra 〉 1 = 40 and N 2 = 32000, 〈 k intra 〉 2 = 14 and find nearly identical results to those in Fig. 5(b). Because N 1 < N 2 , the enhanced intra-layer interactions on the layer 2 induce the strong network reciprocity. As a result the C-state appears on both layers. We generally confirm that the intra-layer interactions on the layer with the larger size affect the steady-state considerably.

Discussion
In summary, we study the SEPGG on DRN. When two CG's of the same size interact through the inter-coupling with k inter = 1 and p = 1, the steady-state density, ρ α S of each strategy on each layer α can be exactly described by the mean-field theory. If the 〈 k intra 〉 α decreases then ρ α S 's on each layer slightly deviates from the mean-field expectation. Such deviation is relatively small when the multiplication factor r is small. While if 〈 k intra 〉 1 < 〈 k intra 〉 2 (asymmetric DRN), then ρ α S 's are determined by ρ α S 's on the layer 1 and the network reciprocity can be reinforced through the inter-layer coupling for 0 < p < 1. For 0 < p < 1 on the symmetric DRN, ρ α S 's show exactly the same behavior with those on the corresponding single network, and p only makes the time-delay (see Supplementary  Information). The schematic diagrams of the non-vanishing ρ α S 's are also provided in Supplementary Information. Furthermore, we also investigate the behavior of ρ α S 's on the DRN with 〈 k intra 〉 < 1 or N 1 ≠ N 2 . On the DRN with 〈 k intra 〉 < 1, we find that ρ α S 's are nearly the same as those for 〈 k intra 〉 = 1 with p = 〈 k intra 〉 . Finally, if N 1 < N 2 then the steady-state density is determined by the state of layer 2, thus the network reciprocity of the entire network can be enhanced when 〈 k intra 〉 2 is small enough.
Furthermore, since each individual in real world interacts to each other through several different channels of interactions, such interaction topology sometimes can be well described by the multiplex networks in which all the layers have the same set of nodes in general. The DRN with N 1 = N 2 and k inter = 1 could be related to a kind of multiplex networks. In addition, DRN with N 1 ≠ N 2 or k inter < 1 corresponds to the interconnected networks in which each layer has different set of nodes. Therefore, our SEPGG on DRN model would provide more general framework to study the emergence of cooperation in more realistic systems.
Finally, we want to make some remarks on important open questions. Although we do not assume any detailed topological properties, many studies have revealed that some topological properties of a network such as degree heterogeneity 22,[49][50][51][52] , degree-degree correlation 53 , and clustering coefficient 54 can significantly change the evolution of cooperation. The cost heterogeneity is also known to play a nontrivial role in the emergence of cooperation 22 . As β changes many interesting phenomena related to the phase transition has been reported 55,56 . Thus it would be very important to study how such topological properties, cost heterogeneity, and noise level affect the evolution of cooperation in interlinked networks.