1 Introduction↩︎

Modeling infectious disease transmission has long been a central topic in mathematical epidemiology and complex systems science. Early compartmental models [1], [2], based on homogeneous mixing assumptions, employed systems of differential equations to describe the time evolution of classical susceptible–infected–recovered (SIR) and susceptible–infected–susceptible (SIS) dynamics [3], [4]. These frameworks laid the mathematical foundations for key epidemiological indicators such as the basic reproduction number \(R_0\) [5] and have guided public health interventions ranging from smallpox eradication to seasonal influenza control [6], [7]. However, the assumption of uniform mixing fails to capture the heterogeneous contact structures of real populations, where interaction networks shape epidemic outcomes.

The emergence of network epidemiology [8] revealed that network topology profoundly modulates disease dynamics. Short path lengths in small-world networks accelerate transmission [9], while highly connected hubs in scale-free networks can sustain infection even at infinitesimal transmission rates [10], [11]. Subsequent advances incorporated multilayer [12], [13] and metapopulation [14], [15] frameworks, enabling the modeling of spatial heterogeneity, mobility patterns, and interdependent spreading processes [16]. These developments established a structural foundation for understanding how network organization controls epidemic thresholds, prevalence levels, and intervention efficiency.

Most network-based models originally assumed fixed topologies, ignoring the feedback between individual behavior and epidemic progression. In reality, humans adapt their social interactions dynamically in response to perceived risk, information awareness, or behavioral constraints. This coevolutionary perspective gave rise to adaptive networks [17], [18], where link formation and removal depend on the evolving disease state of nodes [19]. Such feedback can induce complex dynamical phenomena absent in static structures, including bistability, hysteresis, and self-organized oscillations [20]. For example, the rewiring of high-risk contacts has been shown to delay or even suppress outbreaks [21], but excessive adaptation can also fragment the network, altering percolation and resilience properties [22]. These findings underscore that epidemic control cannot be fully understood without accounting for the joint evolution of topology and dynamics.

Network structure fundamentally determines epidemic controllability, motivating extensive studies of immunization and intervention strategies [23], [24]. At the microscopic level, individuals may spontaneously self-isolate or reduce contact frequency, while organized interventions at the macroscopic level (e.g., targeted vaccination or contact restrictions) remove nodes or edges to interrupt transmission pathways. [25]. The fraction of immunized nodes or edges required to eliminate endemic persistence defines the Herd Immunity Threshold (HIT) [26]. Due to the inefficiency of random immunization in heterogeneous networks, structural information has been leveraged to develop targeted methods, including degree-based vaccination [27], acquaintance immunization [28], and influence-based strategies [29]. Beyond node-level interventions, recent studies have demonstrated that selectively removing or immunizing edges based on network topology can substantially enhance containment efficiency. For instance, Liu et al. [30] proposed an explosive percolation-based edge removal strategy that selectively disconnects low-degree contacts to curb misinformation spread, while a subsequent work [31] combined greedy optimization with topological features to develop computationally scalable edge immunization schemes, showing that jointly optimizing local connectivity patterns and global structure yields superior performance. Recent work further emphasizes group-level targeting in realistic contact contexts such as households, schools, and workplaces, highlighting the need for intervention models that operate beyond pairwise links [32]. In parallel, the microscopic Markov chain approach has been applied to analyze multilayer behavioral contagion [33].

Traditional network representations capture only pairwise interactions, yet real contagion processes often unfold within groups, such as meetings, classrooms, or social gatherings, where multiple individuals interact simultaneously [34], [35]. Hypergraphs and simplicial complexes provide natural mathematical formalisms for such higher-order structures [36]. In the \(d\)-uniform hypergraph, each hyperedge connects exactly \(d\) nodes, with \(d=2\) reducing to a standard network. Recent studies employing the simplicial SIS (\(s\)-SIS) model [37], [38] have demonstrated that group interactions introduce nonlinear reinforcement effects. These effects fundamentally alter epidemic dynamics, giving rise to discontinuous phase transitions, bistability, and hysteresis phenomena [39], [40]. Moreover, higher-order coupling lowers effective epidemic thresholds [41], [42], implying that neglecting group effects can severely underestimate outbreak risks. Beyond static higher-order descriptions, recent studies have extended simplicial epidemic models to temporal networks and spatial settings, incorporating control strategies based on critical simplex identification [43] and analyzing pattern dynamics through reaction–diffusion frameworks on higher-order topologies [44], [45]. However, real-world group interactions are neither static nor continuous, as they vary in frequency and intensity, often adapting to behavioral or environmental feedback.

Despite progress in adaptive networks and higher-order epidemic modeling, a unified framework integrating adaptive group interactions and strategic interventions remains limited. Most existing works treat adaptation as either individual-level rewiring or static group suppression, and even those that consider hyperedge removal typically rely on fixed structural rankings rather than dynamically updated risk measures. The intervention therefore does not respond to the shifting landscape of transmission risk during an ongoing outbreak. To bridge this gap, we develop an adaptive \(s\)-SIS model on \(d\)-uniform hypergraphs, where the evolution of node states and hyperedge activity are mutually coupled. Within this framework, we introduce two classes of hyperedge-level interventions: (i) Risk-based immunization, including spontaneous isolation and external deactivation; (ii) Structural rewiring, implemented through random or degree-preferential strategies that reshape the topology of the network. We extend the microscopic Markov chain approximation (MMCA) [46] to analyze this model, deriving analytical conditions for stationary states and nonlinear epidemic thresholds. The theoretical predictions are further validated by Monte Carlo (MC) simulations, allowing us to quantify how the interplay of topology, behavior, and interventions governs epidemic resilience.

To summarize, the principal contributions of our work are as follows.

