Introduction

In the last two decades network science has largely contributed to understanding how collective behaviors emerge in complex systems. Representing and characterizing the intricate pattern of interactions among the constituents of a complex system as a graph1,2 has allowed to investigate how the system’s structure affects its dynamics3. A large variety of dynamical processes, ranging from percolation4 to synchronization5, epidemics6 and social cooperation7, has been considered and explored. This burst of activity in the so-called structure-function problem has led to the discovery of novel microscopic mechanisms that trigger and control collective behaviors.

One of the most significant phenomena studied in this line of research is the emergence of explosive transitions. Here, explosive refers to abrupt phase transitions resulting from microscopic dynamical mechanisms that depend on the underlying topology in such a way to hinder the formation of a macroscopic state8,9. Initially introduced for sharp (yet continuous) percolation processes10,11,12 where, for instance, the probability of linking two nodes depends on the size of the components present in the systems, explosive transitions have since been observed in synchronization dynamics13,14,15, e.g., in networks of phase oscillators in which the natural frequency of a node increases with its degree, and contagion dynamics with transmission probability between an infected and a susceptible depends on the states of the surrounding nodes16,17,18. These studies have unveiled different scenarios of abrupt transitions, covering discontinuous behavior with13,17 or without9,19,20 bistability (hysteresis), but always relying on the synergistic coupling between network structure and dynamics.

Recently, explosive transitions have attracted increasing interest in the context of higher-order networks21,22,23,24, a framework going beyond (pairwise) networks to capture how interactions in groups of more than two units influence collective behaviors25,26,27,28,29,30,31,32,33. Within this framework, first-order transitions with an associated bistable region are often observed, with examples spanning social contagion25,26,27, synchronization28,29,30, and game theory33,34,35. Battiston et al. highlighted this ubiquity in their perspective paper, asserting that higher-order interactions provide a general pathway to explosive transitions36. Moreover, Kuehn and Bick provided a general mathematical framework to study the conditions under which the nature of the transition changes from first-order to second-order, or viceversa37. In particular, they showed that the emergence of first order transition, when classical models are generalized by incorporating additional effects, such as non-pairwise interactions, is not surprising but rather a universally expected effect.

In the early stages, studies of dynamical processes on higher-order networks focused on understanding how the nature of the transitions changes by tuning the number and strength of higher-order interactions. Firstly, models of social contagion in homogeneous systems with uncorrelated higher-order interactions were considered. It was found that discontinuous transitions occur exclusively when the strength of higher-order interactions exceeds a critical threshold25,27,28. However, subsequent research38,39,40 revealed that the strength of the interactions is not the sole determinant of behavioral changes in social systems. Remarkable examples include the dynamical influence of heterogeneity in the node participation within groups (analyzed through group size distributions)40, and the role of giant higher-order components (HOC)41 in the stability of social contagion processes, even in absence of pairwise interactions42.

Higher-order networks have usually been modeled either as hypergraphs with uncorrelated hyperedges, or as simplicial complexes. These models are at the two opposite ends of a continuous spectrum of possible structures, as a simplicial complex can be seen as a hypergraph with the constraint, known as inclusion (or downward closure), that, when an interaction exists between a set of nodes, all interactions involving all possible subsets of nodes also exist. For example, if in a simplicial complex we have a simplex of dimension 2, that is the three nodes are involved in a three-body interaction, then, the simplicial complex also contains the three underlying pairwise interactions (links). However, in empirical data the inclusion property is satisfied for some of the three-body interactions but not for all of them43,44. Various works have explored how higher-order interactions shape collective dynamics differently in hypergraphs and simplicial complexes45,46,47,48,48,49. In particular, Burgio et al. have studied how the dynamics of a contagion process depends on the extent to which interactions satisfy the downward closure. They found that when two-body interactions are highly nested in three-body interactions, the epidemic threshold is reduced, but smaller outbreaks occur.

To better understand the intricate interplay between structure and dynamics, in this work, we focus on another important, but still unexplored, aspect of the microscopic organization of the nodes in groups. Namely, we will explore how the dynamics of a system is influenced by the way in which group interactions of the same order are correlated, e.g., in the particular case of triads, by how the three-body interactions in a network overlap between each others. To quantify such correlations, we introduce a novel metric, that we call the intra-order hyperedge overlap of order m, measuring the number of nodes shared among hyperedges of the same order m. Equipped with this metric, our work aims to address three main questions: (i) whether overlaps between interactions of the same order are indeed present in real complex systems, (ii) how different levels of intra-order hyperedge overlap drive the emergence of explosive transitions, and (iii) to what extent the effects of overlap are general, in the sense that they apply to different types of dynamical processes, e.g., synchronization and contagion, which are two of the most studied dynamics in complex systems.

