Evolution of Cooperation in the Presence of Higher-Order Interactions: From Networks to Hypergraphs

Abstract

Many real systems are strongly characterized by collective cooperative phenomena whose existence and properties still need a satisfactory explanation. Coherently with their collective nature, they call for new and more accurate descriptions going beyond pairwise models, such as graphs, in which all the interactions are considered as involving only two individuals at a time. Hypergraphs respond to this need, providing a mathematical representation of a system allowing from pairs to larger groups. In this work, through the use of different hypergraphs, we study how group interactions influence the evolution of cooperation in a structured population, by analyzing the evolutionary dynamics of the public goods game. Here we show that, likewise to network reciprocity, group interactions also promote cooperation. More importantly, by means of an invasion analysis in which the conditions for a strategy to survive are studied, we show how, in heterogeneously-structured populations, reciprocity among players is expected to grow with the increasing of the order of the interactions. This is due to the heterogeneity of connections and, particularly, to the presence of individuals standing out as hubs in the population. Our analysis represents a first step towards the study of evolutionary dynamics through higher-order interactions, and gives insights into why cooperation in heterogeneous higher-order structures is enhanced. Lastly, it also gives clues about the co-existence of cooperative and non-cooperative behaviors related to the structural properties of the interaction patterns.

1. Introduction

In behavioral terms, cooperation is the providing of a benefit to another individual at some cost for the provider. In well-mixed populations, where all individuals interact with each other, evolutionary game theory predicts that cooperation cannot survive due to the existence of selfish behaviors [1,2]. Nevertheless, we do observe cooperation in many different real systems, ranging from genomes to human societies [3,4,5,6], and not as a marginal phenomenon: rather, it often strongly shapes and makes even possible the existence of those systems. For example, without cooperative behaviors like the work division and the establishment and respect of social norms, the humankind could not have developed even the simplest forms of society. Among microbes, some viruses and cells have been observed to be involved in (involuntary) cooperative behaviors that are essential for their reproduction and diffusion [4,6]. The existence of cooperation in biological and social systems has led to hypothesize several mechanisms to explain its ubiquity: inclusive fitness, direct and indirect reciprocity, punishment, cultural group selection, etc. (clearly, those requiring substantial cognitive demands, apply only to human’s social systems) [7,8,9,10,11,12,13]. Among those mechanisms, we focus on the so-called network reciprocity [14,15,16,17,18,19,20,21,22,23,24,25,26].

Generally speaking, an act of reciprocity takes place when some individual returns a beneficial act received by someone else, establishing a fruitful mutual exchange. Therefore, it is a cooperative act. In particular, one refers to network reciprocity when the mutual cooperative exchange is sustained by the existence of a suitable structure of interactions, represented through a network, among the individuals of a population. Two conditions lie at the base of network reciprocity: a limited, recurrent set of interacting game opponents (thus avoiding the anonymity between them); and a local adaptation mechanism through which an individual can change (update) its behavior (strategy) depending on how its neighbors in the network behave and perform. Apart from its fascinating theoretical character, network reciprocity is regarded as one of the most important and interesting of the mechanisms from an application point of view. Indeed, the easiness of its requirements to be fulfilled, makes it a transverse mechanism applicable not only to those system made of rational units (e.g., social, economical), but also to those whose units are capable of only very simple behaviors (e.g., many biological systems) [22].

The use of networks to formalize the topology of interactions within a system, however, tacitly assumes that all the interactions that occur in the system are essentially pairwise. This may not be a satisfactory description for those real systems exhibiting higher-order interactions to which more than two entities participate together and which are not reducible to combinations of lower-order ones. Examples of such interactions include triadic [27] and higher-order closures in social networks [28], co-authorship networks in science [29], spatial coexistence relations among species in an ecosystem [30], and trigenic interactions in genetic networks [31]. In all such cases, representing the interaction patterns through a network (that is, approximating those higher-order interactions with pairwise ones) implies the loss of essential information to understand their dynamics [32].

This work tries to bring network reciprocity closer to reality. We investigate the effect of considering higher-order interactions, using hypergraphs, on the evolution of cooperation when the population of individuals possesses structural correlations. Hypergraphs allow to distinguish the dynamics on a clique of m individuals (an m-clique), in which each of the individuals (vertices) interacts separately, in pair, with each of the other m−1 individuals, from the dynamics in which the m individuals interact all together, as a group, in an intrinsically different way (represented as a hyperedge of cardinality m). See Figure 1 for illustration.

Our work is certainly not the first one devoted to the study of cooperation in the presence of group interactions. As reported in the review by Perc et al. [25], many efforts have been made in this direction, enlightening several ways in which the increase of the number of participants to the interactions can lead to a better sustenance of cooperation. Interestingly, the heterogeneity of the underlying network (also when coupled to the distribution of payoffs and/or contributions of players within a group) turns out to be usually a very effective feature for cooperation to thrive. All these works, however, are mainly based on pairwise structures, and thus are unable to solve the lack of uniqueness in defining group interactions from pairwise information only. When higher-order data are available, instead, higher-order structures such as hypergraphs (or, equivalently, bipartite graphs [20]) represent the data without suffering any loss of information. Nevertheless, as minutely explained in Section 6.1, our approach, allowing to smoothly transition from a pairwise structure to a higher-order one, is the first one enabling to investigate the effect of group interactions in a systematic way.

We show here that the mechanism of network reciprocity can be successfully extended to higher-order structures and that, in fact, it even becomes reinforced.

2. Dynamics

2.1. Public Goods Dilemma

To study cooperation in presence of higher-order interactions we need a game definable for any number of players (the size of the group) and belonging to the class of social dilemmas. In such games, the equilibrium solution under unilateral changes of strategy (i.e., considering changes of strategy of only one player at a time), called Nash equilibrium, is not the most efficient solution for all the players (i.e., it is not Pareto-efficient). In particular, in the type of game we are interested in, the equilibrium solution is the one in which all the players choose to defect, although the Pareto-efficient solution is the one in which they all cooperate, hence the social dilemma. The prototypical m-players dilemma of this kind is the public goods game (PGG). Each player, independently, chooses whether to put (cooperating) or not (defecting) a certain amount b (fixed for all players, in our setup) into a public pot. The collected sum is re-scaled by a synergy factor $\alpha \in \mathbb{R}$, $\alpha >1$, and equally redistributed among the m players disregarding of their strategy. When $n_c\le m-1$ neighbors cooperate, a player with strategy $\sigma$ (set to 1 if cooperates, 0 if defects) gets a payoff $f\left(\sigma \right)$ given by

