Introduction

Network science has registered tremendous breakthroughs in the last years. Among these findings are the development of analytical and computational tools to study large real networks, as well as the development of more realistic models. These studies interpreted real networks mostly as single, isolated entities. However, the deeper understanding of complex networks is showing that in fact they are generally organized as networks of networks and therefore the main focus of scientific research is shifting to multiplex and interdependent networks1.

Interdependent networks are organized in two or more layers (subnetworks), in which the functioning of a node in one layer depends on neighbor nodes in other layers. Multiplex networks2,3, on the other hand, are composed of a single type of node and several types (“colors”) of links. One example of a real multiplex network is the transportation network of a city: nodes represent locations in the city and links are of different type depending on the transportation system(s) that connect two given locations (metro, autobus, etc). Interdependent networks have already been applied in studies in diverse areas of science4. There is significant progress in the understanding of their percolation properties and robustness1,5,6,7,8,9,10,11,12,13. Furthermore, interdependent networks were also applied in the evaluation of seismic risk14,15, in the emergence of creativity16 and on the impact of such population structure on voting outcome17.

Network science18,19,20, together with evolutionary game theory21, helps to understand how cooperative behavior can emerge and evolve in structured populations of selfish individuals. Cooperators contribute with a cost so that another individual can receive a certain benefit. Defectors on the other hand do not contribute, yet reap the benefit (which we assume to be larger than the cost). Interactions between individuals are traditionally modelled by one-shot symmetric two-person games, including the Prisoner's Dilemma (PD)22,23, the Snowdrift Game (SG)24 or the Stag-Hunt game (SH)25. In structured populations individuals only interact with their nearest neighbors and it becomes possible for cooperators to escape exploitation by defectors by forming clusters in which they support each other. This process is known as network reciprocity26,27,28,29,30,31,32,33.

Population structure for games was traditionally modelled as a single, isolated network. Only very recently research has started to evaluate the impact of more complex, interdependent networks on the evolution of cooperation. To the best of our knowledge, all previous works have implemented the same game in all layers34,35,36,37,38,39,40. Different layers however can represent distinct environments, characterized by their game rules, which can interact or be coupled. Consider the example of two companies and the joint network of contacts between their employees. The code of conduct in each company may be different, depending on its internal organization and goals. Nevertheless, employees from one company may interact with employees from the other company, because they know each other personally or for instance when establishing a business transaction. By interacting with an acquaintance from the other company, each employee can acquire new working strategies, which can afterwards be imitated by their contacts inside each company. Another natural example is interaction between subpopulations of different cultural backgrounds.

Here we investigate how imitation and interaction between individuals playing different games influence the final cooperation levels. For that, we consider a population organized in two layers of equal size, in which one distinct game is played in each of the layers, as schematized in Fig. 1 (see also Methods for details). Nodes establish intralinks with neighbors of the same layer and interlinks with neighbors of the opposite layer. In the following, we start by studying by computer simulations the case in which both layers are regular random graphs. Afterwards, we also consider the limit case of a well-mixed population of infinite size, also organized in two layers. We introduce bias in imitation: individuals imitate a neighbor from the same layer with probability p, or a neighbor of the opposite layer with probability (1 − p). We investigate the impact of varying p on the final level of cooperation reached in each of the layers (see Methods).

Figure 1
figure 1

Scheme of the population structure.

The population is organized in two layers. Individuals establish intralinks with neighbors of the same layer and interlinks with neighbors of the opposite layer. Depending on the layer in which they are located, individuals compute their payoff taking into account the PD payoff matrix or the SG payoff matrix, indicated on the left of the figure. On the right, we schematize the relation between the parameters of the two games, which we use in this paper. We assume that TPD = TSG = T, SSG = S and SPD = −S, with S [0, 1], T [1, 2].

Results

Regular random layers