Our results section is organized as follows. In the first subsection, we introduce our metric and we use it to explore real-world systems from various domains, including human face-to-face interactions and brain connectomes, which will be particular relevant for the dynamical processes considered in the following sections. We find that hypergraphs describing real-world systems exhibit a wide range of values of intra-order hyperedge overlap. For this reason, in the second and third subsections we investigate the effects of a tunable level of intra-order hyperedge overlap on the emergence and properties of the collective behavior of a system. In the second subsection, we focus on one of the most studied social dynamics involving groups, namely complex contagion. Here, we show that a low intra-order hyperedge overlap leads to explosive transitions characterized by a bistable region, where both an active and an absorbent state coexist. Conversely, hypergraphs with an intra-order hyperedge overlap larger than a critical value can only exhibit continuous transitions. To show that these results extend beyond complex spreading to other dynamics, in the third subsection we study synchronization of coupled dynamical systems with higher-order interactions. Synchronization is a collective behavior that emerges in various social and biological complex systems. E.g. it is relevant to explain brain functioning, both in healthy (e.g., vision50, movement51, and memory52,53) and in pathological (e.g., epilepsy54,55,56,57,58,59,60) conditions. We consider here an extension, with higher-order interactions, of the Kuramoto model61, a system of coupled oscillators that has also been used to study synchronization of the brain at the cortical level62,63,64. Similarly to the case of complex contagion, we also find that the nature of the transition from an incoherent to a synchronized state in the Kuramoto model can be tuned by the value of the intra-order hyperedge overlap. Namely, we switch from explosive to smooth transitions as the intra-order hyperedge overlap of the underlying hypergraph increases. Overall, our results show that, not only the presence of higher-order interactions, but also the intricate way in which they coordinate, shapes the emergence of collective behaviors in complex systems.

Results

Quantifying intra-order correlations in higher-order networks

We model a system with higher-order interactions as a hypergraph \(H=({{{\mathcal{N}}}},{{{\mathcal{E}}}})\), where \({{{\mathcal{N}}}}\) is the set of \(N=| {{{\mathcal{N}}}}|\) nodes and \({{{\mathcal{E}}}}\) is the set of hyperedges, describing the \(E=| {{{\mathcal{E}}}}|\) interactions in groups of two or more nodes. Each hyperedge \(e\in {{{\mathcal{E}}}}\) is a subset of the nodes in \({{{\mathcal{N}}}}\) and is characterized by its order, m, which is defined in terms of its cardinality, ∣e∣, as m = ∣e∣ − 1. So hyperedges of order 1 represent pairwise interactions, interactions in group of 3 nodes correspond to hyperdges of order 2, and so on. For each node i, \({{{{\mathcal{E}}}}}_{i}^{(m)}\) is the set of hyperedges ej of order m such that i ∈ ej, and \({k}_{i}^{(m)}\) is its generalized degree of order m, defined as the cardinality of the subset \({{{{\mathcal{E}}}}}_{i}^{(m)}\), i.e., \({k}_{i}^{(m)}=| {{{{\mathcal{E}}}}}_{i}^{(m)}|\)65. For each order m, we can also define an adjacency tensor A(m) whose generic element \({a}_{{i}_{0},\ldots,{i}_{m}}^{(m)}=1\) if the m-hyperedge containing nodes i0, …, im exists, and zero otherwise. From the definition of the adjacency tensor of order m, we can define the generalized density of m-hyperedges, δ(m) as the ratio between the number of 2-hyperedges that are present in the hypergraph and the number of possible 2-hyperedges:

$${\delta }^{(m)}=\frac{\frac{1}{(m+1)!}{\sum }_{{i}_{0},{i}_{1},\ldots,{i}_{m}}{a}_{{i}_{0},{i}_{1},\ldots,{i}_{m}}^{(m)}}{\left(\begin{array}{c}N\\ m+1\end{array}\right)}\ .$$
(1)

The hyperedges to which node i belongs determine which nodes interact or not with i, and thus which nodes influence or not its dynamics. In the case of a simple graph, each link of node i connects it to a distinct neighbor. Conversely, in a hypergraph, a node i with two or more hyperedges of order m ≥ 2 can share one or more neighboring nodes with different hyperedges. Figure 1a–c shows three different configurations for a node with four 2-hyperedges. The extent of overlap of the hyperedges of i varies, resulting in diverse microscopic structures ranging from the non-overlapping case in Fig. 1a to the extreme case of maximal intra-order hyperedge overlap in Fig. 1c. To measure the local overlap among hyperedges of order m we introduce the node intra-order hyperedge overlap \({T}_{i}^{(m)}\) of a node i, defined for \({k}_{i}^{(m)} > 1\) as:

$${T}_{i}^{(m)}=1-\frac{{{{{\mathcal{S}}}}}_{i}^{(m)}-{{{{\mathcal{S}}}}}_{i}^{(m),-}}{{{{{\mathcal{S}}}}}_{i}^{(m),+}-{{{{\mathcal{S}}}}}_{i}^{(m),-}},$$
(2)

