Dynamics of the model and the phase diagram for n = 2 public resources

As shown in the section Materials and methods, the model can be solved analytically in terms of the replicator-mutator equation. Analysis of the model by simulations and numerical solution of the replicator dynamics reveals that the dynamics can settle in a fixed point or a periodic orbit depending on the parameter values and the initial conditions. We begin by studying cyclic behavior. In Fig 1A, we plot the density of cooperators (solid blue) and defectors (dashed red) in the game i, denoted as, respectively, $$ {\rho }_C^i$$ and $$ {\rho }_D^i$$, for i = 1 and i = 2, as a function of time. The density of cooperators (defectors) in the game i is defined as the fraction of individuals in the population who prefer public resource i and cooperate (defect). In Fig 1B the total density of cooperators $$ {\rho }_C={\rho }_C^1+{\rho }_C^2$$ (solid blue), and the density of the individuals who prefer public resource 1, $$ {\rho }^1={\rho }_C^1+{\rho }_D^1$$, (dashed red) are plotted. The total density of defectors, $$ {\rho }_D={\rho }_D^1+{\rho }_D^2$$, and the density of individuals who prefer public resource 2, $$ {\rho }^2={\rho }_C^2+{\rho }_D^2$$, can be inferred from these by noting ρ_D = 1 − ρ_C, and ρ² = 1 − ρ¹. Here, g = 10, r₁ = 4.4, r₂ = 2.8, ν = 10⁻³, and the initial condition is a random initial condition in which both the strategy and the preferred game of the individuals are assigned at random. Here, and in the following, we fix c = π₀ = 1. The top panels show numerical solutions of the replicator dynamics, and the bottom panels show the results of a simulation in a population of size N = 10000.

Fig 1

Dynamics of the model for n = 2 public resources.

A: The density of cooperators (solid blue), $$ {\rho }_C^i$$, and defectors (dashed red), $$ {\rho }_D^i$$, in a game i, for i = 1 and 2, as a function of time. The two top panels result from the replicator dynamics, and the two bottom panels result from a simulation. B: The total density of cooperators, $$ {\rho }_C={\rho }_C^1+{\rho }_C^2$$ (solid blue), and the density of the individuals who prefer public resource 1, $$ {\rho }^1={\rho }_C^1+{\rho }_D^1$$ (dashed red), as a function of time. The top panel results from replicator dynamics and, the bottom panel results from a simulation. Parameter values: g = 10, ν = 10⁻³, r₁ = 4.4, and r₂ = 2.8. Simulations are performed in a population of size N = 10000.

A comparison shows a high level of agreement between simulations in a finite population and solutions of replicator dynamics, which is an exact solution of the model in the infinite population limit. Both simulation and analytical solutions show that cooperators are preserved and cyclically dominate the population. When cooperators’ density in a public resource i increases, defectors in the game i obtain a high payoff and increase in frequency. This decreases the payoff of cooperators in the game i, and consequently, they obtain a higher payoff by switching to the other public resource. This, in turn, decreases the payoff of defectors in the game i, and thus $$ {\rho }_D^i$$ starts to decrease when $$ {\rho }_C^i$$ decreases enough.

As mentioned before, cyclic behavior is not the only attractor of the dynamics. To see this, in Fig 2A and 2B, the phase diagrams of the system in the r₁−r₂ plane for g = 10, and respectively, ν = 10⁻⁵ and ν = 10⁻³ are plotted. For too small enhancement factors, having a choice between public resources can not alleviate the strong disadvantage cooperators experience due to a small enhancement factor. Consequently, the dynamics settle into a fixed point in which defectors dominate the population. As either enhancement factors increase beyond a threshold (marked by blue triangles), cooperators can survive. In this region, the dynamics become bistable, and depending on the initial condition, settle either in a defective fixed point or a periodic orbit. By further increasing the enhancement factors, at a second bifurcation line (denoted by red crosses), the defective fixed point loses stability. The dynamics become mono-stable and settle in a periodic orbit, starting from all the initial conditions. The phase boundary between the defective and cyclic phases in the r₁ − r₂ plane can be defined as the line where the transition between the two phases occurs starting from a homogeneous initial condition (that is $$ {\rho }_C^1={\rho }_D^1={\rho }_C^2={\rho }_D^2=0.25$$) [39]. Green circles in the figure denote this boundary.

Fig 2

The phase diagram of the model for n = 2 public resources.