For reference, we present in Fig. 2 the final fraction of cooperators for the PD and the SG played in isolated populations. Panel A) shows the results of computer simulations in regular random networks, while panel B) shows the stationary solutions of the replicator equation for an infinite well-mixed population, that is a population in which every individual has the same probability of interacting with anyone else41. In the replicator equation, x represents the fraction of cooperators in the population, fC(x) and fD(x) stand for the mean fitness of cooperators and defectors respectively (with fC(x) = xR + (1 − x)S and fD = xT + (1 − x)P). The PD is the harshest social dilemma for cooperation and therefore only a small region characterized by S close to 0 and T close to 1 is able to escape full defection, in regular random networks. The SG, on the other hand, is a coexistence game, in which the final fraction of cooperators x in an infinite well-mixed population is given by x = (PS)/(RST + P).

Figure 2
figure 2

Final fraction of cooperators in single, isolated populations.

The contour plot of panel A) shows the final average fraction of cooperators for the PD and SG when these social dilemmas are played in an isolated regular random network of size N = 103 and node degree k = 4. Panel B) shows the corresponding results for an infinite well-mixed population, solution of the replicator equation , where fC(x) and fD(x) stand for the mean fitness of cooperators and defectors respectively41. Red corresponds to full defection, blue corresponds to full cooperation. The parameter β is set to 1.0 for panel A).

To evaluate the impact of an interdependence between the PD and the SG on the evolution of cooperation, we start by considering a population in which each of the layers is modeled by a regular random graph with nodes of degree k = 4. This graph has no finite loops in the infinite size limit. We assume that each node establishes one interlink with a uniformly randomly chosen node from the opposite layer. That is, individuals engage in four interactions with neighbors of the same layer and one interaction with a neighbor from the other layer.

The contour plots in Fig. 3 show the result of our simulations, namely, the final fraction of cooperators in each of the layers as a function of the game parameters T and S, for several biased imitation probabilities p (0 ≤ p ≤ 1). When p = 1 individuals can only imitate neighbors from the same layer and the results obtained are qualitatively similar to those known for isolated populations (see panel A) of Fig. 2).

Figure 3
figure 3

Final fraction of cooperators when layers are modelled as random regular graphs, for various values of probability p.

Contour plots show the final average fraction of cooperators on both the PD (S < 0) and SG (S > 0) layers as a function of the game parameters T and S, starting from 50% of cooperators and defectors randomly distributed in each of the layers. The dashed line separates the results for the SG layer (upper quadrants) from the results for the PD layer (lower quadrants). Red corresponds to full defection, blue corresponds to full cooperation (x stands for the fraction of cooperators in a given layer).

When p < 1 strategy transfer between layers becomes possible. Comparing the results for the PD layer for p = 0.9, p = 0.95 and p = 1, we reach an interesting conclusion: although the fraction of cooperators is still low, even a small probability of imitating a neighbor playing a distinct social dilemma (in this case, the SG) is sufficient to prevent full defection in the majority of the PD quadrant.

As p further decreases, for a given pair of game parameters T, S the fraction of cooperators in the stationary state becomes approximately equal in both layers (note that SSG = S, SPD = −S with S [0, 1]). However, a careful inspection of the contour plots presented in Fig. 3 for 0.1 ≤ p ≤ 0.5 shows that the fraction of cooperators in both layers is not exactly identical (i.e. the contours are not exactly symmetrical): for a fixed pair of T and S values, the PD layer tends to have a slightly lower fraction of cooperators.

Finally, the results for p = 0 are surprisingly different from the ones for higher p. This contrast indicates the different structures of absorbing states at different values of p. When p = 1 four absorbing states exist and both layers reach homogeneous states: both reaching full defection (full cooperation), or one layer reaching full cooperation and the other full defection (and vice versa). When 0 < p < 1, individuals can imitate neighbors from both layers and therefore only the two absorbing states in which both layers reach the same homogeneous state are possible. In contrast, when p = 0, each individual can only imitate his neighbor connected through the interlink. As a result, there is an infinite number of absorbing states. When all interlinks connect pairs of individuals with equal strategies, evolution stops, as schematized in panel A) of Fig. 4. Since strategies and interlinks are distributed randomly in the population, on average less than N/2 interlinks will actually connect individuals of different strategies and the population will rapidly evolve to an absorbing state.

Figure 4
figure 4

Absorbing states for p = 0.