1. We propose and formalize an adaptive \(s\)-SIS model on uniform hypergraphs, establishing a unified framework for higher-order epidemic dynamics with feedback-driven hyperedge activity.

2. We design and systematically compare multiple hyperedge-level immunization and rewiring strategies, quantifying their effects on the stationary-state infection density and network structure.

3. We demonstrate the emergence of discontinuous phase transitions, nonlinear thresholds, and hysteresis phenomena through analytical derivations and extensive simulations, revealing how adaptive feedback and initial conditions jointly shape epidemic outcomes.

The remainder of this paper is structured as follows. Section 2 develops the adaptive simplicial SIS model, integrating both intervention strategies and rewiring mechanisms within a unified adaptive framework, and presents the corresponding MMCA analytical formulation. Section 3 conducts numerical simulations and Monte Carlo validations. Finally, in Section 4, we summarize the conclusions and discuss their broader implications.

2 Modelling↩︎

This section introduces the adaptive simplicial SIS model on \(d\)-uniform hypergraphs and its associated hyperedge-level intervention mechanisms. We first describe the adaptive structure that couples node-level epidemic dynamics with hyperedge activity in Section 2.1, followed by the formulation of two classes of intervention strategies, immunization and rewiring, which act at the group level to mitigate epidemic spread in Section 2.2. Finally, in Section 2.3, we provide a theoretical analysis of the model using a homogeneous MMCA framework to derive nonlinear thresholds and equilibrium conditions. Table 1 summarizes the notations used throughout this paper.

2.1 Simplicial SIS model with adaptive hyperedges↩︎

We establish an adaptive epidemic modelling framework that explicitly integrates the coevolution between disease transmission and higher-order social interactions. In many real-world systems, such as group gatherings, collaborative teams, or biological complexes, contagion does not spread solely through pairwise contacts but through group-based interactions involving multiple individuals simultaneously. To capture these higher-order effects, we represent the population as a \(d\)-uniform hypergraph, which naturally extends classical pairwise networks by allowing each hyperedge to connect \(d\) distinct nodes instead of two.

Formally, we generalize the Erdős–Rényi random graph \(G(n,m)\) to its higher-order analogue \(H(n,m,d)\), where \(n\) denotes the number of nodes (individuals), \(m\) the number of hyperedges (group interactions), and \(d\) the hyperedge dimension representing group size. The vertex set is defined as \(\mathcal{V}=\{v_1,v_2,\dots,v_n\}\), and the hyperedge set as \(\mathcal{E}=\{e_1,e_2,\dots,e_m\}\), where each hyperedge \(e\subseteq\mathcal{V}\) satisfies \(|e|=d\). In this construction, all hyperedges are treated as statistically equivalent, implying uniform group sizes and random composition of members.

In the sparse regime where \(m\ll n^d\), the probability that two nodes share multiple hyperedges becomes negligible, and the hyperdegree distribution referred to the number of hyperedges incident to a node is approximately Poisson with mean \[\langle k \rangle = \frac{md}{n} = \frac{1}{n}\sum_{i=1}^{n} |\mathcal{E}_i|,\] where \(\mathcal{E}_i\) denotes the set of hyperedges containing node \(i\). The average hyperdegree \(\langle k\rangle\) characterizes the expected number of group interactions per individual and thus determines the effective connectivity of the system.

In real-world social and biological systems, contagion processes are rarely independent of the network structure on which they occur. Individuals tend to alter their group participation in response to perceived infection risk, while the collective organization of social groups evolves as a consequence of behavioral adaptation. Such adaptive feedback between disease dynamics and network structure has been widely observed during epidemics, where people spontaneously reduce high-risk interactions but gradually restore connectivity once perceived danger declines. To capture this adaptive interplay, our model introduces a dynamical mechanism in which epidemic dynamics and hyperedge activity evolve simultaneously.

Specifically, the model integrates two coupled processes, as illustrated in Fig. 1 (a):

1. Node-level epidemic dynamics: Each node \(i\) exhibits a binary state \(X_i(t) \in \{0,1\}\), representing susceptible (\(S\)) and infected (\(I\)) states. Infection transmission occurs through higher-order interactions, requiring a critical configuration in which all other \(d-1\) members within a hyperedge are infected. Each such hyperedge independently transmits infection with probability \(\beta_d\) per time step, while infected nodes recover with probability \(\mu\).

2. Hyperedge-level adaptive dynamics: Each hyperedge \(e\) is endowed with a dynamic activity level \(\Theta_e(t) \in [0,1]\), representing its effective interaction intensity. This variable captures behavioral adaptation at the group level: hyperedges with high infection pressure tend to reduce activity, while inactive ones gradually regain it over time. The temporal evolution of \(\Theta_e(t)\) follows \[\label{eq:Theta95update} \Theta_e(t+1) = \Theta_e(t)\left[1-\eta\Phi_e(t)\right] + \gamma\left[1-\Theta_e(t)\right],\tag{1}\] where \(\eta \in [0,1]\) quantifies sensitivity to local infection pressure, and \(\gamma \in [0,1]\) characterizes spontaneous recovery of interaction activity. The first term drives the activity downward when the hyperedge is exposed to high infection pressure \(\Phi_e(t)\), while the second term drives a spontaneous return toward full activity irrespective of the current state. Together they produce a dynamic balance between risk-driven suppression and intrinsic recovery.