$$f\left(\sigma \right)=b\left[\left(\frac{\alpha }{m}-1\right)\sigma +\frac{\alpha \left(m-1\right)}{m}q\right]$$

(1)

where $q\equiv n_c/(m-1)$ is the fraction of neighbors that cooperate. The first term accounts for the cost and proportional gain from the contribution to the pot of the player, while the second corresponds to the gain received from the cooperating neighbors. Since the second term is always non-negative, the best response for a player depends on whether $\alpha$ is smaller or greater than m. The dilemma exists for $\alpha <m$ because, in this case, according to Equation (1), the best strategy for any player is to defect. However, each of them gets a null payoff, whereas it would be positive if all of them chose to cooperate.

2.2. Structural Connectivity

We give here a brief introduction to the connectivity structures we use in this work: simple hypergraphs. Given a finite set V of vertices and a family E of subsets ${\left\{e_I\right\}}_{I\subseteq V}$ of V, the couple $H=(V,E)$ defines a hypergraph, and each $e_I$ is called hyperedge, representing a relation among all the vertices in it. The maximum and minimum cardinality (degree) of the hyperedges in H are called the rank and co-rank of the hypergraph, respectively: $\mathrm{rank}\left(H\right)=m_{max}$ and $\mathrm{co}-\mathrm{rank}\left(H\right)=m_{min}$. If they both are equal, with value $s\in \mathbb{N}$, the hypergraph is said to be s-uniform [33].

Another important concept is the 2-section of H. We can define it as the graph whose vertices are the vertices of H, and where two distinct vertices form an edge if and only if they belong to, at least, one common hyperedge in H [33] (see Figure 1).

Additionally, a hypergraph H is called simple if $e_I\subseteq e_J\Rightarrow I=J$. This means that it has no repeated hyperedges, and can be obtained from a generic hypergraph $(V,E)$ removing from E all the hyperedges which are subsets of at least another hyperedge.

Finally, given the family $S\left(i\right)$ of hyperedges ${\left\{e_J\right\}}_{J\subseteq V}$ containing vertex i, and assuming no repeated hyperedge, we define the m-degree $k_m^{\left(i\right)}$ of i as the number of hyperedges of degree m in $S\left(i\right)$. The generalized degree of node i is the number of hyperedges incident on it, given by

$$k^{\left(i\right)}=\sum _{m=m_{min}}^{m_{max}}{k}_m^{\left(i\right)}$$

(2)

We construct a hypergraph starting from a certain base network and converting in hyperedges a randomly selected fraction p of its 3-cliques, see Section 6.1 for further details. From now on, we will call these selected hyperedges as triangles.

2.3. Public Goods Dynamics in Hypergraphs

At the base of our model lie the following microscopic rules, imposed at each time step (round): (i) each player plays with its strategy, either cooperate or defect, within each of the hyperedges incident on it; (ii) each strategist chooses one of its neighbors uniformly at random in the 2-section of the hypergraph, and performs the strategy update. These steps are performed synchronously by all the players.

The payoff $f^{\left(i\right)}\left(t\right)$ that player i with m-degree $k_m^{\left(i\right)}$ and playing with strategy ${\sigma }_i\left(t\right)$ got at the end of a time step t, obeys the following dynamic equation:

$$f^{\left(i\right)}\left(t\right)=b\sum _{m=m_{min}}^{m_{max}}{k}_m^{\left(i\right)}\left[\left(\frac{\alpha }{m}-1\right){\sigma }_i\left(t\right)+\frac{\alpha (m-1)}{m}{q}_m^{\left(i\right)}\left(t\right)\right]$$

(3)

where $q_m^{\left(i\right)}\left(t\right)$ is the probability that a neighbor of node i, in a hyperedge of degree m, cooperates at time t.

The players’ perfect rationality (and this being common knowledge among them) assumed in classical game theory, according to which a player always chooses the best response to its opponents’ strategies, results to be often unrealistic. Indeed, individuals could have non-negligible cognitive limitations or could have non-consistent preferences, gathering information about possible outcomes and payoffs could be costly, or the common knowledge of the players’ rationality could fail to hold. All these limiting factors can be collected under the term bounded rationality. In an evolutionary process, this means that what drives selection is often not just the game-based fitness, but also other, uncontrollable factors intervene significantly, making selection weaker and noisier. We therefore consider smoothed imitation. To this end, we take as probability function for the strategy update, the Fermi distribution function:

$$\mathcal{F}(x;\beta )=\frac{1}{1+e^{-\beta x}}$$

(4)

where, when making the strategy update step for a player i, the value of variable x would be the difference between the payoff $f^{\left(j\right)}$ of a neighboring player j, randomly selected by player i, and the payoff $f^{\left(i\right)}$ of i itself. Then, $\mathcal{F}(x>0;\beta )>0.5$ and $\mathcal{F}(x<0;\beta )<0.5$, making higher and lower than $0.5$, respectively, the probability for i to imitate j. The noise parameter $\beta >0$ allows to regulate the strength of the selection process. The selection pressure gets stronger and more deterministic when increasing $\beta$ (implying best response dynamics for $\beta \to \infty$), whereas it gets weaker and noisier when decreasing it, becoming completely random for $\beta \to 0$. In particular, for $\beta \ll 1$, $\mathcal{F}(x;\beta )\approx 1/2+\beta x/4$ and we get a linear dependence on x, i.e., a replicator-like rule, in weak selection regime.

The global state of a population of size N is defined by the density of cooperators:

$$c\left(t\right)=\frac{1}{N}\sum _{i=1}^N{\sigma }_i\left(t\right)$$

(5)

We initialize the system with a fraction $c_0\equiv c(t=0)$ of randomly selected cooperators.

From Equation (3) it immediately follows that $c=0$ and $c=1$ are absorbing states of the dynamics, for no change of state is possible from them.

3. Results