Each panel represents a snapshot of the two layers seen in profile. Panel A) schematizes one possible absorbing state when layers are modelled as regular random networks. Since each individual can only imitate his neighbor connected through the interlink, the population reaches an absorbing state organized as an arbitrary set of independent pairs (cooperator–cooperator and defector–defector) in which one member is in one layer and the other is in the other layer. Panel B) represents one of the two possible absorbing states in a well-mixed population. Since each node establishes interlinks to every node from the opposite layer, evolution ceases when both layers reach full cooperation or full defection.

Figure 5 shows the cross section of the contour plots of Fig. 3 at fixed S = 0.5 for several p values. This illustrates the different behavior of the level of cooperation in each of the layers for decreasing p. Excluding the results for p = 0 (which are a direct consequence of allowing exactly one interlink for each node) we observe that the fraction of cooperators monotonously decreases for decreasing p in the SG layer, while in the PD layer it monotonously increases. As we will see in the following, this behavior will not be maintained when this population structure is replaced by a well-mixed population.

Figure 5
figure 5

Final fraction of cooperators versus T at fixed S = 0.5, for both layers.

Results correspond to those presented in Fig. 3 when S = 0.5. The fraction of cooperators for the SG layer is shown in the upper plot, while the corresponding curves for the PD layer are presented in the bottom plot. For simplicity, we include only some of the p values studied in Fig. 3.

Well-mixed population

We now proceed to analyze the general case of a population modelled by a fully connected network organized in two layers. Each layer has N/2 nodes and intralinks. Contrary to the previous case of regular random layers, in which each node established one interlink, we now assume that each node establishes N/2 interlinks to the N/2 nodes of the opposite layer, in accordance with the structure of a fully connected network. Therefore, each node has N − 1 neighbors: N/2 − 1 in the same layer as itself and N/2 in the opposite layer. nPD and nSG stand for the number of cooperators in the PD (SG) layer, respectively. In this case, mean-field theory is exact and the evolution of population can be expressed completely by nPD (t) and nSG (t).

We obtain the solutions for the stationary state of this model in the infinite population limit, N → ∞. For that, we adopted the variables xPD = nPD/(N/2) and xSG = nSG/(N/2) (xPD [0, 1] and xSG [0, 1]). The mean fitness of individuals in each of the layers is expressed in terms of the game parameters and the fraction of cooperators and defectors. In the PD layer, the fitnesses of cooperators (fC,PD) and defectors (fD,PD) are expressed as:

Similarly, in the SG layer the fitnesses of cooperators (fC,SG) and defectors (fD,SG) are

We obtain the following rate equations for the evolution of the fraction of cooperators in the PD and SG layers:

For instance, in Eq. (5) the factor (1 − xPD) stands for the probability of selecting a defector from the PD layer. The second factor, in the square brackets, accounts for the two possible scenarios: with probability p the defector can imitate a cooperator from the PD layer, while with complementary probability (1 − p) the defector can imitate a cooperator that belongs to the SG layer. The interpretation is analogous for the second line of this equation, as well as for Eq. (6). The stationary state solution corresponds to that obtained from Eqs. (5) and (6) when and .

Figure 6 shows the final stationary states of this mean-field model for different values of probability p, provided by the numerical solutions of Eqs. (5) and (6) for β = 1.0. These plots actually explain our simulations in the previous section for coupled regular random networks (see Fig. 3). One can clearly see that Figs. 3 and 6 show qualitative similarities, with only one dramatic exception, p = 0, which will be discussed below. In particular, in the case of p = 1 both Figs. 3 and 6, as well as Figs. 5 and 7, are close to the corresponding results for isolated populations (random regular and well-mixed, respectively, compare with Fig. 2).

Figure 6
figure 6

Final fraction of cooperators on each layer of the well-mixed population, for various values of probability p.

Contour plots show the final fraction of cooperators on both the PD (S < 0) and SG (S > 0) layers (solution of Eqs. (5) and (6)) as a function of the game parameters T and S. The dashed line separates the results for the SG layer (upper quadrants) from the results for the PD layer (lower quadrants). Red corresponds to full defection, blue corresponds to full cooperation (x stands for the fraction of cooperators in a given layer). Parameters: β = 1.0.

Figure 7
figure 7