The phase diagram of the model for ν = 10⁻⁵, A, and ν = 10⁻³, B. The model shows two distinct mono-stable phases where the dynamics settle into a fixed point (denoted by FP) or a periodic orbit (denoted by PO). In between, for small enhancement factors, the system possesses a bi-stable region where both the periodic orbit and the fixed point are stable (the region between blue (triangles) and red (crosses) lines, denoted by PO+FP). The phase boundary (green line marked by circles) can be defined as the line where the transition between the two phases occurs starting from a random assignment of strategies (that is $$ {\rho }_C^1={\rho }_C^2={\rho }_D^1={\rho }_D^2=0.25$$). Here, g = 10 and the replicator dynamics is used to derive the phase diagrams.

By further increasing the enhancement factors, beyond a final transition line (the line where different marker types coalesce), the periodic orbit becomes unstable, and the dynamics settle in a fixed point in which cooperators dominate. We note that in contrast to the transition from the defective fixed point to the cyclic orbit, the transition between the cyclic orbit and the cooperative fixed point shows no bi-stability and happens continuously (in the sense that the amplitude of fluctuations decreases gradually by increasing the enhancement factors until fluctuations vanish at the transition line). Consequently, this transition line occurs at the same value of r₁ and r₂ for all initial conditions. A comparison of the phase diagram for different mutation rates shows that the same scenario occurs for all the mutation rates. Quantitatively, however, larger mutation rates decrease the size of the region in the r₁ − r₂ plane where cyclic behavior can occur.

Cooperation level of the population for n = 2 public resources

An interesting question is how the level of cooperation in the population depends on the enhancement factors? To see this, in Fig 3A, we plot the time average of the total density of cooperators, $$ {\rho }_C={\rho }_C^1+{\rho }_C^2$$, in the r₁ − r₂ plane. In all panels of Fig 3, g = 10, ν = 10⁻³, and a homogeneous initial condition is used. The top panels result from the replicator dynamics, and the bottom panels show the result of a simulation in a population of size N = 20000. As can be seen, the density of cooperators is the highest close to the diagonal, where r₁ = r₂, and decreases with increasing the distance from the diagonal. This shows that increasing the enhancement factor of one of the games beyond its optimal value can have a surprisingly detrimental effect on the evolution of cooperation.

Fig 3

Cooperation level of the population in the model with n = 2 public resources.

A: The contour plot of the time average density of cooperators, $$ {\langle {\rho }_C\rangle }_t={\langle {\rho }_C^2+{\rho }_C^1\rangle }_t$$. Cooperation evolves for large enough enhancement factors, and it is maximized when the enhancement factors of the two resources are similar. B and C: The time average density of cooperators B, and defectors C, in resource 1. Compared to cooperators, defectors are more likely to prefer the highest quality resource. Consequently, cooperators can use the low-quality resource as a shelter when defectors dominate. D: Contour plot of $$ \gamma ={\langle {\rho }_C^1/({\rho }_C^1+{\rho }_C^2)-{\rho }_D^1/({\rho }_D^1+{\rho }_D^2)\rangle }_t$$. This is a measure of how much the cooperators are more likely to prefer resource 1, compared to defectors. Cooperators are more likely to prefer resource 1, when it is the lower quality resource, and defectors are more likely to prefer resource 1 when it is the higher quality resource. Here, g = 10 and ν = 10⁻³. Top panels result from the replicator dynamics, and bottom panels show the result of a simulation in a population of size N = 20000.

To shed more light on the mechanism by which freedom to choose between different public resources promotes cooperation, and to see why cooperation is maximized when the two competing resources have similar enhancement factors, in Fig 3B and 3C, we plot the time average density of, respectively, cooperators and defectors in resource 1 (Due to the symmetry of the two games, the corresponding densities in resource 2 result from this by reflection with respect to the diagonal). As can be seen, most individuals prefer resource 1, as long as it has a higher enhancement factor. Interestingly, the low-quality resource is often immune to invasion by defectors: while the inferior resource can attract cooperators, defectors strongly prefer the higher quality resource. This phenomenon can be made more quantitative by calculating the fraction of cooperators who prefer resource 1 as $$ {\rho }_C^1/({\rho }_C^1+{\rho }_C^2)$$, and the fraction of defectors who prefer resource 1, $$ {\rho }_D^1/({\rho }_D^1+{\rho }_D^2)$$. These quantities measure how likely a cooperator or a defector is to prefer resource 1. The time average difference between these two quantities, $$ \gamma ={\langle {\rho }_C^1/({\rho }_C^1+{\rho }_C^2)-{\rho }_D^1/({\rho }_D^1+{\rho }_D^2)\rangle }_t$$ is plotted in Fig 3D. As can be seen, cooperators are always more likely to prefer the lower quality resource than the defectors are. On the other hand, defectors are more likely to prefer the higher quality resource than the cooperators are. This phenomenon is crucial for freedom to choose between the two resources to promote cooperation: When the density of defectors increases in the population, cooperators can shelter in the lower quality resource and make it the more profitable one, and thus, by enjoying a higher growth rate, out-compete defectors. However, when the inferior resource’s quality is much lower than that of the superior resource, it becomes more difficult for cooperators to make the inferior resource the more profitable one. This, in turn, leads to a reduction in the cooperation level of the population by increasing the difference between the quality of the resources.