For every set of hypergraphs we report the effective asymptotic cooperators density $c_{\infty }$ and the convergence time t as a function of $\alpha$ and for different values of the conversion fraction p, as found in the performed Monte Carlo (MC) simulations; see Section 6.2 for further details. To systematically reveal the effect on the dynamics of converting 3-cliques in triangles, we report the results found for hypergraphs generated through the Holme–Kim model (HK) [34] and the Dorogotsev–Mendes model (DM) [35]. They both grow heterogeneous networks having topological properties typical of many real complex systems, and in particular social networks, like power-law degree distributions, the small-world property, and high transitivity (i.e., high clustering coefficient).

Let us call ${\alpha }_{\mathrm{surv}}\left(p\right)$ the value of $\alpha$ at which, for a given p, the asymptotic cooperator density $c_{\infty }$ starts to grow above zero, to be compared with the corresponding critical value ${\alpha }_{\mathrm{cr}}\left(p\right)$ for a well-mixed population. To this purpose, we define

$$r\left(p\right)\equiv \frac{{\alpha }_{\mathrm{surv}}\left(p\right)}{{\alpha }_{\mathrm{cr}}\left(p\right)}$$

(6)

where

$${\alpha }_{\mathrm{cr}}\left(p\right)=\sum _{m=m_{min}}^{m_{max}}m{p}_m$$

(7)

being $p_m$ the fraction of hyperedges of cardinality m in the hypergraph. Values of $r\left(p\right)$ smaller than 1 mean some level of structural reciprocity favoring the survival of cooperation. As a particular case, $r\left(0\right)$ measures the network reciprocity.

In Figure 2 we show the results for MC simulations with $c_0=0.5$ and two values of $\beta$, on hypergraphs stemming from HK networks with $N=500$ vertices, $m=3$ (average degree $\langle k\rangle \approx 6$ on the 2-section), and $P_t=1$. Parameter $P_t$, the probability of making a triad-formation step after a preferential-attachment step in the HK model, has been set to 1 to get the maximum possible transitivity, thus making more evident the effects of increasing the fraction p of converted 3-cliques. We find that, as we increase p from 0 to 1, $r\left(p\right)$ decreases monotonically from 0.78 to 0.60 when $\beta =1$, and from 0.79 to 0.63 when $\beta =0.1$. Thus, reciprocity is enhanced when we increase p, the fraction of 3-cliques promoted to triangular hyperedges. We also note how, by decreasing $\beta$ (weaker selection), the structural reciprocity decreases as well.

Looking at the convergence times in Figure 2, we see the typical peaked shape characterizing the critical region. Dealing with small finite systems, we expect the convergence time to show a peak near the transition, but still to remain finite. Since the maximum number of allowed time steps $t_{max}$ is fixed to $2\times {10}^4$, and we average over 100 realizations, we can state that, if the average convergence time is smaller than 200, then all the trials converged and the final populations are perfectly homogeneous; note that this is only a sufficient condition. In this case, an effective cooperation density $c_{\infty }$ other than 0 or 1 represents an average of homogeneous final states. This is the case for r values out of the transition region.

Counting how many times the system has not converged $n_{NC}$ in 100 realizations, and knowing the value of the average convergence time t, the average time $t_C$ the system took to converge in the remaining $100-n_{NC}$ trials is $t_C=100[t-n_{NC}(t_{max}/100)]/(100-n_{NC})$. For the previous hypergraphs we find $n_{NC}\le 12.3$ for $\beta =1.0$ and $n_{NC}\le 2.4$ for $\beta =0.1$, in the critical regions. The maximum values of $t_C$ are found to be around 1490 and 710, respectively. Being $t_C$ well below $t_{max}$, we expect that the non-convergence of some realizations does not depend on the value chosen for the latter. This suggests that the lack of convergence is not a purely dynamical effect—like at the critical point for an infinite well-mixed population—but it is due to the topological constraints of the structure. In fact, we find the presence of configurations of local states preventing the system to escape from them, acting as topological traps, see Section 4.3.

We also note that with the increasing of p, the critical region gets wider, whatever is the value of $\beta$. In Section 4.2.2, we give an explanation for this based on the degree heterogeneity and the local clustering of the structure.

Quite similar results are found for hypergraphs stemming from DM networks. In Figure 3 we show the results obtained for base networks with $N=500$ vertices ($\langle k\rangle =4$). Going from $p=0$ to $p=1$, we get the following monotonically decreasing values of $r\left(p\right)$: from 0.75 to 0.58 for $\beta =1$, from 0.78 to 0.62 for $\beta =0.1$. Furthermore, we find, respectively, $n_{NC}\le 39.7$ and $n_{NC}=0$ in the critical regions. The maximum values of $t_C$ are found to be around 2660 and 380. As for HK hypergraphs, setting $t_{max}=5\times {10}^4$, we find no significant differences. In Appendix A we report the results obtained also for $c_0=0.3$ and $c_0=0.7$.

In the case of hypergraphs derived from Erdős–Rényi networks (ER) [36], the transitions are always sharp and just below $r=1$ (see Appendix B), given the high structural homogeneity and low clustering of ER networks. The average convergence time is found not over 300 steps, along with the absence of traps.

4. Mathematical Analysis

4.1. Mean-Field Approximation

The mean-field (MF) approach relies on two approximations: (i) all the vertices have the same m-degree sequence $(\langle k_{m_{min}}\rangle ,\dots ,\langle k_{m_{max}}\rangle )$, averaged over the original hypergraph with m-degree sequence distribution $P(k_{m_{min}},\dots ,k_{m_{max}})$; (ii) the local states of different vertices are not correlated. Then, a simple master equation for the density of cooperators $c\left(t\right)$ follows:

$$c(t+1)=c\left(t\right)+c\left(t\right)\left(1-c\left(t\right)\right)\left(\mathcal{F}(\Delta f)-\mathcal{F}(-\Delta f)\right)$$

(8)

where

$$\Delta f=\langle f^{\left(c\right)}\rangle -\langle f^{\left(d\right)}\rangle =b\sum _{m=m_{min}}^{m_{max}}\left(\frac{\alpha }{m}-1\right)\langle k_m\rangle$$

(9)

Imposing $c(t+1)=c\left(t\right)$ we get, besides the absorbing state solutions $c\left(t\right)=0$ and $c\left(t\right)=1$, $\forall \alpha$, the stationary solution corresponding to $\mathcal{F}(\Delta f)=\mathcal{F}(-\Delta f)$, that is, $\Delta f=0$, given by