while, for \({k}_{i}^{(m)}\le 1\), we set \({T}_{i}^{(m)}=0\) since no neighbors can overlap. In Eq. (2) \({{{{\mathcal{S}}}}}_{i}^{(m)}\) is the number of unique neighbors of node i in interactions of order m, whereas \({{{{\mathcal{S}}}}}_{i}^{(m),-}\) (\({{{{\mathcal{S}}}}}_{i}^{(m),+}\)) accounts for the minimum (maximum) number of unique neighbors, in interactions of order m, that i can have. Note that the expressions for the quantities \({{{{\mathcal{S}}}}}_{i}^{(m),+}\) and \({{{{\mathcal{S}}}}}_{i}^{(m),-}\) depend on the node generalized degree \({k}_{i}^{(m)}\), and are discussed in Methods. The node intra-order hyperedge overlap in Eq. (2) takes its minimum, i.e., \({T}_{i}^{(m)}=0\), when none of the neighbors is shared with other hyperdeges of order m, while it takes its maximum, i.e., \({T}_{i}^{(m)}=1\), when each neighbor is shared by two or more hyperedges of order m. For the three configurations in Fig. 1a–c, where m = 2, Eq. (2) gives \({T}_{i}^{(2)}=1\), \({T}_{i}^{(2)}=0.5\), and \({T}_{i}^{(2)}=0\), respectively. Finally, to characterize the degree of intra-order hyperedge overlap at the global scale of the whole hypergraph we average \({T}_{i}^{(m)}\) over all the nodes:

$${{{{\mathcal{T}}}}}^{(m)}=\frac{{\sum }_{i}{k}_{i}^{(m)}{T}_{i}^{(m)}}{{\sum }_{i}{k}_{i}^{(m)}}.$$
(3)

Note that, due to its normalization, also the quantity \({{{{\mathcal{T}}}}}^{(m)}\) takes values in the range [0, 1].

Fig. 1: Intra-order hyperedge overlap in synthetic and real-world systems.
figure 1

a–c Examples of three configurations with different levels of node intra-order hyperedge overlap for a node i with generalized 2-degree \({k}_{i}^{(2)}=4\): a \({T}_{i}^{(2)}=1\); b \({T}_{i}^{(2)}=0.5\); c \({T}_{i}^{(2)}=0\). d Intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\) and generalized density δ(2) of real-world hypernetworks. The dashed vertical line corresponds to the density of the synthetic hypernetworks with tunable \({{{{\mathcal{T}}}}}^{(2)}\) that we have used to study the effects of overlap on dynamical processes. Diamonds, circles and square refer respectively to hypernetworks constructed from face-to-face contacts (CD), brain connectomes (BC), and other systems (RH). Information on each hypernetwork can be found through its label’s code in Supplementary Table I–III. The inset shows the relation between the intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\) and the average generalized degree 〈k(2)〉.

Full size image

To show that real-world systems can exhibit different levels of intra-order hyperedge overlap, we have analyzed 39 hypernetworks describing interactions in social and biological systems. Of these, 39 hypernetworks 7 have been constructed from high-resolution face-to-face contact data collected in various social contexts66,67,68,69,70,71,72, whereas 12 have been built upon biological data from the brain connectomes of different animals73,74,75,76,77,78,79 (more details on the datasets and how they are processed are described in the Methods section). The remaining 20 hypernetworks describe higher-order relations in different systems80,81,82,83,84,85,86,87, ranging from coauthorship in scientific papers80,81 to relations derived from users’ actions on very well-known websites83,85,87 (more details on these datasets can also be found in the Methods section). Figure 1d reports the average generalized density δ(2) and the intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\) for each of these hypernetworks. The results show that the values of \({{{{\mathcal{T}}}}}^{(2)}\) span the whole range [0, 1]. In particular, we find that hypernetworks with the same low value of δ(2) can display very different values of the intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\). This is an indication of the independence of the two metrics δ(2) and \({{{{\mathcal{T}}}}}^{(2)}\) in sparse networks. E.g., for δ(2) = 0.2 × 10−4 (see vertical dashed line), we observe that the intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\) of real-world systems can vary from values close to 0 to values as large as 0.8. Similar conclusions can be drawn from the inset in Fig. 1d, where we report the average generalized degree \(\langle {k}^{(2)}\rangle={\sum }_{i}{k}_{i}^{(2)}/N\) and the intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\) for each of the hypernetworks.

We will now show that this diversity of \({{{{\mathcal{T}}}}}^{(2)}\) plays a crucial role in shaping the overall collective behavior of a system. To investigate in a systematic way the effect of the intra-order hyperedge overlap on dynamical processes, we introduce a method to construct hypergraphs with a tunable value of \({{{{\mathcal{T}}}}}^{(2)}\), given a fixed regular degree distribution. We start from a configuration with maximum intra-order hyperedge overlap, i.e., \({{{{\mathcal{T}}}}}^{(2)}=1\), obtained by letting each node to share the same number of neighbors at each order of hyperedges (in particular, 1-hyperedges and 2-hyperedges). Then, the value of \({{{{\mathcal{T}}}}}^{(2)}\) is tuned by rewiring a fraction of the existing 2 − hyperedges, without changing neither the nodes generalized degree \({k}_{i}^{(2)}\), nor the average generalized degree 〈k(2)〉, nor the generalized density δ(2) (more details on the construction are provided in the Methods section). In this way, we obtain hypergraphs with a fixed regular degree distribution, and thus fixed 〈k(2)〉 and δ(2), but with different degrees of intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\), as illustrated by the vertical dashed line in Fig. 1d.

Effect of intra-order overlap on complex contagion

