A unified framework for identifying influential nodes in hypergraphs

In this section, we will utilize two representative hypergraph-based propagation models to elucidate the complex effects of dynamical processes on node influence, a finding that motivates this study. After that, we propose a node influence evaluation framework that concisely and effectively integrates the underlying dynamics.

A hypergraph is generally given as $\mathcal{H}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}=\{v_{1},v_{2},...,v_{N}\}$ represents the set of $N$ nodes and $\mathcal{E}=\{h_{1},h_{2},...,h_{M}\}$ is the set of hyperedges. Here, a hyperedge $h$ contains a set of nodes, and $m=|h|$ is its size. Intuitively, a hyperedge denotes the interaction among a set of nodes (such as a chat group), which is a general extension of the normal pairwise edges where $m=2$ Boccaletti et al. (2023).

To reveal the variance of node influence ranking under different dynamics, we consider two widespread propagation processes, the higher-order contagion model with power-law infection kernel (HCP model) St-Onge et al. (2022); Iacopini et al. (2019) and the higher-order threshold contagion model (HTC model) de Arruda et al. (2020); Ferraz de Arruda et al. (2023) (see details in Methods). Using a widely used hypergraph constructed from the Email-Enron dataset as the underlying topology, we identify the ground-truth influential nodes under each of the two hypergraph dynamics models by measuring the total infection size through extensive simulations, following standard practice in prior works Kitsak et al. (2010); Mancastroppa et al. (2023). The details of the simulations and the dataset can be seen in Methods. We evaluate three groups of propagation parameter settings for each model, denoted as CP₁–CP₃ for the HCP model and TP₁–TP₃ for the HTC model. Figure 1a-c highlight the top 25% influential nodes under different dynamical mechanisms and parameters: Fig.1a corresponds to the HCP mechanism with parameter setting CP₁, while Fig.1b and c correspond to the HTC dynamics with parameter settings TP₁ and TP₂. Notably, the influential nodes and their importance significantly change under different dynamical settings. This indicates that both dynamical mechanisms and parameters play a pivotal role in identifying key nodes. For a more comprehensive and intuitive illustration, we arranged all nodes in a single sequence, using different colors to represent their influence. In Fig. 1d and Fig. 1e, CP₁ and TP₁ respectively serve as the baselines, where the nodes are ordered according to their influence scores. The ordering-similarity among CP₁, CP₂, and CP₃ reveals stable node influence patterns with no significant fluctuations. However, notable differences emerge when comparing CP₁ with TP₁, as well as among TP₁, TP₂, and TP₃. Such observations indicate that the sensitivity of node influence to parameter changes varies across different higher-order dynamical processes. A natural question, then, is how to define a general measure of node influence that integrates dynamical information and remains effective under diverse higher-order interactions.

To begin with, we consider the standard $S$–$I$–$R$ notation, which has been widely used to describe various real-world spreading dynamics Boccaletti et al. (2023); Mancastroppa et al. (2023). Generally, nodes that have not yet undergone propagation are denoted as $S$ (susceptible), nodes that have been activated and can propagate information, disease, or opinions are marked as $I$ (infected), and nodes that are immune to the propagation process and no longer participate in it are labeled as $R$ (recovered or removed). Denote the number of nodes in state $S$, $I$ and $R$ at time $t$ as $S_{t}$, $I_{t}$, and $R_{t}$, respectively. The essence of evaluating node influence lies in establishing a connection between the seed node and its eventual propagation outcome, i.e., the expected number of removed nodes at the end of the process, denoted as $\mathbb{E}[R_{\infty}|s]$, where $s$ is a given seed. In fact, regardless of the specific dynamical mechanism, the final spreading range can be viewed as the accumulation of infected nodes over time, and is progressively approximated by the sequence $\{\mathbb{E}[I_{t}+R_{t}|s]\}$ (for $\lim_{t\to\infty}\mathbb{E}[I_{t}|s]=0$ in SIR models). For convenience of analysis, we suppose the recovery rate from $I$ to $R$ is $\mu=1$, and we have:

See Fig. 1f for an instantiation. For cases when $\mu\neq 1$, we only need to refine $I_{t}$ as the number of newly infected nodes at time $t$. Note that each propagation step relies on the former step and ultimately traces back to the original seed node. Intuitively, the early-stage propagation inherently encodes both the seed’s topological position and the system’s dynamical characteristics, which in turn, largely determine the final propagation outcome. Therefore, we approximate the ultimate influence of each seed using its early-stage propagation range, forming the basis of our IPS (Initial Propagation Score) framework. Formally, given a dynamical process $P$ and a target time step $t$, the IPS measure is denoted as $IPS_{t}^{P}=\mathbb{E}[I_{t}+R_{t}\mid s]$ (also called the $t$-th order IPS).

In this section, we adopt the HCP model to demonstrate the calculation process and show that even the first-order IPS ($t=1$) could exhibit strong predictive power in real-world systems. At $t=1$, propagation can only occur within the first-order neighborhood of the seed node $s$, and the expected number of newly infected nodes is the sum of infection probabilities over its neighbors. In HCP, the infection rate in a hyperedge with $i$ infected nodes is given as $\lambda i^{\nu}$, where $\lambda$ is a global propagation rate parameter and $\nu$ reflects the reinforcement effect. Since only the seed node is infected at $t=0$ and all other nodes are susceptible, the probability that a neighboring node becomes infected via a shared hyperedge is $1-e^{-{\lambda}}$ (see details in Methods). Due to the independence of infection events across hyperedges, the probability that a neighbor $i$ of $s$ becomes infected at $t=1$ is $p_{1}(i\mid s)=1-e^{-k_{s}(i){\lambda}}\approx k_{s}(i)\lambda$ ($\lambda$ is typically small), where $k_{s}(i)=|\{h\mid s,i\in h,h\in\mathcal{E}\}|$ is the number of hyperedges adjacent to $s$ that contain $i$. Then the expected propagation range (number of infected nodes and recovered nodes) of the seed node $s$ at $t=1$ is:

where $N_{v}^{1}(s)=\bigcup_{s\in h,h\in\mathcal{E}}h\setminus\{s\}$ is the set of first-order neighbors of $s$, $N_{h}^{1}(s)=\{h|s\in h,h\in\mathcal{E}\}$ is the set of adjacent hyperedges of $s$. Eq.(2) indicates that the ranking of $IPS_{1}^{HCP}$ is solely determined by the sum of the sizes of its neighboring hyperedges, leading to a computational complexity of $\mathcal{O}(N\langle k\rangle$ that is highly efficient for large-scale applications (Here $\langle k\rangle$ is the average hyper-degree of the hypergraph, and the details for time complexity are provided in Supplementary Section S1.3). We also notice that the ranking of $IPS_{1}^{HCP}$ is independent of the dynamical parameters, which is consistent with the previously observed robustness of influence rankings across different HCP parameters in Fig. 1d.

We evaluate the predictive performance of $IPS_{1}^{HCP}$ on 20 real-world higher-order datasets from various domains, including human contact networks, email communications, legislative collaboration, online platforms, and ecological systems. These datasets vary in size, density, and functional context. A complete summary of the datasets is provided in Methods.

For comprehensive benchmarking, we include 10 existing centrality measures widely used in higher-order interactions, including the basic hyper-degree, the recently proposed hyper-coreness-$R$, hyper-coreness-$R_{w}$ Mancastroppa et al. (2023), the Node-and-edge-nonlinear-eigenvector centrality measures (abbreviated as eigenvector-linear and eigenvector-max), and the natural extension of five classic measures to the 2-projection networks of hypergraphs, containing degree, k-coreness, eigenvector, betweenness, and closeness. To distinguish these projection-based measures from their standard counterparts in pairwise networks, we prefix them with ‘2-’ (e.g., 2-degree, 2-betweenness, etc.).

Here we use different indicators to evaluate and compare measures’ performance. Specifically, we first use Kendall’s $\tau$ to quantify the consistency between the node influence ranking based on different measures and the ground truth derived from over 300 times simulations, with a larger $\tau$ (closer to 1) indicating better performance of the measure. Figure 2a clearly shows that the $IPS_{1}^{HCP}$ method, corresponding to the red box, has the best overall performance, which even outperforms the second-best algorithm by around 10%. We then unfold the results on each hypergraph in Fig. 2b. Notably, the $IPS_{1}^{HCP}$ method consistently performs the best in almost all datasets (red points), while other measures show large fluctuations across different hypergraphs. We also employ the Jaccard similarity, which quantifies the overlap between the predicted top $\lfloor rN\rfloor$ influential nodes by centrality measures and the ground truth, to further evaluate the prediction accuracy. Figure 2c–i present results on several representative hypergraphs (highlighted in gray in Fig. 2b), and Fig. 2j additionally shows a case on a large-scale hypergraph. Again, the $IPS_{1}^{HCP}$ method consistently achieves the highest similarity across nearly all settings. See the results for other hypergraphs in Supplementary Figures S3, S4, and S6.

Taken together, these results validate the superior stability and effectiveness of the $IPS_{1}^{HCP}$ method across a wide range of real-world hypergraphs. Other methods, by contrast, exhibit considerable variability, with performance often depending on both the network structure and the proportion $r$. For example, hyper-coreness-$R_{w}$ and eigenvector-linear occasionally achieve high performance within certain ranges of $r$, but display sharp declines as $r$ increases. Moreover, while they often appear as the second-best performers (e.g., Fig. 2e, f, g), they are outperformed by hyper-degree in several hypergraphs (e.g., Fig. 2c, h, i). Furthermore, the general trend that hyperstructure-based methods outperform projection-based counterparts is reversed in the large-scale threads-math-sx hypergraph, where 2-k-coreness and 2-degree match the performance of $IPS_{1}^{HCP}$. Notably, several methods could not even be evaluated on this hypergraph due to the excessive computational cost.

In Fig. 3, we further analyze the consistency between the $IPS_{1}^{HCP}$ method and other baseline methods. We plot each method’s Kendall correlation with $IPS_{1}^{HCP}$ ($\tau_{m}$) against its predictive performance ($\tau$, as defined in Fig. 2a, b). Each point represents a specific method evaluated on a certain hypergraph, with color-coding consistent with the legend in Fig. 2. We find a pronounced positive correlation, whereby methods that exhibit stronger agreement with $IPS_{1}^{HCP}$ tend to achieve higher accuracy. Such a clear alignment, which is not observed in other measures (see Supplementary Figure S5, offers a unified explanation for the performance fluctuations of different algorithms across different higher-order structures. For instance, the highlighted nodes and the subgraphs in Fig. 3 illustrate why some baseline methods outperform others in Fig. 2e, i, and j. More importantly, this suggests that our method may have captured the fundamental features governing node importance, namely, the coupling between local topology and early dynamics.

We incorporate the second-order structural and dynamical information, and extend $IPS_{1}^{HCP}$ to $IPS_{2}^{HCP}$ under the HCP model, assessing whether such enrichment can substantially improve predictive accuracy. Furthermore, we will systematically analyze the robustness of the IPS methods under varying dynamical parameters.

Following Eq. (1), when considering second-order information ($t=2$), we have:

Here $\mathbb{E}[R_{1}+I_{1}\mid s]$ is given by Eq.(2), and $\mathbb{E}[I_{2}\mid s]$ denotes the number of newly infected nodes at time $t=2$. Let $p_{t,n}(i|s)$ represents the probability that the $n$-th order neighbor $i$ of the seed $s$ is infected at time $t$. Then, $\mathbb{E}[I_{2}\mid s]$ can be derived by

where the first item is the expected number of first-order neighbors that are infected at time $t=2$, and the remaining item records the expected number of second-order infected neighbors.

We then need to compute $p_{t,n}(i|s)$. Let $N_{h}^{1}(s)$ and $N_{h}^{2}(s)=\bigcup_{i\in N_{v}^{1}(s)}N_{h}^{1}(i)\setminus N_{h}^{1}(s)$ represent the set of first-order adjacent hyperedges and second-order hyperedges, respectively. Moreover, the expected number of infected nodes within hyperedge $h$ can be approximated as follows: If $i\in N_{v}^{1}(s)$, $I_{1}(i,h)=\sum_{j\in h\cap N_{v}^{1}(s),j\neq i}p_{1,1}(j|s)$; Otherwise, $I_{2}(h)=\sum_{j\in h\cap N_{v}^{1}(s)}p_{1,1}(j|s)$. Finally, through enormous computation (see details in Supplementary Section S1.1), we can obtain the expression of $p_{t,n}(i|s)$, i.e.,

with:

See Supplementary Section S1.1 for an example in a hypergraph.

It is noteworthy that $IPS_{2}^{HCP}$ exhibits a complexity of $\mathcal{O}(N\langle k\rangle k_{\max}m_{\max}^{2})$, where $k_{\max}$ is the maximal hyper-degree and $m_{\max}=\max_{h\in\mathcal{E}}|h|$ is the maximal size of hyperedges. Similar to Eq. (2), through infinitesimal substitution that transforms node-level calculations into the hyperedge-level, we can make additional approximations and obtain a faster version with time complexity reduced to $\mathcal{O}(N\langle k\rangle k_{\max}m_{\max})$, denoted as $IPS_{2r}^{HCP}$. Details are provided in Supplementary Section S1.2 and S1.3. When $k_{\max}$ and $m_{\max}$ are bounded (e.g., sparse hypergraphs), the complexity of both versions scales linearly with node number $N$.

We conduct extensive experiments with a large range of dynamical parameters to validate the performance and robustness of IPS methods systematically. We consider cases where the recovery probability $\mu=0.1,0.5,1$, and the reinforcement parameter $\nu=1,2,3,4$, respectively. For the contagion parameter $\lambda$, we experimentally calculate the propagation threshold $\lambda_{c}$ for each network with given $\mu$ and $\nu$ (see Methods) and consider $0.1\lambda_{c}$, $\lambda_{c}$, and $10\lambda_{c}$ respectively. In summary, we examine $20$ real hypergraphs, each with $36$ groups of different dynamical parameters. We employ Kendall’s $\tau$ as before to evaluate the effectiveness of different methods.

Figure 4 presents the results for the case of $\mu=1$. Overall, the IPS methods outperform all other baselines and exhibit strong robustness to parameter variations. Among them, the approximation $IPS_{2r}^{HCP}$ achieves performance comparable to $IPS_{2}^{HCP}$, while the relative ranking between second-order IPSs and $IPS_{1}^{HCP}$ depends on the propagation rate $\lambda$. Not surprisingly, when $\lambda$ is set near the threshold, where the propagation process is highly sensitive to the choice of seed nodes and the accuracy of node ranking is therefore most crucial, $IPS_{2}^{HCP}$, which incorporates richer information, achieves the best performance (Fig. 4e–h). For small $\lambda$ (Fig. 4a–d, where the propagation efficiency is limited), $IPS_{2}^{HCP}$ and $IPS_{1}^{HCP}$ exhibit similar performance, implying that one-step information is already sufficient to accurately assess node influence. As for cases with extremely large $\lambda$ (Fig. 4i–l), which are rare in the real world, $IPS_{2}^{HCP}$ shows no improvement over $IPS_{1}^{HCP}$. This is because, at such extremely high spread rates, propagation largely depends on the first step. In addition, when $\mu\neq 1$, our IPS methods maintain superior performance across all parameter settings. See results in Supplementary Section S2.2.

To sum up, IPS methods consistently surpass other centrality measures while maintaining strong robustness towards different parameter regions. Its dynamics-informed formulation also provides clear physical interpretability, thereby enhancing its reliability. Moreover, our results show that although the second-order IPSs perform better in certain scenarios, particularly near the critical threshold, the first-order IPS performs comparably and remains substantially superior to existing approaches. Considering the trade-off between computational efficiency and predictive accuracy, $IPS_{1}^{P}$ is more practical for real-world applications. We therefore focus on $IPS_{1}^{P}$ in the subsequent analyses.

IPS builds on the intuition that a node’s early propagation range in a hypergraph serves as an accurate proxy for its eventual influence, thereby enabling natural generalization across diverse higher-order contagion dynamics. In this section, we further assess its transferability under a size-attenuation variant of the HCP model and the HTC model.

We consider a phenomenon where the interaction intensity and communication efficiency within groups may decline as the group size increases, particularly in face-to-face interactions St-Onge et al. (2022); Lowry et al. (2006). We model this phenomenon by introducing an attenuation term: define the infection rate for susceptible nodes within a hyperedge $h$ containing $i$ infected nodes as $\frac{\lambda}{|h|}i^{\nu}$. This model is a variant of the basic HCP model, and we call it the higher-order contagion model with size attenuation (HCSA). Following the computational methodology of $IPS^{HCP}_{1}$, we derive:

where HD$(s)$ represents the hyper-degree value of node $s$.

We evaluate the performance of $IPS_{1}^{HCSA}$ in the same way as before. Results in Fig. 5b–g show that $IPS_{1}^{HCSA}$ exhibits robust and superior performance across all scenarios, while the hyper-degree method, which is formally similar to $IPS_{1}^{HCSA}$, exhibits comparable performance. However, this arises from the “weak monotonicity” of hyper-degree, which assigns identical values to many nodes and thereby induces a spurious bias in Kendall’s $\tau$. Once the spurious gain is neutralized, $IPS^{HCSA}_{1}$ again maintains a consistent advantage over hyper-degree (see details in Supplementary Section S4.1.2). Additionally, contrary to the case of the previous HCP model (Fig. 2a), eigenvector-max performs better than eigenvector-linear in the HCSA model. This is because eigenvector-linear assigns higher centrality to nodes belonging to larger hyperedges and possessing more connections, while, eigenvector-max disregards the size of hyperedges and focuses more on the number of connections Tudisco and Higham (2021). That is, methods that are inherently more similar to the corresponding IPS approach tend to achieve better predictive performance, which is in line with the earlier observations. Above results are corresponding to the cases when $\mu=1$, we also provide more results for cases when $\mu=0.1$ in Supplementary Section S4.1.2.

We also consider the higher-order threshold contagion (HTC) model, where infection occurs collectively, i.e., if the number of infectious nodes in a hyperedge exceeds (or equals) a fraction $\theta$ of its size, all susceptible nodes in it get infected with probability $\eta$. Note that the propagation only occurs in hyperedges that satisfy the threshold condition, we have:

where $N_{v,\theta}^{1}(s)=\{i\mid i\in h\in N_{h}^{1}(s),|h|\leq 1/\theta\}$, ${k_{s}^{\theta}(i)}=|\{h\mid i\in h\in N_{h}^{1}(s),|h|\leq 1/\theta\}|$, and $N_{h,\theta}^{1}(s)=\{h\mid s\in h,|h|\leq 1/\theta\}$.

We consider $\theta\in\{0.25,\,0.5\}$, and $\eta\in\{0.1\eta_{c},\,\eta_{c},\,10\eta_{c}\,(\text{or}\,\,2\eta_{c}\,\,\text{instead if}\,\,10\eta_{c}\geq 1)\}$ to evaluate $IPS^{HTC}_{1}$. Here, $\eta_{c}$ is the threshold of $\eta$ for the given $\theta$ and recovery probability $\mu$ (see the calculation of $\eta_{c}$ in Methods). The performance of different measures is shown in Fig. 5i-n. We find that the approximations of $IPS^{HTC}_{1}$ attain superior predictive accuracy compared to all other baselines. Besides, contrary to most cases under HCP and HCSA, the projection-based 2-eigenvector and 2-degree outperform hyper-coreness-$R_{w}$ and eigenvector-linear. The performance fluctuation among baselines, together with the consistent superiority of $IPS_{1}^{HTC}$, once again underscores the importance of incorporating dynamical information into node ranking in the context of complex higher-order interactions. Above results are the cases when $\mu=1$; results when $\mu=0.1$ are in Supplementary Section S4.2.2.

From a broader perspective, the first-order IPS formulation theoretically reveals the parameter sensitivity of node influence ranking. Under the HTC model, $IPS^{HTC}_{1}$ depends explicitly on the dynamical parameters through the threshold condition, i.e, the indicator of whether $|h|\leq 1/\theta$. Consequently, $IPS^{HTC}_{1}$ at $\theta=0.25$ differs from that at $\theta=0.5$, accounting for the substantial variation in influence rankings across HTC parameter settings (Fig. 1e). By contrast, for the HCP model, $IPS^{HCP}_{1}$ is independent of the dynamical parameters, which explains the high concordance of rankings across different HCP parameter choices (Fig. 1d). Taken together, the outstanding predictive performance and clear interpretability of IPS demonstrate its strong transferability across modeling dynamics and parameter regimes.

Finally, we assess cross-paradigm transferability of IPS by moving beyond contagion to an opinion dynamics—the naming game (NG) model. The NG model simulates the emergence of convention under competitive opinion evolution, which is widely used to investigate the critical mass of committed minority opinion to overturn the majority Dall’Asta et al. (2006); Centola et al. (2018). Here, we adopt the higher-order naming game model (HNG) in the hypergraph framework Iacopini et al. (2022), and select committed minority agents based on different influence ranking algorithms to examine which approach most effectively promotes the reversal of minority opinions.

In the HNG model, there are two names (denoted as $A$ and $B$), representing different social norms. Initially, a fraction $p$ of committed nodes hold opinion vocabulary $\{A\}$, whereas the remaining (uncommitted) majority holds $\{B\}$. At each following step, a hyperedge is sampled, and two processes occur: (1) broadcast: a random speaker proposes a name from its vocabulary and shares it with other members; (2) negotiation: all members update their vocabulary (except committed nodes) based on a particular agreement rule and the intensity of social influence $\beta$. See details of the model in Methods. With the different setup of the $\beta$ and $p$, the names evolve to different abundances in steady state. Particularly, as $p$ increases, the minority can take over the system completely, where the minimal $p$ required is the critical mass $p_{c}$ Mancastroppa et al. (2023); Iacopini et al. (2022).

Following the above mechanism, we focus on its broadcast stage, which is similar to the propagation process. For a committed agent $s$, its expected broadcast range at the initial moment is naturally given as:

We evaluate the effectiveness of $IPS^{HNG}_{1}$ on reducing the critical mass $p_{c}$ in real-world hypergraphs and compare it with 11 baselines, including a random selection strategy. We select the top $\lfloor pN\rfloor$ nodes and record the steady-state abundance of name $A$ (the fraction of nodes holding only $A$, denoted as $\rho_{A}$). Figure 6 shows heatmaps of $\rho_{A}$ over the $(\beta,p)$ grid for each measure. Black regions indicate successful minority takeover ($\rho_{A}=1$), whereas colored regions correspond to $\rho_{A}<1$; for each fixed $\beta$, the smallest $p$ at which the heatmap turns black gives the critical mass $p_{c}$. Thus, a larger black region implies a lower $p_{c}$ and hence more efficient reversal. Compared with random selection, centrality-based selection substantially lowers $p_{c}$, although effectiveness varies across measures. Among them, $IPS^{HNG}_{1}$—and, as expected, the formally similar hyper-degree (HD)—outperform the alternatives by a noticeable margin. By contrast, hyper-coreness-$R_{w}$ and eigenvector-linear, which perform well under other dynamics, are less effective under HNG. Here we set the agreement rule as the unanimity rule; see results in more hypergraphs and results for the union rule in Supplementary Section S5. These results highlight the potential of the IPS method for broader applicability across multiple classes of dynamical processes.

In the higher-order naming game (HNG) model, there are two names (denoted as $A$ and $B$), and each agent (i.e., the hypergraph node) can hold one or both names in its vocabulary. For each pair of $\beta$ and $p$, we initialize the system as follows. First, we rank nodes according to a given measure and select the top $\lfloor{pN}\rfloor$ ranked nodes as the committed minority. These nodes have a fixed vocabulary $\{A\}$ and never change it. All remaining nodes start with vocabulary $\{B\}$. The subsequent name evolution proceeds in discrete time steps: At each time step, randomly select a hyperedge, and then randomly select a node in this hyperedge as the speaker. The speaker samples a name uniformly at random from its vocabulary and broadcasts this name to all other nodes in the hyperedge (hearers). The hearers in the hyperedge jointly determine whether an agreement is reached. Under the unanimity rule, agreement is achieved if all hearers already have the broadcast name in their vocabularies. Also, we can adopt a more relaxed rule, called the union rule, under which agreement is achieved if at least one hearer in the group already holds the broadcast name. If agreement is achieved, then with probability $\beta\in[0,1]$ all non-committed nodes in the hyperedge erase all other names and keep only the broadcast name. Committed nodes always keep $\{A\}$. If agreement is not achieved, then every non-committed hearer that does not already have the broadcast name adds it to its vocabulary. Then, the system goes to the next time step. The process continues until it reaches an absorbing state with $\rho_{A}=1$ (all non-committed holds the name $A$ only) or until a maximum time step $t_{\max}=1\times 10^{6}$ is reached, by which time the abundance of names has converged to a stationary regime. In the latter case, the stationary value of $\rho_{A}$ is obtained by averaging over 100 samples taken in the last $1\times 10^{5}$ time steps. For each pair of $\beta$ and $p$ in the main text, the reported $\rho_{A}$ is the median value of 100 simulations (50 simulations for the congress-bills hypergraph).

To numerically estimate the critical threshold in our SIR-type propagation models, we adopt the susceptibility-based method commonly used to detect phase transitions in finite systems Ferraz de Arruda et al. (2023); de Arruda et al. (2018); de Oliveira and Dickman (2005). The susceptibility $\chi$ is defined as $\chi=\frac{\langle R_{\infty}^{2}\rangle-\langle R_{\infty}\rangle^{2}}{\langle R_{\infty}\rangle}$, where $R_{\infty}$ is the final propagation range in one simulation (i.e., the number of recovered nodes at the end of each simulation ). This quantity measures the relative magnitude of fluctuations in outbreak size and typically peaks near the transition point, thus providing a reliable numerical estimate of the propagation threshold when an analytical solution is not available. For the HCP and HCSA model, we fix $\nu$ and $\mu$, and scan the parameter $\lambda$ on a logarithmic grid. For each $\lambda$, we perform 300 realizations when $\mu=0.1$ and 1,000 realizations when $\mu=1$, using randomly selected initial seeds. We then take the value of $\lambda$ that maximizes $\chi(\lambda)$ as the threshold $\lambda_{c}$. Similarly, for the HTC model, we fix $\theta$ and $\mu$ and scan the parameter $\eta$. The critical threshold $\eta_{c}$ is determined as the value of $\eta$ that maximizes $\chi(\eta)$ following the same procedure.

In this work, we adopt several widely used higher-order networks spanning diverse domains, comprising 20 medium-sized hypergraphs and one large-scale hypergraph (Table 1). The medium-sized hypergraphs contain topics of email communications, political affiliations, online interactions, ecological datasets, and face-to-face interactions. We follow the previous work Mancastroppa et al. (2023) to preprocess the medium-sized hypergraphs and follow the previous work Benson et al. (2018) for the large-scale hypergraph threads-math-sx.

Email communication datasets contain email-EU and email-Enron Benson et al. (2018); Paranjape et al. (2017); Klimt and Yang (2004), nodes are the email addresses within a European research and Enron corporation, respectively, while each hyperedge includes the senders and receivers of an email.

Political affiliation datasets contain congress-bills, senate-bills, house-committees, and senate-committees Porter et al. (2005); Fowler (2006); Benson et al. (2018); Chodrow et al. (2021). Nodes in these hypergraphs are members of the U.S. Congress or the U.S. Senate, and each hyperedge includes the members in a committee or cosponsoring a bill.

Online interaction datasets contain music-review, geometry-questions, and algebra-questions Ni et al. (2019); Amburg et al. (2022). Nodes in music-review are Amazon users, and each hyperedge involves nodes who review a specific product belonging to the category of blues music. Nodes in geometry-questions and algebra-questions are users in MathOverflow, and each hyperedge contains users who have answered the same question on the topic of algebra or geometry, respectively.

Ecological systems datasets contain M_PL_062_ins, M_PL_015_ins, M_PL_062_pl, and M_PL_015_pl Robertson (1977); Petanidou (1991); Fortuna et al. (2014). Nodes in M_PL_062_ins and M_PL_015_ins are insects, and each hyperedge contains insects that pollinate the same plant. By contrast, nodes in M_PL_062_pl and M_PL_015_pl are plants, and each hyperedge connects nodes pollinated by the same insect species.

Face-to-face interaction datasets are aggregated from time-resolved face-to-face human interaction databases, containing middle-school (Mid1 Toth et al. (2015)), primary-schools (LyonSchool Stehlé et al. (2011); Benson et al. (2018), Elem1 Toth et al. (2015)), high-school (Thiers Mastrandrea et al. (2015)), conference (SFHH Génois and Barrat (2018)), workplace (InVS15 Génois et al. (2015); Génois and Barrat (2018)), and hospital (LH10 Vanhems et al. (2013)). Nodes in these datasets are humans. The hyperedges are the maximal cliques of the time-resolved networks, which are aggregated over a time window of 15 minutes.

In the large-scale dataset threads-math-sx Benson et al. (2018), nodes represent the users of a forum, and each hyperedge denotes a set of users answering a question.