$$\alpha ={\alpha }_{cr}\equiv \frac{{\displaystyle \sum _{m=m_{min}}^{m_{max}}}\langle k_m\rangle }{{\displaystyle \sum _{m=m_{min}}^{m_{max}}{m}^{-1}\langle k_m\rangle }}=\sum _{m=m_{min}}^{m_{max}}m{p}_m=\langle m\rangle$$

(10)

for all $c\left(0\right)\in (0,1)$, being

$$p_m=\frac{m^{-1}\langle k_m\rangle }{{\displaystyle \sum _{m=m_{min}}^{m_{max}}{m}^{-1}\langle k_m\rangle }}$$

(11)

the fraction of hyperedges of cardinality m or, for a well-mixed population, the probability of playing in a group of size m, and $\langle m\rangle$ their average cardinality. Equation (10) gives the expected abrupt phase transition between full defection and full cooperation, for an infinite well-mixed population, in the presence of either one or more allowed group sizes (corresponding to uniform and non-uniform hypergraphs, respectively).

To be precise, the shown abrupt transition holds in the limit of an infinite population. The finite-size, trivially-structured version of a well-mixed population is an all-connected-to-all structure. Then, considering a m-uniform complete simple hypergraph of N vertices, it can be shown (see Appendix C) that, for $\beta b\gg 1$, the transition is at $\alpha \equiv {\alpha }_{cr}^m=m(N-1)/(N-m)$, greater than ${\alpha }_{cr}=m$ for any finite value of N. Additionally, ${\alpha }_{cr}^m/{\alpha }_{cr}$ grows with m. This finite-size effect slightly raises the structural reciprocity herein reported.

4.2. Invasion Analysis

The MF prediction is clearly reliable only for highly homogeneous structures. In all the other cases, instead, it serves as a reference point for measuring the structural reciprocity due to the complex heterogeneity of connections. In this section, we study the conditions for the invasion of the structured population, or for the resistance to be invaded for a player, using a certain strategy in specific configurations. We give a justification to the higher reciprocity level and the enlargement of the transition region, both found when increasing p in the used heterogeneous hypergraphs.

4.2.1. The Role of Cooperator Hubs

For all the analyzed scale-free networks in Figure 2 and Figure 3, we observe the transition moving towards notable smaller values of r, whatever is the value of p. This is not the case for homogeneous hypergraphs, like the ones derived from ER networks. We show here that the key difference is the presence of hubs in the former networks.

We focus on a cooperator hub with m-degree sequence $(k_{m_{min}}^{hub},\dots ,k_{m_{max}}^{hub})$, surrounded by neighbors with average k-degree sequence $({\overline{k}}_{m_{min}},\dots ,{\overline{k}}_{m_{min}})$, where

$${\overline{k}}_m={\overline{k}}_m\left(\left\{k_m^{hub}\right\}\right)=\sum _{k_m}P\left(k_m\right|k_m^{hub})k_m$$

(12)

Imposing the payoff of the hub playing as cooperator is equal to the average payoff of its neighbors playing as defectors, $f^{\left(hub\right)}={\overline{f}}^{\left(d\right)}$, we get the threshold value ${\alpha }_{\mathrm{th}}^{\mathrm{hub}}$ above which, whenever the hub and any of its defector neighbors match for the strategy update, the survival of the hub is favored. Indicating with ${{\overline{q}}_0}^{\left(d\right)}=c_{0}N/(N-1)$ and $q_0^{\left(hub\right)}={{\overline{q}}_0}^{\left(d\right)}-1/(N-1)$ the initial ($t=0$) average fraction of cooperating neighbors, respectively, of the hub’s defector neighbors and of the hub itself, from Equation (3) we find

$${\alpha }_{\mathrm{th}}^{\mathrm{hub}}=\frac{{\displaystyle \sum _{m=m_{min}}^{m_{max}}}{k}_m^{\mathrm{hub}}}{{\displaystyle \sum _{m=m_{min}}^{m_{max}}}\left[{\displaystyle \frac{k_m^{\mathrm{hub}}}{m}+\frac{m-1}{m}\left(q_0^{\left(hub\right)}{k}_m^{\mathrm{hub}}-{{\overline{q}}_0}^{\left(d\right)}{\overline{k}}_m\right)}\right]}$$

(13)

Considering s-uniform hypergraphs and putting ${\overline{k}}_s={\xi }_s{k}_s^{hub}$, from Equation (13) follows

$$r^{\mathrm{hub}}\equiv \frac{{\alpha }_{\mathrm{th}}^{\mathrm{hub}}}{{\alpha }_{cr}}=\frac{1}{1+(s-1)(q_0^{\left(hub\right)}-{{\overline{q}}_0}^{\left(d\right)}{\xi }_s)}$$

(14)

which is a monotone decreasing function of s and, for ${\xi }_s<1-{\left(N{c}_0\right)}^{-1}$, of $c_0$ as well (left panel of Figure 4). Taking as hubs those vertices with degree equal or larger, respectively, to the 97th and the 99th percentile of the degree distribution, we get, for each value of s, a small interval in which, the values of $r\left(p\right)$ found in the simulations for $c_0=0.5$, fit well. In fact, we expect Equation (13) to work well for intermediate values of $c_0$, as commented in Appendix D. Moreover, we expect the method to work for $\beta$ large enough, when the imitation process is more deterministic.

To evaluate Equation (13), we take $k_s^{\mathrm{hub}}$ as the average degree of those vertices classified as hubs and then compute the corresponding average value of ${\xi }_s$. Setting $c_0=0.5$, for HK hypergraphs we find $r^{\mathrm{hub}}$ between 0.71 and 0.79 for $s=2$, and between 0.55 and 0.65 for $s=3$; for DM hypergraphs we get $r^{\mathrm{hub}}$ between 0.72 and 0.79 for $s=2$, and between 0.56 and 0.65 for $s=3$. In Figure 4 (left panel), to visualize $r^{\mathrm{hub}}$ versus s, we reported the values taken by the former computed considering ${\xi }_s$ as independent of s and fixed equal to $({\xi }_2+{\xi }_3)/2$, taking ${\xi }_2$ and ${\xi }_3$ as the average of the respective values found in the used hypergraphs. In fact, only ${\xi }_2$ and ${\xi }_3$ are accessible in HK and DM hypergraphs. Nevertheless, ${\xi }_s$ only decreases by about $0.01$ going from $s=2$ to $s=3$, therefore keeping it fixed is a sensible choice.