The hyperedge infection pressure \(\Phi_e(t)\), defined as the joint probability that all nodes in \(e\) are simultaneously infected, is given by \[\label{eq:Phi} \Phi_e(t) = \Pr ( \bigcap_{i \in e} \{X_i(t) = I\}) \approx \prod_{i \in e} p_i(t) ,\tag{2}\] where \(p_i(t) = \Pr(X_i(t) = I)\). Equation (3) extends the edge epidemic importance (EI)-based immunization strategy introduced in [47]. It computes the joint probability that all nodes within a hyperedge are simultaneously infected, which serves as the basis for ranking and deactivating or rewiring hyperedges in our intervention strategies. This probability follows from the individual-based mean-field (IBMF) theory on hypergraphs, which tracks each node’s infection probability \(p_i\) independently and neglects dynamical correlations between node states [46]. Hence the approximation \(\Phi_e(t) \approx \prod_{i\in e} p_i(t)\) is a natural consequence of the IBMF closure.

The probability that susceptible node \(i\) avoids infection through all incident hyperedges is \[\label{eq:qi40t41} q_i(t) = \prod_{e\in\mathcal{E}_i}\left(1-\Theta_e(t)\beta_d\prod_{\substack{j\in e \\ j\neq i}} p_j(t)\right),\tag{3}\] leading to the complete node-state evolution as \[\label{eq:mmca95node} p_i(t+1) = \left(1-p_i(t)\right)\left(1-q_i(t)\right) + (1-\mu)p_i(t).\tag{4}\]

Eq. 1 –4 define the microscopic dynamics. They capture the core feedback loop where increased local infection pressure reduces hyperedge activity introduced in Eq. 1 , which in turn lowers transmission probabilities in subsequent steps in Eq. 3 , and vice versa.

For analytical tractability, the model relies on several simplifying assumptions: (i) the hypergraph is \(d\)-uniform and constructed uniformly at random, so that hyperedges are statistically equivalent and the hyperdegree distribution is approximately Poisson; (ii) all hyperedges share the same adaptation parameters \(\eta\) and \(\gamma\), implying uniform responsiveness to infection pressure across groups; (iii) node states within a hyperedge are treated as independent.

2.2 Hyperedge-level intervention strategies↩︎

In realistic social systems, both individual states and group interactions evolve dynamically in response to epidemic risk. To capture these adaptive processes, we integrate intervention strategies directly into the model framework, allowing the hypergraph structure to reorganize under infection pressure. Two complementary classes of interventions are considered: (i) Immunization, which reduces transmission by deactivating highly infectious hyperedges illustrated in Fig. 1 (b), and (ii) Rewiring, which reconstructs inactive hyperedges to preserve global connectivity shown in Fig. 1 (c). These mechanisms jointly define the adaptive feedback governing higher-order epidemic evolution.

2.2.1 Immunization strategies↩︎

Immunization aims to suppress contagion by selectively deactivating risky hyperedges, thereby reducing the effective transmission channels without completely dismantling the network. To capture both decentralized behavioral responses and centralized policy-driven interventions, we design two complementary strategies and describe their key principles and corresponding mathematical formulations below.

(1) Spontaneous Isolation (SI). This mechanism mimics decentralized self-protection behavior commonly observed in real epidemics, where social groups spontaneously reduce interaction when perceiving elevated infection risk. As illustrated in Algorithm 2, each hyperedge \(e\) has an activity level \(\Theta_e(t)\) that reflects its current interaction intensity. At every discrete time step, the system evaluates all hyperedges and identifies those whose activity falls below a critical threshold \(\theta_{\mathrm{min}}\) as \[\mathcal{E}_{\text{isolate}}(t) = \left\{ e \in \mathcal{E} : \Theta_e(t) < \theta_{\mathrm{min}} \right\}. \label{eq:isolate95set}\tag{5}\] These hyperedges are considered likely to suspend contact spontaneously.

However, to model the realistic limitation of resources and compliance, the total number of isolated hyperedges cannot exceed a global immunization budget \(w \in [0,1]\), corresponding to \(m_{\text{max}} = \lfloor w \cdot m \rfloor\) eligible hyperedges. Operationally, the budget \(w\) reflects the fraction of hyperedges that can be simultaneously targeted under limited public health capacity, such as the number of available response teams or isolation facilities, and the same uniform budget constraint applies to every hyperedge. Beyond constraining the intervention scale, \(w\) also serves as the control parameter for identifying the HIT.

Algorithmically, all candidate hyperedges in \(\mathcal{E}_{\text{isolate}}(t)\) are ranked by their current activity level \(\Theta_e(t)\) in descending order, and the top \(m_{\text{max}}\) are selected for isolation. Formally, for each isolated hyperedge \(e\), its activity and effective transmission capacity are set to zero, expressed as \[\Theta_e(t+1) = 0, \qquad \text{active}[e] \leftarrow \text{False},\] ensuring that deactivated hyperedges no longer participate in subsequent infection or rewiring events, effectively removing them from the propagation process. Such local, threshold-based decision rules capture spontaneous behavioral adaptation without requiring any global information about network state.

(2) Targeted Immunization (TI). Unlike SI, the targeted immunization scheme represents a centralized, information-driven strategy designed to allocate limited intervention resources most efficiently. Its logic parallels optimal control approaches, where the system identifies and deactivates hyperedges that contribute most to sustaining infection. As illustrated in Algorithm 3, we employ the hyperedge infection pressure \(\Phi_e(t)\) introduced in Eq. (2 ), defined as the joint probability that all nodes within \(e\) are infected by \[\Phi_e(t) = \prod_{i\in e} p_i(t),\] where the approximation follows from statistical independence under the MMCA framework. This criterion extends the edge epidemic importance measure from network-based immunization to higher-order systems. Unlike previous hyperedge removal strategies that rank group interactions by static structural properties, the infection pressure \(\Phi_e(t)\) is recomputed at every time step as the node infection probabilities evolve, so the set of targeted hyperedges adapts to the real-time epidemic situation rather than being fixed a priori. Intuitively, \(\Phi_e(t)\) serves as a quantitative measure of the overall infection burden within a hyperedge, establishing itself as a principled risk metric for group-level transmission.