In order to investigate how the microscopic organization of group interactions affects the dynamics of complex systems, we focus on two different types of dynamical processes on hypergraphs, namely complex contagion and synchronization of phase oscillators. As the first case study, we consider a Susceptible-Infected-Susceptible (SIS) compartmental model, which has been widely used to mimic the spread of ideas through social contagion in hypergraphs25,27. In this framework, each node of the hypergraph represents an individual that can be either susceptible (S) or infected (I). The transition from the state S to I occurs when a susceptible individual enters in contact with infected ones, via a pairwise or a three-body interaction. More specifically, a susceptible individual can become infected, with probability β(1), through a 1-hyperedge (a link) with an infected individual, and with probability β(2), through a 2-hyperedge (an interaction in a group of three) with two other infected nodes. As in the standard SIS model, an infected individual can recover with probability μ. In absence of higher-order interactions, the SIS model undergoes a continuous transition from an absorbent state with a vanishing stationary fraction of infected individuals, ρ⋆ =  0, to an active state with ρ⋆ ≠ 088. In previous works25,40, it has been shown that the presence of higher-order interactions changes the nature of the transition from continuous to discontinuous, hence leading to the emergence of a region of bistability, where both the absorbent and active states co-exist. Here, we show that this does not occur for any hypergraph given a fixed interaction strength, but strongly depends on the level of intra-order hyperedge overlap.

Figure 2 shows the results of stochastic simulations of the SIS model on hypergraphs with tunable intra-order hyperedge overlap. Figure 2a reports the fraction ρ⋆ of infected individuals in the stationary state as a function of \({{{{\mathcal{T}}}}}^{(2)}\) and of the rescaled transmission probability λ(1) = 〈k(1)〉β(1)/μ, while keeping the rescaled transmission probability through 2-hyperedges fixed to λ(2) = 〈k(2)〉β(2)/μ = 3. The phase diagram of Fig. 2a was obtained by averaging M = 200 simulations, half of which started with a small fraction of infected individuals ρ(0) = 0.01, and the other half with ρ(0) = 0.8. We observe the presence of three phases, with a region of bistability that only appears for small values of \({{{{\mathcal{T}}}}}^{(2)}\). An example of the explosive transition associated to bistability is shown in Fig. 2b, which illustrates ρ⋆ vs. λ(1) for \({{{{\mathcal{T}}}}}^{(2)}=0.1\). The three phases merge at a tricritical point, leading to the disappearance of the bistability. Beyond this point, the transition from the absorbent to the active phase is continuous. As an example, Fig. 2c illustrates the continuous behavior of ρ⋆ vs. λ(1) in a hypergraph with \({{{{\mathcal{T}}}}}^{(2)}=0.8\). In this case, the high overlap of the hyperedges of a node reduces the effectiveness in propagating its infected state.

Fig. 2: Effect of intra-order hyperedge overlap on complex contagion.
figure 2

a Phase diagram for the SIS model on a hypergraph with N = 1000 nodes and average generalized degrees k(1) = 5 and k(2) = 6. The value of 2-hyperedge infectivity is set to λ(2) = 3, while the value of the recovery rate is fixed to μ = 0.05. Three phases emerge as a function of λ(1) and of the hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\): an absorbent phase with ρ⋆ = 0, an active phase with an endemic stationary state ρ⋆ ≠ 0, and a bistability phase, where the stationary state depends on the initial conditions. Two cuts of the diagram: panel b, for \({{{{\mathcal{T}}}}}^{(2)}=0.1\), shows an explosive transition, while panel c, for \({{{{\mathcal{T}}}}}^{(2)}=0.8\), shows a continuous transition. Note that in panel b the solid and dashed lines respectively correspond to ρ(0) = 0.01 and ρ(0) = 0.8 initial conditions.

Full size image

Effect of intra-order overlap on synchronization

To gain a deeper insight on how hyperedge overlap affects the emergence of collective states, we consider synchronization as a second case study. Specifically, we investigate a system of Kuramoto oscillators61 coupled with both pairwise and three-body interactions. The Kuramoto model has been proposed to study synchronization dynamics in biological systems, from pacemaker cells to crowds of fireflies89, but has also been used to analyze synchronization in the brain at the cortical level62,63,64, an emergent phenomenon that plays an important role in many brain functions50,51,52,53,54,55,56,57,58,59,60,90. When implemented on a network the Kuramoto model assumes that each node i (i = 1, …, N) of the network is a phase oscillator and is characterized, at time t, by its phase θi(t). In the case when the oscillators are coupled through pairwise and three-body interactions, the equations governing the time evolution of the oscillator phases read91,92:

$${\dot{\theta }}_{i}= {\omega }_{i}+\frac{{\sigma }^{(1)}}{\langle {k}^{(1)}\rangle }{\sum}_{j=1}^{N}{a}_{ij}^{(1)}\sin ({\theta }_{j}-{\theta }_{i}) \\ +\frac{{\sigma }^{(2)}}{2!\langle {k}^{(2)}\rangle }{\sum}_{j=1}^{N}{\sum}_{k=1}^{N}{a}_{ijk}^{(2)}\sin ({\theta }_{j}+{\theta }_{k}-2{\theta }_{i}),\qquad i=1,\ldots,N$$
(4)