The case of n = 3 public resources

As mentioned in the section The model, the framework can be extended to a general case where individuals can choose among an arbitrary number, n, of public resources. In this section, we study the case of n = 3 public resources. The phase diagram of the model with three resources is plotted in Fig 4A and 4B, for respectively, ν = 10⁻⁵ and ν = 10⁻³. Here, the replicator dynamics is used to derive the phase diagrams in the r₁ − r₂ plane, fixing the enhancement factor of game 3 to r₃ = 3, and setting g = 10. The model shows a similar phase diagram to that observed for the case of n = 2 resources: For small enhancement factors, the dynamics settle into a defective fixed point. As the enhancement factors increase, the system becomes bistable and can settle in either the defective fixed point or a periodic orbit in which cooperators survive and cyclically dominate the population. The blue line (marked by blue triangles) shows the lower boundary of the coexistence region above which the periodic orbit becomes stable. The red line (marked by crosses) shows the upper boundary of the coexistence region, above which the defective fixed point becomes unstable, and the dynamics settle in the periodic orbit. In between the two bifurcation lines and in the coexistence region, both the defective fixed point and the periodic orbit are stable, and the dynamics settle into one of these attractors depending on the initial conditions. The phase boundary, defined as the boundary above which the transition from the fixed point to periodic orbit occurs starting from a homogeneous initial condition [39], is indicated by the green line (marked by circles). Finally, for large enhancement factors, the system shows a final bifurcation line above which the dynamics settle in a cooperative fixed point, in which the vast majority of the individuals cooperate. This boundary is indicated by the line in large enhancement factors where different marker types coincide.

Fig 4

The phase diagram of the model with n = 3 public resources.

The phase diagram of the model with n = 3 public resources, for ν = 10⁻⁵, A, and ν = 10⁻³, B. For small enhancement factors, the model shows bistability and can settle in either a fixed point (FP) where defectors dominate or a periodic orbit (PO) in which cooperators survive. The blue line (marked by triangles) shows the lower boundary of the coexistence region above which the cyclic orbit becomes stable, and the red line (marked by crosses) shows the upper boundary of the coexistence region, above which the defective fixed point becomes unstable. The green line (marked by circles) shows the phase boundary resulted from an unbiased (homogeneous) initial condition. Here, g = 10, r₃ = 3, and the replicator dynamics are used to derive the phase diagrams.

To see how the cooperation level of the population changes in different regions of the phase diagram, in Fig 5A and 5B, we plot the time average population’s cooperation level, $$ {\langle {\rho }_C\rangle }_t={\langle {\rho }_C^1+{\rho }_C^2+{\rho }_C^3\rangle }_t$$. Here, and in Fig 6 in the following, ν = 10⁻³, g = 10, and r₃ = 3. The initial condition is a homogeneous initial condition in which the individuals’ strategies and preferred games are randomly assigned. Fig 5A shows the solutions of the replicator dynamics, and Fig 5B shows the result of simulations in a population of N = 10000 individuals. As was the case in the model with n = 2 resources, the cooperation level of the population is maximized close to the diagonal, where the two higher quality games have similar enhancement factors and decreases by going away from the diagonal.

Fig 5

Cooperation level of the population in the model with n = 3 public resources.

The contour plots of the time average cooperation level of the population, $$ {\langle {\rho }_C\rangle }_t={\langle {\rho }_C^1+{\rho }_C^2+{\rho }_C^3\rangle }_t$$, resulting from numerical solutions of the replicator dynamics with n = 3 games, A, and a simulation in population of size N = 10000, B. The population’s cooperation level is the highest when the two highest quality resources have similar enhancement factors. Here, g = 10, r₃ = 3, and ν = 10⁻³. The initial condition is a homogeneous initial condition in which the strategies and the preferred game of the individuals are assigned at random.

Fig 6

The densities of different strategies in the model with n = 3 public resources.

The time average density of cooperators in resource 1 A, the time average density of defectors in resource 1 B, the time average density of cooperators in resource 3 C, and the time average density of defectors in resource 3 D, are plotted. The top panels show the replicator dynamics results, and the bottom panels show the result of a simulation in a population of size N = 10000. Here, g = 10, r₃ = 3, and ν = 10⁻³. The initial condition is a homogeneous initial condition in which the strategies and the preferred games of the individuals are assigned at random.