When $\alpha >{\alpha }_{\mathrm{th}}^{\mathrm{hub}}$, a cooperator hub is resilient to the invasion by any defector in its neighborhood and, by persisting in the cooperative state, gives rise to a cooperative community whose members sustain each other (for a similar dynamics on networks see [17]). In this way, this community is able to enlarge itself further or, at least, to resist to be invaded by the defectors in the neighborhoods, thanks to the local topological constraints. Remarkably, in HK and DM networks, by construction, the hubs are frequently very close to each other (first or second neighbors), making it easier for the cooperative community to thrive when more than one of them happens to be in the cooperative state. When, instead, the hub is a defector, its local proliferation is favored whatever is the value of $\alpha$. However, once a defector hub easily succeeded in converting all its lower-degree cooperator neighbors in defectors, none of the players in this defective community provide/receive any benefit playing within it. This prevents the defective community to cover the entire population, unless $\alpha$ is enough small, outside the critical region.

4.2.2. The Role of Isolated Defectors

Let us now consider the limit case in which a ${\sigma }_1$-strategist is totally surrounded—up to its $z^{th}$ neighborhood, with z large enough—by ${\sigma }_2$-strategists, with ${\sigma }_1\ne {\sigma }_2$, to calculate the joint probability that the ${\sigma }_1$-strategist survives and some of the ${\sigma }_2$-strategists adopt strategy ${\sigma }_1$, as a function of $\alpha$. The probability $p_i\left(t\right)$ that a player of strategy ${\sigma }_i$ and neighbors $\mathsf{\Omega }\left(i\right)$ changes its strategy, reads:

$$p_i\left(t\right)=\frac{1}{\left|\mathsf{\Omega }\left(i\right)\right|}\sum _{\underset{{\sigma }_j\ne {\sigma }_i}{j\in \mathsf{\Omega }\left(i\right)}}\mathcal{F}\left(f^{\left(j\right)}\left(t\right)-f^{\left(i\right)}\left(t\right)\right)$$

(15)

In particular, considering a defector d of degree pair $(k_2^{\left(d\right)},k_3^{\left(d\right)})$ as the isolated strategist, the probability for it to survive as a defector at time t is equal to $1-p_d\left(t\right)$, in which the difference between the payoff $f^{\left(j\right)}\left(t\right)$ of a cooperator neighbor j and its payoff $f^{\left(d\right)}\left(t\right)$ is given by:

$$f^{\left(j\right)}\left(t\right)-f^{\left(d\right)}\left(t\right)=b\left[\alpha \left(k_2^{\left(j\right)}-\frac{1}{2}{k}_2^{\left(d\right)}+k_3^{\left(j\right)}-\frac{2}{3}{k}_3^{\left(d\right)}-a^{\left(j\right)}\right)-k_2^{\left(j\right)}-k_3^{\left(j\right)}\right]$$

(16)

where $a^{\left(j\right)}$ is equal to $1/2$ if j is a neighbor through a link and to ${\kappa }_d^{\left(j\right)}/3$ if it is a neighbor through a triangle, where ${\kappa }_d^{\left(j\right)}$ is the number of triangles that j shares with vertex d. Similarly, the probability that a cooperator neighbor j becomes a defector is equal to $\mathcal{F}\left(f^{\left(d\right)}\left(t\right)-f^{\left(j\right)}\left(t\right)\right)/\left|\mathsf{\Omega }\left(j\right)\right|$.

To get interesting insights, it is sufficient to consider all the defector’s neighbors (all cooperators) as equivalent to the average vertex in the neighborhood. Therefore, they have the same degree pair $(k_2^{\left(c\right)},k_3^{\left(c\right)})$, the same number of neighbors $\left|\mathsf{\Omega }\right|$, the same ${\kappa }_d$, and, consequently, the same payoff $f^{\left(c\right)}\left(t\right)$; it makes sense because the neighbors are selected uniformly at random when updating the strategies. Then, the probability for d of surviving as a defector and invading all at once $w\le \left|\mathsf{\Omega }\left(d\right)\right|$ neighboring sites reads

$$(1-p_d\left(t\right))·{\left[\frac{1}{\left|\mathsf{\Omega }\right|}\mathcal{F}\left(f^{\left(d\right)}\left(t\right)-f^{\left(c\right)}\left(t\right)\right)\right]}^w$$

(17)

In Figure 4 (right panel) we show this combined probability computed for $w=1$, $\left|\mathsf{\Omega }\left(d\right)\right|=8$, different values of $\left|\mathsf{\Omega }\right|$, ${\kappa }_d=1.75$, and in the presence of either only links ($p=0$) or only triangles ($p=1$). These values correspond to neighborhoods compatible with DM and HK hypergraphs (see Appendix E for details). Interestingly, the defector should prefer interacting through triangles whenever $\left|\mathsf{\Omega }\left(d\right)\right|>\left|\mathsf{\Omega }\right|$ and through links otherwise. This is simply understood observing that, thanks to the percolation of 3-cliques (or triangles) present in HK and DM networks (hypergraphs), ${\left.k_3^{\left(d\right)}\right|}_{p=1}$ is close to ${\left.k_2^{\left(d\right)}\right|}_{p=0}$; in particular, ${\left.k_3^{\left(i\right)}\right|}_{p=1}={\left.k_2^{\left(i\right)}\right|}_{p=0}-1$, $\forall i$, for DM structures. Since when d plays within a triangle obtains two times what it gets playing within a link, it is advantageous playing through triangles whenever it possesses a sufficient number of neighbors. In Figure 4 we also note that, for $\left|\mathsf{\Omega }\right|$ sufficiently smaller than $\left|\mathsf{\Omega }\left(d\right)\right|$, the probability of proliferating for the defective strategy counter-intuitively slightly increases with $\alpha$ (or r), even above the mean-field critical point ($r=1$).