where ωi is the natural frequency associated to oscillator i, which is randomly sampled from a uniform distribution. The quantities σ(1) and σ(2) represent the coupling strength for pairwise (order 1) and three-body (order 2) interactions, respectively, while \({a}_{ij}^{(1)}=1\) if there is a link between node i and j, and \({a}_{ijk}^{(2)}=1\) if the hypergraph has a hyperedge containing i, j and k. We remark that various extensions of the Kuramoto model to hypernetworks, mainly differing for the functional form of the higher-order interactions, have been presented29,30,91. Note that the version of the Kuramoto model with three-body interactions in Eq. (4) that we adopt in our study is the one introduced in ref. 91. In a similar way to the original Kuramoto model on networks, the synchronization transition is captured by the order parameter defined as:

$$r(t)={N}^{-1}\left| {\sum}_{j=1}^{N}{e}^{{{{\rm{i}}}}{\theta }_{j}(t)}\right| \,$$
(5)

ranging from 0 (incoherence) to 1 (full synchronization). With only pairwise interactions, and in the absence of additional ingredients9, a continuous synchronization transition takes place when the coupling strength σ(1) exceeds a critical value, \({\sigma }_{c}^{(1)}\)93.

As in the case of the SIS dynamics, we have studied this model on hypergraphs with tunable intra-order hyperedge overlap. Figure 3 shows the results of the numerical integration of Eq. (4) for σ(2) = 3 and for various values of σ(1) and \({{{{\mathcal{T}}}}}^{(2)}\). For each set of parameters, we perform M = 200 simulations, half with a narrow (Δθ = 0.02) and half with a wider (Δθ  = 2π) distribution of initial phases. For each simulation, we let the system evolve a transient time tr and calculate \(\langle r\rangle=\frac{1}{\Delta t}\int_{{t}_{r}}^{{t}_{r}+\Delta t}r(t)dt\), where Δt → ∞. We then consider the median of 〈r〉 over all runs.

Fig. 3: Effect of intra-order hyperedge overlap on synchronization.
figure 3

a Phase diagram for the Kuramoto model on a hypergraph with N = 1000 nodes and average generalized degrees k(1) = 5 and k(2) = 6. The coupling strength for three-body interactions is set to σ(2) = 3. Three phases emerge as a function of σ(1) and \({{{{\mathcal{T}}}}}^{(2)}\): an incoherent phase with low values of 〈r〉, a synchronized phase with large 〈r〉, and a bistability phase where the system can be synchronized or not depending on the initial conditions. b–g Two cuts of the diagram: panels b–d, for \({{{{\mathcal{T}}}}}^{(2)}=0.06\), show an explosive transition, while panels e–g, for \({{{{\mathcal{T}}}}}^{(2)}=0.6\), show a continuous transition. Panels b and e illustrate the order parameter 〈r〉, while panels c and f the local synchronization parameter 〈rloc〉. Panels d and g show the average effective frequencies \(\langle {\dot{\theta }}_{i}\rangle\) of the hypergraph nodes in the forward  continuation branch. Note that in panels b, c and e, f the solid (dashed) branch stands for the forward (backwards) continuation.

Full size image

The phase diagram in Fig. 3a reporting 〈r〉 as a function of \({{{{\mathcal{T}}}}}^{(2)}\) and σ(2) resembles that of the SIS model in Fig. 2a. Again, the region of bistability is present for small \({{{{\mathcal{T}}}}}^{(2)}\), but disappears for large values. Figure 3b, e show 〈r〉 as a function of σ(1), respectively for \({{{{\mathcal{T}}}}}^{(2)}=0.06\), where the transition to synchronization is explosive, and for \({{{{\mathcal{T}}}}}^{(2)}=0.6\), where the transition is continuous. To gain a better insight into the microscopic mechanisms underlying the two different transitions to synchronization of Eq. (4), we introduce the following measure of higher-order local synchronization:

$$\langle {r}_{loc}\rangle= \frac{1}{{\sum}_{i=1}\left({\sum}_{j=1}{a}_{ij}^{(1)}+\frac{1}{2!}{\sum}_{j=1}{\sum}_{k=1}{a}_{ijk}^{(2)}\right)} \times {\sum}_{i=1}\left({\sum}_{j=1}\left| {\lim}_{\Delta t\to \infty }\frac{{a}_{ij}^{(1)}}{\Delta t}{\int}_{{t}_{r}}^{{t}_{r}+\Delta t}{e}^{{{{\rm{i}}}}[{\theta}_{j}(t)-{\theta }_{i}(t)]}dt\right| \right.\\ \left.+\frac{1}{2!}{\sum}_{j=1}{\sum}_{k=1}\left| {\lim}_{\Delta t\to \infty }\frac{{a}_{ijk}^{(2)}}{\Delta t}{\int}_{{t}_{r}}^{{t}_{r}+\Delta t}{e}^{{{{\rm{i}}}}[{\theta }_{j}(t)+{\theta }_{k}(t)-2{\theta }_{i}(t)]}dt\right| \phantom{\frac{1}{2}}\right).$$
(6)

This definition of 〈rloc〉 extends the one used for network edges94, since it gathers the phase coherence between pairs of nodes connected by links, as well as among triplets of nodes connected by 2-hyperedges. As for the order parameter 〈r〉, full synchronization corresponds to 〈rloc〉 = 1, while complete incoherence yields 〈rloc〉 = 0.