The algorithm ranks all active hyperedges in decreasing order of their infection pressure \(\Phi_e(t)\) and selects the top \(m_{\text{target}} = \lfloor w \cdot m \rfloor\) hyperedges for immunization. Formally, let \(\varepsilon \subseteq \mathcal{E}\) denote a candidate subset of hyperedges with cardinality \(|\varepsilon| = m_{\text{target}}\). The optimal target set is determined by maximizing the cumulative infection pressure over all such subsets, \[\mathcal{E}_{\text{target}} = \arg\max_{\varepsilon \subseteq \mathcal{E}} \sum_{e \in \varepsilon} \Phi_e(t).\] Each selected hyperedge \(e \in \mathcal{E}_{\text{target}}\) is subsequently deactivated by setting its activity \(\Theta_e(t+1)=0\) and flagging it as inactive, consistent with the procedure used in spontaneous isolation. This approach concentrates intervention effort on the most influential hyperedges, which act as key transmission hubs, thus achieving a larger reduction in global infection density under the same resource constraints.

2.2.2 Rewiring mechanism↩︎

While both SI and TI suppress transmission by removing high-risk connections, such removal inevitably alters the network’s connectivity pattern and may lead to fragmentation. In real social systems, however, individuals rarely remain permanently disconnected. Instead, they form new groups or reorganize their interactions after previous connections are disrupted. To capture this inherent adaptive behavior, we introduce a rewiring process that regenerates new hyperedges in place of the immunized ones. This mechanism maintains the total number of hyperedges \(m\) constant, ensuring the system preserves its size while allowing the higher-order structure to evolve dynamically. Hence, immunization and rewiring jointly constitute a feedback loop—one that suppresses infection through selective deactivation while simultaneously restoring connectivity through adaptive reorganization.

(1) Random Rewiring. This mechanism represents unbiased social mixing, where each new hyperedge is sampled uniformly at random from the set of all possible \(d\)-node combinations of \(\mathcal{V}\): \[P\left \{e_{\mathrm{new}} = e\right \} = \frac{1}{|\mathcal{U}_d|}, \qquad \forall\, e \in \mathcal{U}_d,\] where \(\mathcal{U}_d\) denotes the collection of all distinct \(d\)-element subsets of the node set. To model cautious behavioral recovery, all new hyperedges are initialized with a small baseline activity \(\Theta_0 = 0.1\), reflecting the gradual restoration of social contacts after isolation.

(2) Degree-Preferential Rewiring. In many social systems, individuals with higher connectivity are more likely to form new groups after disruptions. To incorporate this empirically observed bias, we define the probability that node \(i\) is selected during rewiring as proportional to its current hyperdegree \(k_i\), adjusted by a small offset \(\alpha>0\) to prevent isolation of low-degree nodes as \[P_{\text{pref}}(i) = \frac{k_i + \alpha}{\sum_{j=1}^n (k_j + \alpha)} = \frac{k_i + \alpha}{md + n\alpha}.\] Newly generated hyperedges under this rule are also initialized with \(\Theta_0=0.1\), ensuring consistency across different rewiring strategies.

It is worth noting that the proposed intervention strategies differ substantially in their information requirements and operational practicality. Spontaneous isolation (SI) relies only on local threshold judgments, since each hyperedge simply compares its own activity \(\Theta_e\) against \(\theta_{\min}\), and can therefore be implemented in a decentralized manner with minimal information cost. Targeted immunization (TI), by contrast, requires global sorting of all active hyperedges by their infection pressure \(\Phi_e\), demanding centralized surveillance and coordination. Similarly, random rewiring requires no network information, whereas degree-preferential rewiring needs the current hyperdegrees of all nodes. These differences imply a trade-off between efficacy and practicality, since strategies that use more information tend to perform better but incur higher implementation costs.

Together, these two intervention components constitute a closed adaptive control loop linking microscopic behavioral responses to macroscopic network resilience. The former dynamically suppresses risky interactions, while the latter continuously reorganizes higher-order connectivity, yielding a flexible framework for analyzing the coevolution of epidemic spreading and social adaptation. In this framework, each new hyperedge is sampled uniformly from the set of all distinct \(d\)-node subsets, structurally excluding self-loops and repeated hyperedges. Moreover, when a hyperedge is deactivated, either through adaptive feedback or external intervention, its member nodes remain fully active in all other hyperedges they belong to. The deactivation only removes the transmission pathway associated with that particular group, reflecting a targeted, partial disruption of social contacts rather than a global isolation of the involved individuals.

2.3 Theoretical analysis↩︎

To gain analytical insight into the adaptive dynamics of the simplicial SIS model, we perform a homogeneous reduction of the MMCA equations, treating nodes and hyperedges as statistically equivalent. This scalar approximation replaces microscopic variables with global averages, yielding closed-form self-consistency relations for the mean infection density and average hyperedge activity, which are convenient for threshold and bifurcation analysis.

Proof. Let \(\rho(t)\) denote the average infection probability of nodes at time \(t\) and \(\Theta(t)\) the mean hyperedge activity: \[\rho(t) = \frac{1}{n}\sum_{i=1}^{n} p_i(t), \qquad \Theta(t) = \frac{1}{m}\sum_{e \in \mathcal{E}} \Theta_e(t).\]