Final fraction of cooperators versus T at fixed S = 0.5 for both layers of a fully connected network.

Results correspond to those presented in Fig. 6 when S = 0.5. The fraction of cooperators for the SG layer is shown in the upper plot, while the corresponding curves for the PD layer are shown in the bottom plot.

One can see from Fig. 6 that a small probability of imitating a neighbor from the SG layer is beneficial for the levels of cooperation reached in the PD layer. The final fraction of cooperators in each of the layers becomes approximately equal as p decreases, yielding almost symmetrical contours for 0.1 ≤ p ≤ 0.5. However, the fraction of cooperators is in general slightly lower in the PD layer than in the SG layer for a fixed pair of game parameters T and S.

Note some difference between Figs. 3 and 6 at sufficiently low, still non-zero, p. With diminishing p, the region of full defection in coupled regular random populations (Fig. 3) approaches some limit. In contrast, for decreasing p, coupled well-mixed populations (Fig. 6) do not show any area with full defection, but the region with high level of defection becomes wider and wider as p decreases. Finally, at p = 0, the full defection state is realized in the entire TS plane. This difference is also observed between Figs. 5 and 7. Especially it manifests itself in two different behaviors of fractions of cooperators in the PD layer in Figs. 5 and 7, monotonous with p and non-monotonous, respectively.

At p = 0 the difference between the coupled regular random networks and well-mixed population is dramatic. While for the coupled well-mixed populations, the final state is full defection for the whole TS plane (Fig. 6), for regular random networks we observed coexistence of defection and cooperation for the whole TS plane. The reason for this controversy is the fundamental difference between the structures of the absorbing states in these two systems. Clearly, in the mean-field system (well-mixed) there are two absorbing states: full cooperation and full defection. For p = 0, the full defection scenario is realized. In contrast, in the coupled regular random populations, at p = 0 the system has an infinite number of absorbing states. We already explained that this state is organized as an arbitrary set of independent pairs (cooperator-cooperator and defector-defector) in which one member is in one layer and the other is in the other layer. The final fraction of cooperators and defectors in the population for a given point on the TS plane is determined by a specific initial condition.

It is interesting to relate our findings with observations in other works39,40 which studied the role of biased fitness (utility) on the level of cooperation in interdependent structured populations. Although these works explored a bias and network architectures very different from our study, they reported the existence of an optimal control parameter value for cooperation, similarly to what we notice for a well-mixed population.

We can also relate our model to the multi-network model42,43, where links between nodes are of two sorts: the interaction links through which individuals play and the imitation links through which individuals revise their strategies. All links in our model are used both for interaction and imitation. The principal difference is that interaction occurs with probability 1 for every link, while imitation takes place with probability p along intralinks and the complementary probability 1 − p along interlinks.

Discussion

We explored the evolution of cooperation in two games—the Prisoner's Dilemma and the Snowdrift Game—played within an interdependent network. We demonstrated that as the probability to imitate neighbors from the opposite layer increases, the final level of cooperation in each of the layers behaves differently. A small probability 1 − p is always detrimental for the level of cooperation in the SG layer, regardless of the population structure (well-mixed or when layers are modelled as regular random networks). In the case of a well-mixed population we observed an intermediate optimal value of p for which the fraction of cooperators in the steady state is maximal in the PD quadrant. On the other hand, when layers are regular random graphs, the level of cooperation in the PD layer registers a significant increase from full defection as p decreases from 1, as is shown in Fig. 3. We have also shown that different arrangements of interconnections between layers result in fundamentally different structures of absorbing states for p = 0, that is the limit case in which individuals can only imitate neighbors through the interlinks.

The constraint between the game parameters of the PD and the SG implemented in this work is one of possible options. Based on our mean-field approach we have verified that our conclusions remain qualitatively valid for other relationships between game parameters44. In this work we studied a representative case of β = 1. In the Supplementary Information we consider other values of β.

The Prisoner's Dilemma and the Snowdrift Game are two of the possible two-person games. Future research directions include establishing an interdependence between other two- and N-person games, as well as extending the biased imitation process here reported to the case in which layers are strongly heterogeneous, such as scale-free networks45.

Methods

Evolutionary games