Taken together, these observations suggest that such dependence on $\alpha$ may contribute substantially to the reported enlargement of the critical region when p is increased. Specifically, on one hand, the increasing of $\alpha$ always brings the cooperators’ performance closer to the defectors’ in any local instance of the game, making it easier for cooperators to invade (with certainty, eventually); on the other hand, defectors with high degrees are increasingly able to invade their neighborhoods, containing the diffusion of the cooperative strategy. This competitive dynamics, existing only for heterogeneous structures, is thus expected to delay the point at which cooperators succeed in fully invading the population, stretching out the transition, which is almost abrupt in homogeneous random hypergraphs.

Next, we will explain why this competition gets fiercer when p is increased, enlarging even more the transition region. From Equation (16) we can calculate the condition that the degrees of the defector and of its cooperator neighbors have to fulfill in order that $f^{\left(c\right)}-f^{\left(d\right)}$ decreases when $\alpha$ increases, thus raising the competition. Taking, as done before, all the defector’s neighbors (all cooperators) as equivalent, by imposing $\partial (f^{\left(c\right)}-f^{\left(d\right)})/\partial \alpha <0$, we get

$$\frac{\left|\mathsf{\Omega }\left(d\right)\right|+1}{\left|\mathsf{\Omega }\left(d\right)\right|}\left(\frac{1}{2}{k}_2^{\left(d\right)}+\frac{2}{3}{k}_3^{\left(d\right)}\right)>k_2^{\left(c\right)}+k_3^{\left(c\right)}$$

(18)

where $\left|\mathsf{\Omega }\left(d\right)\right|=k_2^{\left(d\right)}+2k_3^{\left(d\right)}/{\kappa }_d$ is the number of neighbors of the isolated defector, and ${\kappa }_d$ is the number of triangles that each of the cooperator neighbors shares with it. Since $\left(\right|\mathsf{\Omega }\left(d\right)|+1)/\left|\mathsf{\Omega }\right(d\left)\right|\le 3/2$, Equation (18) requires some overall degree heterogeneity in favor of the isolated defector. Instead, as anticipated, for a regular hypergraph it always holds $\partial (f^{\left(c\right)}-f^{\left(d\right)})/\partial \alpha >0$. Thus, when in a heterogeneous structure a defector possesses enough higher degrees with respect to the surrounding cooperators, its survival is favored when $\alpha$ is increased, smearing the transition.

To get some understanding about how this phenomenon depends on p, i.e., on the order of the interactions, we refer to hypergraphs that become 3-uniform when $p=1$, and 2-uniform when $p=0$, like the ones we have used. Then, Equation (18) splits into the following conditions for $p=0$ and $p=1$, respectively:

$$\begin{array}{cc}\hfill \frac{k_2^{\left(d\right)}}{k_2^{\left(c\right)}}& >2-\frac{1}{k_2^{\left(c\right)}}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{k_3^{\left(d\right)}}{k_3^{\left(c\right)}}& >\frac{3}{2}-\frac{{\kappa }_d}{2k_3^{\left(c\right)}}\hfill \end{array}$$

(20)

Equating the expressions taken by $\left|\mathsf{\Omega }\right(\sigma \left)\right|$ for $p=0$ and $p=1$, the identity ${\left.k_3^{\left(\sigma \right)}\right|}_{p=1}={\left.k_2^{\left(\sigma \right)}\right|}_{p=0}·{\kappa }_{\sigma }/2$, to be read as $k_3^{\left(\sigma \right)}$ computed when $p=1$, given $k_2^{\left(\sigma \right)}$ when $p=0$, follows. We can then rewrite Equation (20) in terms of $k_2^{\left(d\right)}$ and $k_2^{\left(c\right)}$ and directly compare it with Equation (19). Equation (20) becomes

$$\frac{k_2^{\left(d\right)}}{k_2^{\left(c\right)}}>\frac{3}{2}\frac{{\kappa }_c}{{\kappa }_d}-\frac{1}{k_2^{\left(c\right)}}$$

(21)

where the degrees are calculated for $p=0$. It follows that Equation (21) is a weaker condition than Equation (19), i.e., Equation (19) ⇒ Equation (21), if and only if ${\kappa }_d/{\kappa }_c>3/4$. Remarkably, whenever the latter inequality is fulfilled, $\partial (f^{\left(c\right)}-f^{\left(d\right)})/\partial \alpha <0$ is more negative for $p=1$ than for $p=0$. Expressing ${\kappa }_{\sigma }$ in terms of the local clustering coefficient $C\left(\sigma \right)$ and the link-degree $k_2^{\left(\sigma \right)}$ as ${\kappa }_{\sigma }={[C\left(\sigma \right)·(k_2^{\left(\sigma \right)}-1)]}_{p=0}$ (see Appendix E for the derivation), we finally get

$$\mathrm{Equation}\left(19\right)\Rightarrow \mathrm{Equation}\left(21\right)⟺{\left[\frac{C\left(d\right)·(k_2^{\left(d\right)}-1)}{C\left(c\right)·(k_2^{\left(c\right)}-1)}>\frac{3}{4}\right]}_{p=0}$$

(22)

According to the inequalities in Equations (19) and (21), the most interesting case is obtained when $k_2^{\left(d\right)}>k_2^{\left(c\right)}$. Then, Equation (22) can be satisfied even when the isolated defector possesses a smaller clustering coefficient than its cooperators neighbors. Specifically, HK and DM models grow networks whose local clustering coefficient C scales with the 2-degree $k_2$ as $C\left(k_2\right)\sim k_2^{-\gamma }$, where $\gamma$ is around 0.8 for HK ($P_t=1$) and exactly 1 for DM (see Appendix E). Substituting in Equation (22), one finds that the above implication holds whenever $k_2^{\left(d\right)}>k_2^{\left(c\right)}$, given $\gamma \le 1$.

If we now ask for $k_2^{\left(d\right)}$ such that, for a given $k_2^{\left(c\right)}$, Equation (21) is satisfied but Equation (19) is not, we find the condition

$$\frac{3}{2}\frac{{\kappa }_c}{{\kappa }_d}{k}_2^{\left(c\right)}-1<k_2^{\left(d\right)}\le 2k_2^{\left(c\right)}-1$$

(23)