Figure 3c, f show that the intra-order hyperedge overlap radically changes the behavior of 〈rloc〉 as a function of σ(1). In agreement with 〈r〉, also 〈rloc〉 indicates that the transition is explosive for hypergraphs with low intra-order hyperedge overlap, while it is continuous for hypergraphs with high intra-order hyperedge overlap. However, at variance with 〈r〉, for large values of σ(1), 〈rloc〉 reaches the value 〈rloc〉 = 1, meaning that, even if the system is not globally synchronized, there is still a high degree of local synchronization. The complete phase diagram for 〈rloc〉 as function of σ(1) and \({{{{\mathcal{T}}}}}^{(2)}\) is reported in Supplementary Fig. 1.

The high degree of local synchronization is a consequence of the formation of synchronous clusters in the system, and depends on the level of intra-order hyperedge overlap. This effect becomes evident from the analysis of the average effective frequencies of the oscillators, defined as:

$$\langle {\dot{\theta }}_{i}\rangle=\frac{1}{\Delta t}\int_{{t}_{r}}^{{t}_{r}+\Delta t}{\dot{\theta }}_{i}(t)dt\,$$
(7)

where Δt → ∞. In Fig. 3d, g, we compare the effective frequencies \(\langle {\dot{\theta }}_{i}\rangle\) in hypergraphs with different levels of intra-order hyperedge overlap. When \({{{{\mathcal{T}}}}}^{(2)}\) is small, as in Fig. 3d, the 2-hyperedges have few nodes in commmon and this makes difficult the nucleation of synchronous clusters. Consequently, the synchronization onset is abrupt and, as we increase σ(1), the effective frequencies \(\langle {\dot{\theta }}_{i}\rangle\) converge suddenly to their mean value. Conversely, when \({{{{\mathcal{T}}}}}^{(2)}\) is large, as in Fig. 3g, the 2-hyperedges share multiple nodes, which promotes the formation of synchronous clusters. The result is that, as σ(1) increases from zero, we find several groups of nodes locked at the same effective frequency. As σ(1) further increases, these groups smoothly merge together, until the whole system becomes synchronized.

Discussion

In this work, we have investigated how the microscopic organization of higher-order interactions affects the dynamical properties of a complex system. To this aim we have purposely introduced a way to measure the overlap of the group interactions to which a node participates. Such an overlap is quantified by counting and properly normalizing the number of shared nodes among hyperedges of the same order. In this way, our metric allows for the differentiation of hypergraphs having the same number of hyperedges of any size, the same generalized node-connectivity, but different microscopic organization. Our intra-order overlap metric ranges from its minimal value of 0, in structures where nodes have non-coincident neighbors in their groups interactions, to its maximal value of 1, when the different groups of the same size share the maximum possible number of nodes. At variance with the metric proposed by Yin et al.95, which, for pairwise networks, extends the concept of clustering coefficient beyond three-body motifs, our metric applies to structures with truly higher-order interactions. Moreover, our metric can be used in all types of hypergraphs, unlike the framework by Kartun-Giles et al.96, which is specifically tailored to the study of geometrical and topological properties of simplicial complexes. Additionally, we focus on each order (size) of hyperedges in a separate way, in contrast to the approach by Lee et al.43 in which a single measure is evaluated for the whole structure considering all hyperedges, regardless of their size. Our metric also differs from the approach of refs. 97,98,99, since counting the number of pairs of nodes appearing in two different hyperedges corresponds to counting the number of cycles of the possible shortest length, namely four, in a bipartite network connecting nodes and hyperedges. Finally, our metric is complementary to other microscopic/local metrics of inter-order hyperedge overlap and/or downward closure45,46,47, and to global approaches that instead try to capture macroscopic features such as the presence of a giant higher-order component42.

We have first used our measure to study hypergraphs describing higher-order interactions in real-world complex systems. We found that, together with a broad range of values of measures of correlations among hyperedges, such as nestedness and inter-order hyperedge overlap (inclusion or simpliciality)43,44,100, complex systems also show a diverse range of intra-order hyperedge overlap. This indicates that in the real world the likelihood of two distinct hyperedges of the same order sharing more than one node is not negligible and, therefore, the way e.g., three-body interactions are organized and correlated can play an important role for the function and dynamics of a complex system.

Traditional random models for hypergraphs, however, do not allow to take this aspect into account. In a synthetic random hypergraph, hyperedges are randomly assigned and the likelihood of two distinct hyperedges of order two sharing more than one node becomes negligible. This arises from the fact that the number of hyperedges scales with the size of the structure \({{{\mathcal{O}}}}(N)\), while the cardinality of these hyperedges remains orders of magnitude smaller, i.e., \({{{\mathcal{O}}}}(1)\)46,99. To overcome this limitation, we developed a model to construct synthetic hypergraphs with fixed node hyperdegrees and a tunable value of the intra-order hyperedge overlap. This allowed us to examine, in a controlled and systematic way, the effects of intra-order hyperedge overlap in shaping the onset of collective phenomena in systems with higher-order interactions. We have found that the hyperedge overlap has the same effect on two different types of dynamical processes, namely on complex contagion and on synchronization. In both processes, it drives the emergence of collective behaviors and also determines the nature of the transition, with lower values of overlap in general facilitating explosivity and bistability. Examples of simulations of these two dynamics on hypergraphs describing the structure of real-world systems can be found in Supplementary Figs. 2 and 3. These results open the way for further research to investigate the role of intra-order hyperedge overlap in other dynamical processes where higher-order interactions can play a role, such as in the dynamics of ecosystems101, in evolutionary game theory33,34,35 in cluster synchronization102,103 or in interacting contagion dynamics104,105. The spatial nature of most of those systems suggests that interactions of the same order can inherently overlap, thus providing a particularly relevant case study for the analysis of the effect of correlations.