We implement standard two-person games. Each individual can act as an unconditional cooperator (C) or as an unconditional defector (D). The possible payoffs resulting from the interaction between two players are resumed in the payoff matrix

in which the entries represent the payoff earned by the row player. When both individuals opt for cooperation (defection), they both receive a Reward, R (Punishment, P). When a cooperator meets a defector, the cooperative player received the Sucker's payoff (S), while the defective player receives the Temptation (T) for defecting.

The relative ordering of the game parameters R, S, T and P defines four different two-person games46. In this study, we consider the Prisoner's Dilemma (PD) (characterized by T > R > P > S)23 and the Snowdrift Game (SG) (characterized by T > R > S > P)24. We fix R = 1, P = 0 as is conventional26.

Population structure

We consider a population organized in two layers (schematized in Fig. 1), in which each node represents a different individual. In each layer all individuals play the same game, the PD in one layer and the SG in the other. We refer to links between nodes of opposite layers as interlinks and links between nodes of the same layer as intralinks. In the absence of interlinks, each layer stands as an isolated population in which one single game is played. When layers are modelled as regular random graphs, we assume that each node establishes one interlink to a randomly chosen node from the opposite layer. On the other hand, when studying a two-layer fully connected network, we assume that each node establishes interlinks to each of the nodes in the opposite layer.

Relation between game parameters

We assume that individuals in the PD layer adopt the game parameters RPD, SPD, TPD and PPD (with TPD > RPD > PPD > SPD). Similarly, individuals in the SG layer consider RSG, SSG, TSG and PSG (with TSG > RSG > SSG > PSG). The game parameters for the PD and the SG can be related in several ways. We opt for RPD = RSG = 1, PPD = PSG = 0 and TPD = TSG = T with T [1, 2]. Furthermore, we constrained SSG = S and SPD = −S, with S [0, 1]. The resulting payoff matrix for the SG is

while the payoff matrix for the PD layer is given by

The difference between matrix (8) for the SG layer and matrix (9) for the PD layer is in the payoff obtained by a cooperator when interacting with a defector, the sucker's payoff. When interacting with a defector, a cooperator located in the SG layer earns S, while a cooperator located in the PD layer earns −S. We stress that this is only one of possible relations between the parameters of the two games.

Individuals play the same game with all neighbors, including the neighbor(s) from the other layer. An individual's fitness represents a measure of his success and corresponds to the payoff accumulated over all interactions in which he engages in one time-step. Through the interlinks, individuals obtain their payoff by taking into account the strategy of their co-player and the payoff matrix used in their own layer.

Update dynamics

Individuals revise their strategies through a biased imitation process: with probability p they imitate a neighbor from the same layer, while with probability (1 − p) individuals imitate a neighbor from the opposite layer. Update is asynchronous: at each time-step we randomly choose one of the layers and from it one randomly chosen individual, A, to revise his strategy. Then we choose with probability p a neighbor from the same layer and with complimentary probability (1 − p), a neighbor from the opposite layer. If A and his chosen neighbor, B, have different strategies, A imitates the strategy of B with a probability F (fBfA) proportional to their fitness difference and given by the Fermi distribution47,48,

fA (fB) stands for the fitness of A (B) and β represents the intensity of selection, which regulates the accuracy of the imitation process. In the limit β → 0 evolution proceeds by random drift, while β → ∞ corresponds to the deterministic limit of pure copying dynamics: A imitates B iff fB > fA, with probability 1. We consider the respresentative case of β = 1.0.

Simulations

The results of Figs. 3 and 5 were obtained for a population of total size N = 2000, in which both layers had an equal size, N/2 = 1000. Each point in the contour plots corresponds to an average over 103 trials. In each trial, we start by randomly selecting one over 102 realizations of a regular random network to model each layer and establish N/2 interlinks between the layers. Each interlink is established by uniformly randomly choosing two end nodes, one from each layer, allowing each node to establish exactly one interlink. Each trial starts by randomly distributing 50% of cooperators and defectors in each of the layers. We allow the population to evolve for 104 generations (where 1 generation corresponds to N/2 strategy revisions), after which we average the fraction of cooperators over the next 103 generations.