where the degrees are calculated for $p=0$. Therefore, if $k_2^{\left(d\right)}$ takes an intermediate value in the sense of Equation (23), with respect to $k_2^{\left(c\right)}$, Equation (21) is satisfied but Equation (19) is not, allowing the competitive dynamics to take place for $p=1$ but not for $p=0$. As shown in Appendix E, such intermediate neighborhoods actually exist in the structures we made use of. It must be noted that, for $k_2^{\left(d\right)}>2k_2^{\left(c\right)}-1$, the competitive dynamics exists for $p=0$ as well. However, as already observed, due to the linear dependence of $f^{\left(c\right)}-f^{\left(d\right)}$ on $\alpha$, $\partial (f^{\left(c\right)}-f^{\left(d\right)})/\partial \alpha <0$ is more negative for $p=1$ than for $p=0$ whenever $k_2^{\left(d\right)}>k_2^{\left(c\right)}$ (corresponding to ${\kappa }_d/{\kappa }_c>3/4$), making the competitive dynamics always fiercer in the former case.

Ultimately, this analysis shows that in heterogeneous and clustered structures, like the ones we used here, the competitive dynamics takes place more frequently and intensively in the presence of 3-way ($p>0$) than 2-way ($p=0$) interactions, in accordance with the observed enlargement of the transition region when increasing p. It must be said that our analysis does not exclude the existence of other mechanisms that could further contribute to the observed enlargement.

4.3. Heterogeneous Populations: The Role of Topological Traps

When running a dynamical process on top of a connected structure, there could be microscopic and mesoscopic configurations which operate as topological traps, blocking the system in a certain state. Traps represent here the only mechanism able to give rise to asymptotic heterogeneous populations in which cooperators and defectors co-exist. A trap can have any size and complexity, i.e., it can involve any number of vertices and develop over many time steps. A particularly simple class of traps is the bottleneck-type’s, consisting of a sub-hypergraph in which some vertices act as bridges between groups of vertices of opposite strategy. In rank-3 hypergraphs, they consist of some bridge-vertices gluing links and/or triangles, as shown in Figure 5.

A trap comes into play when the strategies are updated. According to our setup, it is only the 2-section of the hypergraph that matters for the strategy update, and the presence of higher-order interactions solely affects the ranges of values of $\alpha$ for which a local structure can act as a trap. Then, the effect of a bottleneck trap is fully specified by the comparison between the chosen value of $\alpha$ and the threshold values ${\alpha }_{\mathrm{th}}$s, defined as those values at which a cooperator and a defector, each with its m-degree sequence, get the same final payoff. Any trap other than bottleneck-type is specified simultaneously by a set of coupled ${\alpha }_{\mathrm{th}}$s, one for any link of the 2-section contributing to the trap.

Whatever is the complexity of a trap, its effectiveness is strictly related to the value of the parameter $\beta$. A trap is able to act if the selection is strong enough, i.e., for high enough values of $\beta$. Indeed, $\beta$ weighs the proportional contribution of the difference of the payoffs to the update rule. The smaller is $\beta$, the more random the strategy update is, being increasingly disentangled from the payoffs resulting from the interactions. Accordingly, the MC simulations made for $\beta =0.1$ provide convergence times notably shorter than the ones found for $\beta =1$. Furthermore, we see also a slight reduction of structural reciprocity when lowering $\beta$, meaning that the structures develop some traps able to further anticipate the survival of cooperation towards slightly lower values of r.

A critical value of $\alpha$ can either satisfy or not some of the local conditions for producing a trap. In other words, the found effective densities of cooperators can be the result of an average of both homogeneous final states (the evolutionary stable states) and/or heterogeneous ones (produced by the traps). It was already unveiled in [37] that, for the entire set of 2 × 2 social dilemmas played on different real networks, the presence of bottleneck traps is strictly related to the degree heterogeneity and the lack of redundant paths.

To characterize the generic traps for the PGG played on rank-3 simple hypergraphs, similarly to what is done in [37], we perform two different randomization procedures of the networks (see Section 6 for details). One randomization preserves both the degree and the local clustering coefficient of each vertex, while the other preserves the degree correlations up to the second order, i.e., up to the joint $3K$-distribution $P(k,k^{\prime },k^{″})$. Clustering coefficient and second-order degree correlations are clearly related, but preserving one does not fully imply preserving the other automatically. Now, since the action of a trap is defined in multiple time steps (two at least) in which are involved several vertices with overlapping neighborhoods, the degree correlations of order higher than two should be important in the trapping process. Destroying those correlations we expect the traps to dissolve. Quite surprisingly, we find that what really matters for the presence of traps other than the bottleneck-type, is the high value of the local clustering coefficient (see Figure 6). Preserving the $P(k,k^{\prime },k^{″})$ distribution without caring of the local clustering, unravels the initially present traps. We find this for both DM and HK hypergraphs. On the contrary, the presence of bottleneck traps is clearly negatively correlated with clustering; they are found to be dominant for lowly clustered networks (scarce of redundant paths) exhibiting strong degree heterogeneity, in line with [37].

5. Conclusions

Collective and, in particular, cooperative phenomena, in which many individuals take part as a group, are observed in many real systems, among which human society stands out. In this work we investigated systematically the effect that group (higher-order) interactions have on the evolution of cooperation in structured populations. We studied the dynamics arising from a population of individuals interacting by playing a public goods game on top of hypergraphs. While the theory we developed is general, the numerical simulations specifically regarded rank-3 hypergraphs. In this way, we distinguish whether an element, when a member of a group (clique), interacts separately, in pair, with each of the other members, or interacts simultaneously with all of them.

To get the hypergraphs, we made use of an algorithm preserving the pairwise projection of the structure, which is crucial to satisfy the constraints implicit in the network topology. In fact, although a pairwise projection can correctly describe who interacts with whom, it cannot account for different types of interactions, the temporal order in which the latter occur, and whether interactions involve pairs or larger groups of units. Additionally, the 2-section imposes the maximum possible rank one could reach by generalizing the herein given procedure. The existence of a maximal size for groups is usually related to a saturation restriction due to some unbalanced cost in forming larger groups. This generally depends on the nature of the considered system, and addressing its form goes beyond the scope of this work. Here, we just recognize the existence of such a restriction as encoded in the 2-sections of the used hypergraphs, and preserve it with our generative method.

Monte Carlo simulations clearly indicated that heterogeneous structures, like the ones obtained from HK and DM networks, are able to sustain cooperation for notable ranges of values of the synergy parameter below its critical value. This holds for unstructured populations, whatever the conversion fraction, related to the order of the interactions. Our findings, therefore, extend the phenomenon of network reciprocity to rank-3 hypergraphs.