Under the homogeneous approximation, each node participates in \(\langle k \rangle\) hyperedges, and the average probability that a susceptible node avoids infection from all its incident hyperedges according to Eq. 3 becomes \[q(t) = \left [ 1 - \Theta(t)\, \beta\, \rho(t)^{\,d-1} \right ]^{\langle k \rangle},\] where \(\rho(t)^{d-1}\) approximates the probability that all other \(d-1\) members of a hyperedge are infected. The corresponding evolution of the mean infection density is \[\rho(t+1) = (1-\rho(t))\big(1-q(t)\big) + (1-\mu)\rho(t). \label{eq:rho95t}\tag{6}\]

Averaging the hyperedge activity dynamics gives \[\Theta(t+1) = \Theta(t)\big(1 - \eta \rho(t)^d\big) + \gamma \big(1-\Theta(t)\big),\] with \(\rho(t)^d\) approximating the average infection pressure. At equilibrium where \(\rho(t+1)=\rho(t)=\rho\), \(\Theta(t+1)=\Theta(t)=\Theta^{\star}\), the stationary hyperedge activity is \[\Theta^{\star} = \frac{\gamma}{\gamma + \eta \rho^d}.\]

Substituting \(\Theta^{\star}\) into Eq. 6 , we set up a self-consistent equation for \(\rho\) in the stationary state as \[\rho = \frac{1-\rho}{\mu} \left[ 1 - \left( 1 - \frac{\gamma \beta \rho^{d-1}}{\gamma + \eta \rho^d} \right)^{\langle k \rangle} \right]. \label{eq:self}\tag{7}\]

Using the binomial expansion \((1-\epsilon)^n = 1 - n\epsilon + \tfrac{n(n-1)}{2}\epsilon^2 + O(\epsilon^3)\) to expand near \(\rho\to 0\), we define the self-consistency function \(G(\rho)\) as \[\label{eq:G} G(\rho) = \frac{\langle k \rangle}{\mu} \beta \, \rho^{d-1} + o(\rho^{d-1}),\tag{8}\] since \(\Theta^{\star}\to 1\) as \(\rho\to 0\).

For \(d\ge 3\), the leading-order term in Eq. 8 scales as \(\rho^{d-1}\) with \(d-1\ge 2\), so the linear term vanishes and \(\rho=0\) is stable against infinitesimal perturbations. Non-zero equilibria therefore arise only via a saddle-node bifurcation at finite \(\beta\). Thus, for \(d \geq 3\), the system exhibits bistability and abrupt first-order transitions at \(\beta_n\). In this case, the nonlinear threshold \(\beta_n\) is identified by rewriting the self-consistent equation of Eq. 7 as a parametric function: \[\beta_n = \frac{1 - \left( 1 - \frac{\mu \rho}{1-\rho} \right)^{1/\langle k\rangle}}{\Theta^{\star}\rho^{d-1}}.\] where \(\Theta^{\star} = \gamma / (\gamma + \eta \rho^d)\) and \(\beta_n\) corresponds to the smallest value for which endemic equilibria exist, marking the onset of endemic persistence. \(\square\)

Corollary 1. For pairwise spreading dynamics \((d=2)\), we get the classical epidemic threshold \(\beta_c\) with the linear form \[\beta_c = \frac{\mu}{\Lambda_{\max}(\mathbf{M})},\] where \(\Lambda_{\max}(\mathbf{M})\) denotes the largest eigenvalue of the adjacency matrix of the underlying network.

Proof. When \(d=2\), the self-consistent mapping \(G(\rho)\) in Eq. 8 expands linearly around \(\rho=0\): \[G(\rho) = \frac{\langle k\rangle}{\mu}\,\beta\, \rho + o(\rho).\] The disease-free equilibrium \(\rho=0\) is locally stable if small perturbations decay, which in a one-dimensional map requires \(G'(0) < 1\). Substituting the linear expansion gives \[\frac{\langle k\rangle}{\mu}\,\beta < 1,\] leading to the classical threshold \[\beta_c = \frac{\mu}{\langle k\rangle}. \label{eq:beta95c95linear}\tag{9}\]

Under the homogeneous-degree approximation, we have \(\Lambda_{\max}(\mathbf{M}) \approx \langle k\rangle\), and thus Eq. 9 may be equivalently expressed as \[\beta_c = \frac{\mu}{\Lambda_{\max}(\mathbf{M})}.\] \(\square\)

Remark 1. The distinction between linear and nonlinear epidemic thresholds stems from fundamental differences in early-stage dynamics: pairwise interactions (\(d=2\)) permit linear approximation and spectral characterization, whereas genuine higher-order interactions (\(d\geq 3\)) eliminate linear instability at the disease-free state due to the \(\rho^{d-1}\) scaling. This mathematical foundation explains the emergence of discontinuous transitions and bistability in higher-order contagion processes—phenomena absent in classical pairwise models. Collectively, these results establish that adaptive hyperedge activity fundamentally transforms classical epidemic thresholds and induces nonlinear, bistable regimes in systems with higher-order interactions, providing a theoretical basis for understanding abrupt epidemic outbreaks and hysteresis effects in complex social systems.

3 Numerical simulations↩︎

In this section, we conduct comprehensive numerical simulations of the adaptive dynamics between higher-order network topology and epidemic spreading, building upon the theoretical framework established above. Our experiments systematically examine how the topology-dynamics coupling shapes epidemic outcomes across four complementary dimensions: 1) verification of discontinuous phase transitions; 2) characterization of adaptive hyperedge dynamics; 3) evaluation of immunization strategies; 4) assessment of rewiring mechanisms; and 5) validation on empirical higher-order collaboration networks. In the \(s\)-SIS model the fully susceptible state is absorbing, therefore finite-size simulations tend to prematurely freeze in the disease-free trap. To circumvent this and obtain reliable stationary averages, we adopt the quasistationary (QS) method [48]. Following the original recommendation, and confirmed by our own preliminary runs, the system is kept active by storing a history of 50 previously visited configurations; upon absorption, the current state is replaced with a randomly drawn historical one with probability \(0.2\).

