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

Higher-order dynamics refer to mechanisms where collective mutual or synchronous interactions differ fundamentally from their pairwise counterparts through the concept of many-body interactions. Phenomena absent in pairwise models, such as catastrophic activation, hysteresis, and hybrid transitions, emerge naturally in higher-order interacting systems. Thus, the simulation of contagion dynamics on higher-order structures is algorithmically and computationally challenging due to the complexity of propagation through hyperedges of arbitrary order. To address this issue, optimized Gillespie algorithms were constructed for higher-order structures by means of phantom processes: events that do not change the state of the system but still account for time progression. We investigate the algorithm's performance considering the susceptible-infected-susceptible (SIS) epidemic model with critical mass thresholds on hypergraphs. Optimizations were assessed on networks of different sizes and levels of heterogeneity in both connectivity and order interactions, in a high epidemic prevalence regime. Algorithms with phantom processes are shown to outperform standard approaches by several orders of magnitude in the limit of large sizes. Indeed, a high computational complexity scaling $\mathcal{O}(N^2)$ with system size $N$ of the standard algorithms is improved to low complexity scaling nearly as $\mathcal{O}(N)$. The optimized methods allow for the simulation of highly heterogeneous networks with millions of nodes within affordable computation costs, significantly surpassing the size range and order heterogeneity currently considered.