To take a deeper look into the model’s behavior for n = 3 games, we plot the time average density of cooperators in resource 1, $$ {\langle {\rho }_C^1\rangle }_t$$, and that of defectors in resource 1, $$ {\langle {\rho }_D^1\rangle }_t$$, in respectively, Fig 6A and 6B. The time average density of cooperators in resource 3, $$ {\langle {\rho }_C^3\rangle }_t$$, and the time average density of defectors in resource 3, $$ {\langle {\rho }_D^3\rangle }_t$$ are plotted in, respectively, Fig 6C and 6D. The top panels show the result of the replicator dynamics, and the bottom panels show the result of simulations in a population of size N = 10000. We note that, due to the symmetry of resource 1 and 2, the densities of cooperators and defectors in resource 2 result from that in resource 1, in Fig 6, by reflection with respect to the diagonal.

In the region where both r₁ and r₂ are larger than r₃ = 3, that is when resource 3 is the most inferior resource, the third resource becomes almost obsolete: Very few individuals prefer resource 3, and the dynamic is determined by the interaction of the two higher-quality resources. However, the mere possibility of playing the third game has a positive effect on cooperation, even when it does not attract individuals significantly. We note that, in the case r₁ = r₂, it can happen that cooperators, but not defectors, choose resource 3, provided its enhancement factor is not too small compared to the enhancement factors of the two higher-quality resources. This, once again, shows cooperators are more likely to choose an inferior resource.

On the other hand, in the region where either r₁ or r₂ are smaller than r₃ = 3, that is when resource 3 is the second high-quality resource, the third resource always attracts cooperators, but not defectors (compare Fig 6C and 6D). This observation shows, similarly to the case of n = 2 resources, when individuals can choose between three resources, defectors tend to prefer the highest quality resource. This makes a second high-quality resource a shelter for cooperators to work cooperatively and grow safely when the density of defectors increases in the population. The lowest quality resource, on the other hand, remains almost obsolete (as long as its quality is sufficiently lower than the other two resources). This is because the second high-quality resource, having a low fraction of defectors, is already a good enough shelter for cooperators, and they do not need to switch to the lowest quality resource.

The replicator dynamics

In the limit of infinite population size, the dynamics of the model can be described in terms of the replicator-mutator equation [59]. This equation describes the evolution of the densities of different strategies in an evolutionary context where the strategies grow proportional to their fitness and are subject to mutations. In the general case that individuals can choose among n different public resources, the replicator-mutator dynamics for the model can be written as follows:

The simulations and methods

Simulations used in Figs 3, 5, and 6 are run for T = 5000 time steps, and the time averages are taken over the last 4000 time steps. For the solutions of the replicator dynamics presented in Figs 3, 5, and 6, the replicator dynamics is numerically solved for T = 4000 time steps, and the averages are taken in the time interval t = 2000 to t = 4000. To derive the phase diagram presented in Figs 2 and 4, the following procedure is used. To determine the boundaries of the bistable region where both defective fixed point and cyclic behavior are possible, we note that the most favorable initial condition for the evolution of cooperation is the one in which all the individuals prefer the same resource (for example, $$ {\rho }_D^1=1$$, $$ {\rho }_C^1={\rho }_C^2={\rho }_D^2=0$$ for n = 2 and $$ {\rho }_D^1=1$$, $$ {\rho }_C^1={\rho }_C^2={\rho }_D^2={\rho }_C^3={\rho }_D^3=0$$ for n = 3 resources). Numerical solutions of the replicator dynamics starting from this initial condition give the lower boundary of the coexistence region (the blue line marked by triangles in the figures). The least favorable initial condition for the evolution of cooperation is the one in which all the individuals are defectors and prefer the two games in equal fractions ($$ {\rho }_C^1={\rho }_C^2=0$$ and $$ {\rho }_D^1={\rho }_D^2=0.5$$ for n = 2 and $$ {\rho }_C^1={\rho }_C^2={\rho }_C^3=0$$ and $$ {\rho }_D^1={\rho }_D^2={\rho }_D^3=1/3$$ for n = 3 resources). Numerical solutions of the replicator dynamics starting from this initial condition give the upper boundary of the coexistence region (the red line marked by crosses). Finally, to determine the phase boundary between the defective and cyclic phases (green line marked by circles), the replicator equation with a homogeneous initial condition ($$ {\rho }_C^1={\rho }_C^2={\rho }_D^2={\rho }_D^1=0.25$$ for n = 2, and $$ {\rho }_C^1={\rho }_C^2={\rho }_C^3={\rho }_D^3={\rho }_D^2={\rho }_D^1=1/6$$ for n = 3 resources) is solved. All the three sets are solved for T = 4000 time steps, and the values of r₁ and r₂ where a periodic orbit appears or disappears mark the phase boundaries. As the phase boundary between the cyclic and cooperative phases for large enhancement factors shows no bi-stability, solutions starting from all the initial conditions result in the same line where different types of markers coalesce in the figure.