Remarkably, we show a clear evidence that the sustenance of cooperation is stronger when the fraction of 3-way interactions is increased. Indeed, we attained a monotonous decrease of structural reciprocity with the increasing of conversion fraction, for hypergraphs derived from both HK and DM networks. Moreover, the level of cooperation is always larger for the hypergraph. It is key noting that this enhanced reciprocity is exclusively due to the substitution of some closed triads of first-order interactions (3-cliques) with unique second-order interactions (triangles). In other words, given a heterogeneously-structured population, cooperation thrives more easily if individuals interacting separately, in pairs, within a closed triad are made to interact all together. It should be noticed that, considering the synergy factor as independent of the group sizes, we did not account for explicit synergistic or anti-synergistic effects that could increase or decrease the found levels of reciprocity. The improvement found here, rather, can be referred to as a topological synergistic effect allowed by the evolutionary dynamics.

As could be expected from what is known on networks, an important requirement for a significant reciprocity is the heterogeneity of connections; a separate analysis, not addressed in this work, would be needed for higher-order generalizations of spatial networks. To prove this, we performed an invasion analysis able to provide a good estimate for the critical point of emergence of cooperation, for both hypergraphs with rank equal to 2 and to 3. More generally, the analysis indicated that, the higher the order of the interactions, the stronger the structural reciprocity is. To be more precise, the marked enhancement of reciprocity has been found for heterogeneous structures with a null or slightly negative assortativity, in line with what is reported for networks. According to our analysis, we expect a weaker beneficial effect on hypergraphs derived from assortative networks.

Furthermore, the invasion analysis allowed us to justify, in comparison with the almost abrupt transition found in homogeneous structures, the enlargement of the transition region observed for the used heterogeneous hypergraphs. Relying on their locally-clustered structure, it also allowed to partially explain why the enlargement grows with the increasing of conversion fraction.

Regarding the heterogeneity of strategies in the asymptotic population, we found the presence of topological traps preventing the system to converge to a uniform state. We characterized those traps by means of two randomization procedures. In the used clustered structures, we attended the break-up of the present traps whenever the chosen randomization did not preserve the local clustering, though still preserving second-order degree correlations. We show that traps represent the only way to get a non-uniform asymptotic population.

Finally, by tuning the selection pressure, we found an improvement towards cooperation in conditions of stronger selection. This is precisely what we expect from how topological traps work. Remarkably, this shows up the contribution coming from some traps to enlarge the region where cooperative behaviors can survive.

6. Methods

6.1. Generating Rank-3 Simple Hypergraphs

The algorithm to generate rank-3 simple hypergraphs proceeds as follows:

• A simple network $g=\left(V_0,E_0\right)$ is generated through some model;
• From the set of all the 3-cliques in it, a fraction $0\le p\le 1$ of them is picked uniformly at random;
• To each of the picked 3-cliques a triangle is associated: if $e_{12}=\left\{v_1,v_2\right\},e_{23}=\left\{v_2,v_3\right\},e_{13}=\left\{v_1,v_3\right\}\subset E_0$ are the three edges forming a chosen 3-clique over the subset $\left\{v_1,v_2,v_3\right\}\subset V_0$ of vertices, then the hyperedge (triangle) $e_{123}=\left\{v_1,v_2,v_3\right\}$ is added to $E_0$;
• To obtain a simple hypergraph, after a fraction p of 3-cliques has been converted to triangles, an edge is removed if it is subset of at least one triangle: the three edges $e_{12},e_{23},e_{13}$, being subsets of $e_{123}$, are removed from $E_0$.

This procedure generates a rank-3 simple hypergraph $H=\left(V,E\right)$, where $V=V_0$ and E is the new set of hyperedges constructed from $E_0$ through the steps 2–4. The initial network g is the 2-section of H. See Figure 7 for an illustration of the procedure.

Since the algorithm relies on the amount of 3-cliques present in the base network g, the more clustered is the latter, the stronger will be the effect of increasing the value of p. In HK ($P_t=1$) and DM networks, every link is part of at least a 3-clique. Therefore, by tuning p, they allow us to explore the entire range of hypergraphs from the 2-uniform ($p=0$) to the 3-uniform ($p=1$).

As explained in Appendix E, the above algorithm can be easily extended to get hypergraphs of any rank higher than 2.

6.2. Monte Carlo Simulations

Every simulation starts by picking, uniformly at random, a fraction $c_0\equiv c(t=0)$ of cooperators. For each fixed set of parameters $\beta$, $c_0$, and p, we run 100 simulations over a generated hypergraph and take the average of both the asymptotic cooperators density and the convergence time, i.e., the number of time steps the process took to stop. We do it for 10 realizations of the same kind of hypergraph (i.e., stemming from the same network model), and then average again over them. Whenever the maximum number of time steps (set to $2\times {10}^4$) is reached without converging to one of the absorbing states $c=0$ or $c=1$, the last value of c is considered.

6.3. Randomization Procedures

We briefly summarize the two randomization procedures implemented for the characterization of topological traps. The selected number of randomizing steps was ${10}^5$ for each network.

6.3.1. Preserving Local Clustering Coefficient

Through the following procedure one can randomize a network preserving, besides the degree, also the local clustering coefficient of each vertex. At each step, two vertices are uniformly chosen at random and their local clustering coefficient is computed. Two links, one for each vertex, are randomly selected and interchanged. If the clustering coefficients of the two vertices remain unchanged, the randomizing step is accepted; otherwise, another two vertices are drawn and the randomization is tried again.

6.3.2. Preserving Second-Order Degree Correlations

This procedure allows to obtain a randomized network with exactly the same degree correlations up to the second-order, i.e., preserving the joint $3K$-distribution $P(k,k^{\prime },k^{″})$. At each step, two vertices are uniformly chosen at random and the pair is accepted if only the vertices have equal degree. Then, two links, one for each vertex, are peaked at random (an interchange performed now would ensure the conservation of the $2K$-distribution $P(k,k^{\prime })$). The degrees of the two vertices at the end of the two selected links are compared and, if they match, the interchange is performed (now conserving also $P(k,k^{\prime },k^{″})$); else ways, two new initial vertices are chosen and the randomization is tried again.