Another promising avenue for future research is to explore the combined effect on dynamics of intra-order hyperedege overlap and other structural features both at the local46,47 and at the global scale42. In particular, Burgio et al.47 have demonstrated how correlations between different orders of interaction can also hinder bistability in social contagion dynamics, while Kim et al.42 have shown that real structures with a giant higher-order component, whose presence necessitates some degree of intra-order hyperedge overlap, ensure the stability of the active phase of social contagion processes. Integrating insights from these studies could enhance our understanding of the emergence of collective behavior in complex systems with higher-order interactions.

Regarding the limitations of our study, we must mention the scarcity of data about real-world hypergraphs. This is a common issue that affects the present development of research about dynamical behaviors of coupled systems through higher-order interactions. Accessing data on true higher-order interactions rather than relying on their construction from pairwise contact information poses a challenge that network science must address in the coming years.

Overcoming the above experimental challenge will serve not only to unveil the true higher-order skeleton of interactions of complex systems, but also to be able to translate current theoretical results on the dynamics emerging from group interactions into practical applications to biological, technological and social systems. In this sense, our findings can serve as the starting point for the design of methods to control the observed dynamics by manipulating the structure of group interactions106,107,108,109,110. However, for practical applications, there is still a theoretical gap to be filled, as the number of control actions to be exercised may limit the feasibility of the approach. One possible way to mitigate this limitation could be the formulation of an optimization problem, which explicitly tries to minimize the changes in the structure while meeting the objective of manipulating, for example by minimizing or maximizing, the hyperedge overlap of the system.

Methods

Maximum and minimum number of unique neighbors of a node

The node intra-order hyperedge overlap \({T}_{i}^{(m)}\) is constructed for each node i by counting the number of unique nodes that belong to the set of hyperedges \({{{{\mathcal{E}}}}}_{i}^{(m)}\) to which node i participates. To properly normalize the metric within the range [0,1], it is necessary to consider the minimum, \({S}_{i}^{(m),-}\), and maximum, \({S}_{i}^{(m),+}\), number of unique neighbors that node i can have through its \({k}_{i}^{(m)}\) hyperedges of order m (here \({k}_{i}^{(m)}\) indicates the generalized degree of order m for the node i). Both quantities can be analytically calculated for every generalized degree \({k}_{i}^{(m)}\) of node i and for every order of the interactions m. In particular, the maximum amount of unique nodes is given by

$${{{{\mathcal{S}}}}}_{i}^{(m),+}=m{k}_{i}^{(m)}.$$
(8)

To obtain the minimum number of unique neighbors, \({S}_{i}^{(m),-}\), we follow these steps. Let us, for a moment, indicate this minimum number with n. These n unique neighbors are connected by a number \({k}_{i}^{(m)}\) of hyperedges of order m. At this point, n can be calculated by using the following implicit relation

$$\left(\begin{array}{c}n\\ m\end{array}\right)={k}_{i}^{(m)}.$$
(9)

On the left hand side, we have the number of combinations of neighboring nodes that can be built in groups of size m, and on the right hand side the number of hyperedges of the focal node i. Since the solution could give a non-integrer number, we take the upper integer, as we need to consider the lowest number of nodes that is bigger than the obtained required quantity. Thus, we have that \({{{{\mathcal{S}}}}}_{i}^{(2),-}=\lceil n\rceil\). For the case m = 2, Eq. (9) becomes the following second-order equation in the unknown n:

$${n}^{2}-n-2{k}_{i}^{(2)}=0,$$
(10)

whose only physically valid solution is given by

$$n=\frac{1+\sqrt{1+8{k}_{i}^{(2)}}}{2}.$$
(11)

An illustration of two extreme configurations is presented in Fig. 1a, c, where we consider the local intra-order hyperedge overlap of a node i with a generalized 2-degree \({k}_{i}^{(2)}=4\). In Fig. 1a, we show the configuration in which the number of unique neighbors of node i is at its maximum, i.e., \({{{{\mathcal{S}}}}}_{i}^{(2)}={{{{\mathcal{S}}}}}_{i}^{(2),+}=8\), resulting in \({T}_{i}^{(2)}=0\). In contrast, in Fig. 1c, we consider the configuration where the number of unique neighbors of node i is at its minimum, i.e., \({{{{\mathcal{S}}}}}_{i}^{(2)}={{{{\mathcal{S}}}}}_{i}^{(2),-}=4\), leading to \({T}_{i}^{(2)}=1\).

Hypergraphs of real-world systems