Unless otherwise specified, we set the network size with \(n = 1500\) nodes with hyperedge size \(d = 3\). The epidemic dynamics are parameterized with infection rate \(\beta = 0.2\), average hyperdegree \(\langle k \rangle = 6\), recovery rate \(\mu = 0.1\), while the adaptive hyperedge dynamics are governed by the sensitivity parameter \(\eta = 0.6\), which captures how promptly groups reduce activity upon perceiving infection risk, and the spontaneous recovery rate \(\gamma = 0.02\). Simulations are initialized with an initial infection density \(I(0) = 0.3\). The offset parameter for degree-preferential rewiring is set to \(\alpha = 1.0\). All results obtained through theoretical or MC simulations represent averages over 200 independent simulations.

3.1 Verification of discontinuous phase transitions↩︎

To validate the effectiveness of the proposed MMCA formulation in reproducing epidemic dynamics on adaptive hypergraphs, we conduct a comparative analysis between theoretical predictions and MC simulations. Specifically, we investigate the stationary-state infection density \(\rho\) as a function of the infection rate \(\beta\) under three different average hyperdegrees \(\langle k \rangle = 6, 9, 12\). The results are displayed in Fig. 4. The experimental setup enables us to evaluate the theoretical epidemic threshold and verify the existence of a first-order phase transition predicted by the MMCA analysis.

As illustrated in Fig. 4, the theoretical curves and simulation results exhibit strong agreement across all parameter regimes, confirming the validity of the MMCA in capturing macroscopic epidemic dynamics on hypergraphs. In Fig. 4 (a), both theoretical and numerical results demonstrate a clear discontinuous transition: when \(\beta\) exceeds a critical value \(\beta_n\), the system abruptly shifts from a disease-free to an endemic state. This behavior is fully consistent with our theoretical analysis in Theorem 1, which establishes that higher-order interactions induce a first-order (discontinuous) phase transition and a genuinely nonlinear epidemic threshold arising from a saddle-node bifurcation. Moreover, the stationary infection density increases with the average hyperdegree \(\langle k \rangle\), indicating that enhanced group connectivity facilitates persistent infection through multiple higher-order transmission pathways.

To further elucidate the temporal evolution near the critical region, Figs. 4 (b) and (c) represent the time evolution of each compartmental fraction for \(\langle k \rangle=12\) at \(\beta=0.08\) and \(\beta=0.10\), respectively. For \(\beta=0.08 < \beta_n\), the infection gradually vanishes, and the system converges to a disease-free equilibrium after mild transient oscillations. In contrast, when \(\beta=0.10 > \beta_n\), the infection persists and stabilizes at a finite endemic level after pronounced damped oscillations. These oscillatory patterns originate from the intrinsic feedback loop between epidemic spreading and adaptive topology. An increase in infection density enhances the effective infection pressure \(\Phi_e(t)\) across active hyperedges, which in turn suppresses their activity levels \(\Theta_e(t)\) through the adaptive rule in Eq. 1 . The ensuing reduction in hyperedge activity lowers the overall transmission capacity, leading to a subsequent decline in infection prevalence. As the epidemic subsides, hyperedges gradually regain activity via spontaneous recovery, thereby reopening transmission channels and initiating the next fluctuation cycle. The interplay between node-level infection and group-level adaptation produces self-regulated oscillations around the endemic equilibrium and gives rise to the hysteresis phenomena observed in adaptive higher-order contagion systems. Such behavior highlights that epidemic persistence in hypergraphs is not solely determined by transmission rates but also by the dynamical coupling between contagion and structural adaptation, which jointly shape the resilience of networks and transition pathways.

3.2 Characterization of adaptive hyperedge dynamics↩︎

To systematically examine how adaptive hyperedge behavior shapes epidemic outcomes in hypergraphs, we analyze the joint effects of initial conditions and adaptive parameters on outbreak dynamics and stationary-state prevalence through time-series comparisons and heatmaps in Fig. 5.

Figs. 5 (a) and (d) depict the temporal evolution of infection density \(\rho^I(t)\) under varying initial conditions and adaptive parameters. Interestingly, as shown in panel (a), for high sensitivity \(\eta = 0.8\), a large initial infection density \(I(0) = 0.8\) leads to disease extinction, whereas a moderate initial prevalence \(I(0) = 0.2\) results in endemic persistence. This counter-intuitive behavior manifests as bistability and hysteresis phenomena, where the system exhibits two stable stationary states, disease-free and endemic, for identical parameter values, with the final outcome determined solely by initial conditions. The bistable regime arises from the adaptive feedback between node states and hyperedge activity. Specifically, high initial infection density induces strong infection pressure \(\Phi_e(t)\) on hyperedges, which under high sensitivity rapidly suppresses hyperedge activity \(\Theta_e(t)\), ultimately extinguishing the outbreak through self-induced topological isolation. Moreover, transient oscillations arise from delayed feedback between infection growth and structural adaptation, producing damped oscillations that highlight the rich phenomenology of adaptive dynamics in higher-order networks.

Furthermore, the related parameter analysis in Fig. 5 (b)–(f) quantifies these adaptive mechanisms. Heatmaps (b) and (c) reveal that increasing the sensitivity parameter \(\eta\) and initial infection density \(I(0)\) systematically suppresses endemic prevalence, effectively expanding the disease-free region in the heatmaps. Conversely, heatmaps (e) and (f) demonstrate that a higher spontaneous recovery rate \(\gamma\) leads to increased epidemic prevalence by accelerating hyperedge reactivation. These opposing effects reflect the dual role of adaptive dynamics, as strong infection-induced suppression with high sensitivity \(\eta\) provides natural containment, while rapid structural recovery with \(\gamma\) maintains transmission potential.

