Efficient Gillespie algorithms for spreading phenomena in large and heterogeneous higher-order networks

Higher-order interactions, where groups of nodes interact collectively rather than pairwise, are central to many complex systems, from neural and ecological networks to social contagion. However, simulating dynamical processes on such higher-order structures remains computationally challenging due to the combinatorial growth of possible interactions. Here, we develop efficient and statistically exact Gillespie algorithms for Markovian spreading dynamics on large and heterogeneous hypergraphs. By incorporating phantom processes $-$events that advance time without altering the system’s state$-$, we drastically reduce the computational complexity of standard algorithms (O(N ² )), achieving near-linear scaling with system size (O(N)). Using the susceptible–infected–susceptible model with critical mass thresholds as a benchmark, we show that the optimized algorithms outperform standard approaches by several orders of magnitude, enabling simulations of networks with millions of nodes and broad heterogeneity in both degree and interaction order. These results establish a general framework for scalable, continuous-time simulations of higher-order contagion and related dynamical processes.