Depending on the nature of the data, the 39 hypernetworks of real-world systems we have analyzed can be classified in three main groups: (i) FACE-TO-FACE CONTACTS (a total of 7 hypernetworks) obtained from empirical experiments conducted in different social environments, namely a hospital68, a scientific conference67, a high school69, a workplace66, a primary school70, a Malaui village71 and a university72. (ii) BRAIN CONNECTOMES (a total of 12 hypernetworks) of humans and other animals mapped through different acquisition techniques, namely the human cerebral cortex73, the mouse visual cortex74, the worm hermaphrodite C. elegans anterior section75,76, the worm adult male C. elegans posterior nervous system77, the Norway rat cortex connectome78 and the mouse brain connectome79 (iii) OTHER SYSTEMS (for a total of 20 hypernetworks) collected from the Austin R. Benson’ repository (the link to the repository can be found in the data availability section) and including two coauthorship datasets80,81, three bill co-sponsorship datasets from the US congresspersons and senators80,82,83,84, two datasets of questions answered by users on StackOverflow and mathoverflow85, a dataset of products bought on walmart shopping trips86, a dataset of products reviewed by users in Amazon87, a dataset of hotels clicked on during a browsing session on the trivago website83, a datasets of drugs used by patients80, a dataset of classifications applied to drugs80, a dataset of substances making up drugs80, a dataset of drugs used by patients recorded in emergency room visits80, three datasets of tags applied to question in online websites80, and two datasets of users asking and answering questions on threads in online websites80.

Considering that the original contact data datasets encoded temporal face-to-face contacts and that the original brain connectome datasets encoded weighted pairwise networks, a pre-processing was required to construct the hypernetworks. Regarding the face-to-face contact data, the original datasets have a temporal resolution, whereby all pairwise face-to-face contacts detected by the sensors are recorded at 20-second intervals (or at one-minute intervals, depending on the dataset). To transform the temporal contact datasets into static hypergraphs, following Iacopini et al.25 we initially construct a set of static networks by aggregating all the contacts that occur within a static temporal window of length Δt = 5. That way, the first network is constructed by aggregating all contacts occurring between time t = t0 and t =  t0 + Δt, while the second network is formed by aggregating all contacts occurring between t = t0 + Δt and t = t0 + 2Δt. This process is repeated for each subsequent network. Afterward, we build an static hypergraph for each of the static networks by considering as hyperedges all the cliques of the aggregated network, since they are groups that have been coincident in time. Next, we buid a weighted hypergraph by counting the number of times that each clique is present on each of the static hypergraphs. Finally, to extract the meaningful structure of the interactions as a hypernetwork, we keep the 10% of hyperedges with the highest weights. Regarding the brain connectomes, the networks extracted from the datasets are all very dense pairwise and weighted networks. Therefore, we have first created unweighted networks of the most important relations among brain regions by only keeping the 10% of connections with the highest weights. Afterwards, to build the hypernetwork, following Zhang et al.45 we transformed all the cliques of the unweighted network into hyperedges.

Synthetic hypergraphs

To construct synthetic hypergraphs with tunable intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\), we start from a configuration with maximum intra-order hyperedge overlap, i.e., \({{{{\mathcal{T}}}}}^{(2)}=1\). Then, the value of \({{{{\mathcal{T}}}}}^{(2)}\) in the model is tuned by rewiring a fraction of the existing 2-hyperedges. In more detail, we start from a structure with a natively maximized value of the intra-order hyperedge overlap, which we name the Locally Overlapped Regular Simplicial Complex (LORSC). This is a special type of hypergraph obtained with the following procedure. We first fix a number \(N^{\prime}\) of nodes and the connectivity of the final strucutre k(1). Then, we substitute each of the nodes with a set of k(1) nodes, and we add links between these nodes so that to obtain fully connected 1-cliques within each set. Finally, in order to have a connected structure, we join each node with a randomly selected node belonging to another clique. In this way, we obtain a structure with a number \(N={k}^{(1)}N^{\prime}\) of nodes, each of them having exactly k(1) 1 − hyperedges. Once the backbone of pairwise hyperedges is constructed, we consider as a 2 − hyperedge each of the closed triplets found in the structure. Consequently, each node in the structure will not only have the same number of 1 − hyperedges, i.e., k(1) but also the same number of 2 − hyperedges, i.e., k(2).

The next step in the construction of a structure with tunable global intra-order hyperedge overlap \({{{{\mathcal{T}}}}}^{(2)}\) is to introduce rewiring, such that to reduce the value of \({{{{\mathcal{T}}}}}^{(2)}\). In more detail, at each iteration we randomly select two 2-hyperedges from the set \({{{{\bf{{{{\mathcal{E}}}}}}}}}^{(2)}\). Afterward, we select one node from each of the two previously random selected 2-hyperedges. Then, we swap the selected nodes and calculate \({{{{\mathcal{T}}}}}^{(2)}\) in the original structure and in the structure after swapping. We keep the change only if the global intra-order hyperedge overlap decreases or remains unchanged, otherwise we proceed using the structure before the swapping. By doing so, we are able to reduce \({{{{\mathcal{T}}}}}^{(2)}\) while maintaining the value of the degree \({k}_{i}^{(2)}\) for all nodes. For more details on the synthetic structures, as well as a graphical example of the structures built, we refer the reader to Supplementary Fig. 4.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.