3.3 Evaluation of immunization strategies↩︎

Next, to evaluate the performance of hyperedge-level immunization strategies in adaptive hypergraphs, we compare the stationary-state infection density \(\rho^I\) as a function of the immunization rate \(w\) across structural settings \(\langle k \rangle\) and adaptive parameters \(\gamma\). The results are displayed in Fig. 6. As shown in Fig. 6 (a)–(b), all strategies perform similarly to random immunization at small \(w\), but once \(w\) exceeds a critical threshold, \(\rho^I\) drops sharply, reaching herd immunity at a lower \(w_c\) than the random baseline.

In dense networks with faster recovery shown in Figs. 6 (b) and (d), SI₂ and SI₃ outperform other strategies under limited budgets (\(w < 0.6\)). This effect occurs because these strategies preferentially remove highly active hyperedges, which have relatively low instantaneous infection pressure but play an important role as potential transmission pathways in highly connected systems. Specifically, the higher threshold of SI₃ allows it to quickly block epidemiologically active hyperedges at low budgets, explaining its early advantage. As the budget grows, however, the expansion of immunization to less active hyperedges is delayed by SI₃’s descending activity-based ranking, which postpones intervention on those hyperedges with the highest infection pressure but severely suppressed activity. These residual hyperedges sustain endemic circulation and prevent full elimination. In contrast, SI₁ targets such high-risk hyperedges from the start, achieving lower stationary prevalence at sufficiently large budgets. When \(w\) increases further, TI and SI₁ become more effective since they directly target hyperedges with high infection pressure or lower activity levels, respectively. SI₂ and SI₃ continue to prioritize hyperedges with high activity but low infection pressure, leading to diminishing returns in the later stage. These results reveal a key trade-off in adaptive hypergraphs: proactive isolation of potential transmission hubs enhances resilience under resource constraints, whereas targeted interventions guided by infection pressure are crucial for complete epidemic elimination.

Then, we examine how different immunization strategies influence the herd immunity threshold (HIT) in adaptive hypergraphs. This set of simulations quantifies the way in which transmission intensity, initial epidemic prevalence, and behavioral thresholds jointly determine the minimal intervention effort required to halt an outbreak in higher-order adaptive networks. As shown in Fig. 7, the HIT \(w_c\) increases monotonically with the infection rate \(\beta\) for every strategy, confirming that stronger transmission demands a larger fraction of immunized hyperedges to prevent endemic persistence. More importantly, the performance gap separating the strategic interventions (TI and SI₁) from random immunization widens substantially at higher transmission rates, for instance when \(\beta\) exceeds \(0.31\). This divergence underscores the rapidly diminishing returns of non-strategic approaches once the epidemic becomes more difficult to contain.

Among all methods examined, SI₁ consistently attains the lowest HIT. Its effectiveness stems from its ability to identify and disable hyperedges whose activity has been severely suppressed by high local infection pressure, precisely those groups that act as persistent reservoirs of infection. In Fig. 7 (b), TI delivers a comparable performance, especially at moderate infection rates, while SI₂ and SI₃ require substantially higher thresholds. The weaker performance of these latter strategies reflects their slower responsiveness to early signs of activity suppression, which delays the removal of genuinely high-risk hyperedges. The comparison between Figs. 7 (a) and (b) further reveals that raising the initial infection level from \(I(0)=0.3\) to \(I(0)=0.5\) elevates the HIT across all strategies, yet the relative advantage of SI₁ is preserved. This upward shift occurs because a higher initial prevalence triggers stronger endogenous reductions in hyperedge activity, allowing spontaneous isolation to work in concert with the built-in adaptive feedback. The resulting cooperative coupling amplifies the impact of isolating high-risk hyperedges, delivering improved containment efficiency even when the system starts from a considerably more adverse epidemiological state. This finding reveals that herd immunity in adaptive hypergraphs is not a static network property but an emergent outcome of interacting epidemic, behavioral, and intervention dynamics, and that effective control therefore demands strategies which couple transmission dynamics with adaptive structural interventions.

3.4 Assessment of rewiring mechanisms↩︎

To investigate whether structurally guided resumption of group interactions can sustain the control achieved by prior immunization, we apply two classic rewiring mechanisms, random rewiring and degree-preferential rewiring, after the system has reached a quasi-stationary state under targeted immunization. These rewiring rules model the reopening of social contacts and the re‑establishment of group activities in a post‑outbreak recovery phase, allowing us to ask whether an unstructured return to normal connectivity leads to epidemic resurgence and whether a structured recovery can keep prevalence at a lower level. The assessment is based on both the temporal evolution of infection density and the changes in degree distribution, and the results are presented in Fig. 8.

As illustrated in Fig. 8, the two rewiring mechanisms produce qualitatively different epidemic trajectories. Under random rewiring, shown in panels (a) and (b), the intervention induces only a temporary reduction in infection density. After this brief suppression, the epidemic rapidly returns to its pre‑intervention level. Random reconnection restores overall connectivity without structurally blocking transmission pathways, so alternative spreading routes re‑emerge and hyperedge activity gradually recovers, allowing the infection to re‑establish itself.

In contrast, degree‑preferential rewiring, shown in panels (c) and (d), yields a lasting reduction in infection density. The system remains in a controlled regime well below both the baseline level and the level observed under random rewiring. The degree distribution provides a structural explanation for this persistence. Random rewiring keeps the network close to a Poisson‑like pattern, whereas degree‑preferential rewiring creates a heavy‑tailed distribution that reshapes the transmission landscape. The heavy-tailed degree distribution stems directly from the intrinsic bias of the preferential rewiring kernel, which amplifies existing degree heterogeneity regardless of the infection state. This heterogeneous structure generates high infection pressure on the few most connected nodes, which in turn triggers stronger adaptive suppression of their incident hyperedges and drives the system toward lower endemic prevalence. Hence, even in the recovery phase following an outbreak, strategic organisation of contacts and activities remains essential for preventing epidemic resurgence.

3.5 Validation on empirical higher-order collaboration networks↩︎

To assess the applicability of the proposed model to real-world settings with heterogeneous group sizes, we further test it on an empirical hypergraph constructed from congressional bill cosponsorship data [49]. The dataset contains 1218 nodes, each representing a US congressperson, and 8736 hyperedges corresponding to bills jointly cosponsored in a given year. The hyperedges vary substantially in size, with an average of 10.4 cosponsors per bill and the largest bill involving 52 cosponsors. Each hyperedge thus captures a natural group interaction among the legislators who collaborated on a bill, producing a network with widely varying interaction sizes.

The results, summarized in Fig. 9, confirm that the proposed model and intervention strategies remain effective beyond the synthetic uniform hypergraph setting. As shown in Fig. 9 (a), consistent with the findings on synthetic hypergraphs, higher values of the sensitivity parameter \(\eta\) suppress the endemic state by rapidly reducing hyperedge activity in response to elevated infection pressure, thereby limiting group interactions and expanding the disease‑free region. The comparison of immunization strategies in Fig. 9 (b) reveals that spontaneous isolation (SI\(_1\)) remains highly effective on the general hypergraph with heterogeneous edge sizes, outperforming targeted immunization (TI) and showing a substantial margin over random edge immunization. This indicates that even in realistic higher‑order contact patterns, locally triggered deactivation of highly infected groups can serve as a robust containment mechanism. The contrast between the two post‑intervention rewiring strategies, represented in Fig. 9 (c), shows that degree‑preferential rewiring achieves a modest but consistent reduction in stationary prevalence relative to random rewiring. Under the current setup, new hyperedges are generated with a size equal to the mean edge size of the original network, which is approximately \(10\). The presence of isolated nodes and small hyperedges of size \(2\) or \(3\) in the original network constrains the structural advantage of degree‑preferential rewiring, as these low‑degree nodes provide limited opportunities for the emergence of the heavy‑tailed degree distributions that drive adaptive suppression. This observation suggests that further refinements of the rewiring rules, possibly incorporating infection‑aware or size-dependent mechanisms, may be needed for general hypergraphs.

4 Conclusion↩︎

In this work, we establish an adaptive framework for epidemic dynamics on adaptive hypergraphs by integrating simplicial SIS processes with adaptive hyperedge activity, supported by MMCA analysis and quasistationary Monte Carlo simulations. Our results demonstrate that higher-order interactions fundamentally reshape epidemic spreading through adaptive feedback, giving rise to discontinuous phase transitions, bistability, and hysteresis. The theoretical analysis further reveals nonlinear epidemic thresholds in genuinely higher-order settings, where the disease-free state remains stable against infinitesimal perturbations and endemic states emerge only via saddle-node bifurcations. Moreover, intervention strategies that exploit the adaptive structure exhibit marked improvements in containment performance. Targeted immunization based on infection-pressure ranking and low-threshold spontaneous isolation consistently outperforms random baselines, while degree-preferential rewiring achieves long-term suppression by inducing heterogeneous topological patterns that reinforce adaptive feedback. In contrast, random rewiring provides only short-lived benefits due to the rapid restoration of transmission pathways.

Despite these advances, the present framework assumes uniform adaptive parameters and a single-layer representation, overlooking real-world heterogeneity. In particular, all hyperedges share the same sensitivity and recovery rates, whereas actual groups likely differ in their responsiveness, potentially yielding richer dynamics such as coexisting endemic and disease-free subgroups. However, since the IBMF framework tracks node probabilities independently, our model can naturally be extended to more general structures—for example, by introducing size-dependent transmission rates \(\beta_d\) and a hyperdegree distribution \(P(k_d)\) for hypergraphs with heterogeneous hyperedge sizes, or by accounting for all infection-triggering subsets within simplices in simplicial complexes. Such extensions introduce additional nonlinear terms while preserving the overall MMCA structure.

In summary, our findings highlight the central role of higher-order interactions, adaptive connectivity, and targeted interventions in shaping epidemic dynamics and enhancing controllability on adaptive hypergraphs. Beyond the macroscopic behaviors uncovered in our analysis, the simulations show that highly connected nodes accumulate disproportionately high infection pressure, which in turn triggers stronger adaptive suppression of their incident hyperedges. This localized suppression cascades through the system and significantly amplifies the effectiveness of interventions that specifically target structurally vulnerable or behaviorally responsive groups. Building on these insights, future work could incorporate information-driven decision rules in which awareness, perceived risk, or digital signals dynamically modulate hyperedge activity or integrate optimal intervention scheduling under limited resources. Such developments would deepen the theoretical understanding of adaptive contagion and support the design of more effective public health policies in complex social systems.

References↩︎

---

1. Corresponding author. Email: qinli1022@swu.edu.cn↩︎

2. Corresponding author. Email: myfeng@swu.edu.cn↩︎

3. This work is supported in part by the Natural Science Foundation of Chongqing under Grant CSTB2025YITP-QCRCX0007, in part by the Fundamental Research Funds for the Central Universities under Grant SWU-KT26010, and in part by Slovenian Research and Innovation Agency under Grant P1-0403.↩︎