Multi-disease transmission with coupling and threshold perturbations on hypergraph structure

1 Introduction

Multidisease transmission is a fundamental issue in the dynamics of infectious diseases. Conventional compartmental models, such as SIR and SIS, are predicated on the homogeneous mixing assumption, whereby contacts between individuals are regarded as random and independent [1]. While this assumption simplifies the analysis, it is not able to capture the heterogeneity of contact that is prevalent in real-world populations. Pioneering research proposed by Newman, Pastor-Satorras, and Vespignani demonstrated that incorporating network structure fundamentally alters transmission dynamics. This phenomenon is particularly pronounced in scale-free networks, where transmission thresholds can become undetectable, thereby enabling diseases to persist even at very low transmission rates [2].

Interpersonal contact patterns in real-world settings are characterised by considerable heterogeneity, in addition to the manifestation of pronounced organisational structures. Contact frequencies exhibit significant variation among individuals, with highly connected individuals, such as healthcare workers or those with a high social activity level, often assuming pivotal roles in the transmission process. Concurrently, a considerable volume of social interactions occurs within group settings such as households, classrooms, and workplaces, involving simultaneous participation by multiple individuals. The characterisation of such collective interactions using pair-wise connections is a challenging task. The notion of hypergraph was initially introduced by Berge in 1973, providing a natural tool for the description of higher-order interaction structures [3]. Recent studies have demonstrated that propagation processes based on hypergraph can exhibit dynamical behaviours absent in ordinary networks, such as discontinuous phase transitions and outbreak propagation [4], further revealing the substantial impact of higher-order interactions on propagation dynamics [5, 6].

Infectious diseases frequently coexist and spread concurrently in the real world. In the contemporary context of increasing population mobility and environmental change, the occurrence of simultaneous outbreaks of diverse diseases within the same population has become increasingly prevalent. For instance, during the course of the pandemic caused by the severe acute respiratory syndrome (SARS) virus, influenza viruses and respiratory syncytial viruses continued to circulate, with cases of co-infection being reported [7]. The research indicate that diseases may interact through mechanisms such as immune suppression, symptom overlap, or resource competition [8, 9], whereby the transmission status of one disease alters an individual’s susceptibility or transmissibility to another disease.

The research on multi-disease transmission can be traced back to the early work of Castillo-Chavez et al. [10], and has evolved to include competitive models [11], cooperative models [12], and co-transmission studies [13]. The existing results still predominantly rely on averaged or pairwise contact network frameworks [14–16]. In these models, higher-order group contacts are typically neglected, and interactions between diseases are often introduced only through parameter adjustments. There are few researches incorporating higher-order contact structures and multi-disease coupling mechanisms into transmission dynamics analysis. In their study on the dynamics of epidemic spread in multi-relationship networks, Li Ruiqi et al. extended the study of epidemic spread from single-relationship networks to multi-relationship networks [17], but considered only two types of relationships: work and friends. Real-world social relationships are more complex, and the model did not fully capture them. The construction of models capable of characterising both types of transmission characteristics is crucial for understanding real-world epidemiological [18–20]. The remainder of this paper is organized as follows.

2 Models and methods

2.1 Hypergraph construction

Traditional graph models cannot capture concurrent interactions involving multiple nodes, which can result in the initial outbreak velocity and infection peak being underestimated. Although scale-free networks can reflect hub nodes, they struggle to describe group contacts. Similarly, random hypergraphs possess concurrency features, but lack structural heterogeneity. The ER random graph, proposed by Erdős and Rényi in 1959 [21], is the most fundamental classical model in network science. Given N nodes, this model connects each pair of nodes with a probability p, resulting in a network degree distribution that follows a Poisson distribution. The ER model laid the theoretical foundation for the study of complex networks. The BA scale-free network was proposed by Barabási and Albert in 1999 [22]. This model introduces two mechanisms: growth and preferential attachment. New nodes are continuously added to the network, and these new nodes tend to connect to nodes with higher degrees. The degree distribution of the BA model follows a power-law distribution and can explain the scale-free properties of many real-world networks. In order to better approximate real-world human contact patterns, this model must simultaneously retain group concurrency and individual heterogeneity. Hypergraphs based on ER and BA mechanisms are constructed to characterise random and heterogeneous contact structures, respectively. Furthermore, to reflect variations in population size distribution, ER and BA hypergraphs are categorised as either uniform or non-uniform. This yields four hypergraph structures—BA uniform, BA non-uniform, ER uniform and ER non-uniform—for comparative analysis.

2.1.1 BA hypergraph

The construction of BA supergraphs aims to reflect the scale-free properties and high clustering characteristics of social networks. First, an underlying structure of a BA scale-free network with an average degree of ⟨ K ⟩ = 6 is generated, ensuring a power-law distribution of node degrees. Hyperedges are then generated by identifying triangular structures within the network. For example, if a node v and its neighbours u and w are all connected, a four-node hyperedge { v , u , w , k } is constructed to simulate group interactions, such as those within families, social circles or work teams. Uniform BA hypergraphs are constructed using hyperedges of fixed sizes, ranging from two to 5. Non-uniform BA hypergraphs incorporate hyperedges of sizes ranging from two to five simultaneously. This approach preserves the heterogeneity of BA networks while introducing structures that facilitate higher-order group interactions. As shown in Figure 1.

2.1.2 ER hypergraph

ER hypergraphs are used to model random encounter scenarios. In uniform ER hypergraphs, all superedge sizes are fixed (e.g., k = 4), representing group activities with a consistent participant count. For comparison, this paper also constructs uniform ER hypergraphs ranging from k = 2 to k = 5. Non-uniform ER hypergraphs allow superedges of different sizes to coexist. The construction process ensures that no superedges or inclusion relationships are duplicated, thereby avoiding structural redundancy.

In Figure 2, there are three hyperedges, H 1 { A , B , C }, H 2 { A , C , D }, and H 3 { B , D , E }. Nodes B and C are infected (red), nodes D and E are susceptible (green), and the observed node A is also susceptible (blue). First, we determine the neighbours of node A: in H 1, node A is adjacent to B and C; in H 2, it is adjacent to C and D. Thus, the neighbour multiset of A is { B , C , C , D }. Since B and both occurrences of C are infected, A receives infection attempts from B + C + C, i.e., three times in total, while D does not contribute to infection. Therefore, the infection probability of node A is 0.75.

2.2 Coupling threshold propagation model

2.2.1 Basic framework

In this paper, based on the SIS model, we propose the coupling threshold propagation model (CTPM) to investigate the transmission mechanisms of multiple diseases. In reality, epidemics such as dengue fever, measles and influenza can result in repeated infections due to transient or absent immunity following recovery. Current models struggle to accurately describe these processes. For example, the SIR model assumes permanent immunity, whereas the SIS model effectively captures the characteristics of repeated infections. The CTPM comprehensively accounts for differences in initial transmission rates among diseases, interactions between diseases, and heterogeneity in individual infection thresholds. Specifically, the infection status of any node v and disease i ( i ∈ { 1,2 } ) at time t is denoted as follows:

I i , v t ∈ 0,1 (1)

In this model, I i , v ( t ) = 0 indicates that node v is in the susceptible state for disease i, and I i , v ( t ) = 1 indicates that node v is in the infected state for disease i. The population is divided into the four groups, where S is the susceptible population, I 1 is the population infected with disease D1, I 2 is the population infected with disease D2, and I 3 is the population co-infected with both diseases D1 and D2. Then we would express this as:

S + I 1 + I 2 − I 3 = 1 (2)

The kinetic equation is: M 1 represents nodes infected only with virus D1, while I 1 represents nodes infected only with D1 and those infected with both D1 and D2; similarly, M 2 represents nodes infected only with virus D2, while I 1 represents nodes infected only with D2 and those fected with both D1 and D2. β i is the transmission rate, and γ is the recovery rate.

d M 1 d t = β 1 S I 1 + γ I 3 − β 2 M 1 I 2 − γ M 1 d M 2 d t = β 2 S I 2 + γ I 3 − β 1 M 2 I 1 − γ M 2 d I 3 d t = β 1 I 1 M 2 + β 2 I 2 M 1 − 2 γ − γ 2 I 3 (3)

2.2.2 Infection process

The epidemic propagation process which is simulated in this paper is illustrated by Figure 3. For a susceptible node v and disease i, node v becomes infected with probability β when the proportion of infected neighbours ϕ exceeds its infection threshold θ. This transition occurs from the susceptible state to the infected state. Infected nodes then recover to the susceptible state with a certain probability, after which they can be reinfected. And an initial propagation rate combination was designed, as shown in Table 1.

TABLE 1. Three initial propagation models of β base.

2.2.3 Actual transmission rate

Based on practical considerations, the CTPM model incorporates the time-dependent effects of enhanced human intervention. In the initial phase, the transmission process is often uncontrolled due to constraints in medical resources and awareness levels. However, as monitoring, drug development and non-pharmaceutical interventions improve, transmission risks are proactively reduced. To capture this, the paper incorporates a time-decay term into the transmission rate parameter to reflect the effect of sustained control measures such as isolation, mask usage and vaccination in suppressing transmission. This design reflects the gradual weakening of transmission intensity under intervention, thereby providing a more realistic simulation of transmission dynamics and control effectiveness across different phases.

The time-decaying function in the context described above can be expressed as following: λ is the time decay rate parameter, which is set to 0.05. t is the time step, e is the natural constant (approximately 2.71828).

y t = e − λ t (4)

The actual transmission rate of the corresponding disease β i ( t ) varies with y ( t ) and is expressed as:

β i t = β base ⋅ e − λ t (5)

2.2.4 Probability of contact

To characterize the local infection pressure in the surrounding environment, let φ i , v ( t ) denote the proportion of neighbours of node v that are infected with disease i at time t. This is referred to as the “local infection ratio”, which reflects the infection density around node v. The higher the proportion, the greater the infection risk in the contact environment of node v.

ϕ i , v t = | ∑ v ∈ N v I i , v t = 1 | | N v | (6)

Here, N v represents the set of unique neighbours of node v after removing duplicates. In this paper, a hypergraph is used to represent higher-order contact relationships. Therefore, the neighbours of node v are derived from all hyperedges that contain v. For each hyperedge containing v, we take the set of all nodes within that hyperedge, merge these sets, and then remove node v itself. After eliminating duplicates, the resulting set is N v, which represents the unique set of neighbours of node v.Let e ( v ) denote a hyperedge containing node v, and ε ( v ) represent the set of hyperedges connected to e ( v ), and V ( e ( v ) ) denotes the set of nodes contained in hyperedge e. Based on this, ϕ i , v ( t ) is calculated as the ratio of the number of infected neighbours to the total number of neighbours. Within the propagation mechanism, a node transitions from the susceptible state to the infected state only if both the threshold criterion and the stochastic propagation criterion are satisfied. First, the threshold condition is defined as follows: at time t, if the local infection proportion ϕ i , v ( t ) of node v exceeds its corresponding infection threshold θ i , v ( t ), it indicates that the infection pressure in its contact environment has reached a level capable of triggering infection risk. That is, the condition is expressed as:

N v = ⋃ e v ∈ ε v V e v − v (7)

ϕ i , v t > θ i , v t (8)

When the threshold criteria are met, in order to reflect the randomness and uncertainty of transmission, it is also necessary to introduce random transmission conditions, generate a random number for the interval [0, 1], and compare it with the actual transmission rate β i ( t ) of disease i at time t.and then spread according to the actual transmission rate of the disease β i ( t ), there is a probability of β i ( t ) transmission successfully.

random 0,1 ≤ β i t (9)

Therefore, the update rule of infection status of node v for disease i at time t can be summarized as follows: when and only when the proportion of local infection exceeds the threshold, and the random transmission experiment is successful, the node changes from susceptibility to infection, that is, the value of I v , i ( t ) is 1; Otherwise, the original state remains unchanged, which not only ensures that the transmission is driven by local infection pressure, but also reflects the random transmission controlled by the transmission rate parameter.

I v , i t = 1 , when random 0,1 ≤ β i t and ϕ i , v t > θ i , v t 0 , otherwise (10)

2.2.5 Dynamic threshold normalization

In multi-level hyperedge scenarios, scale correction of thresholds is required to ensure comparability of threshold conditions across different hyperedge levels. The same numerical neighbour infection rate (e.g., 50%) represents different actual infection pressures across contact orders: a binary contact requires only one infected neighbour to reach this rate, whereas a quaternary contact requires more infected individuals, making the triggering condition more stringent. Without adjustment, threshold determination will systematically favour lower-order contacts, leading to an imbalance in transmission mechanisms.

This paper addresses this issue by employing a normalization method based on hyperedge order, which standardises thresholds by adjusting them according to ‘order minus one’. This eliminates bias arising from differences in hyperedge scale. This enables thresholds across different hyperedge orders to be compared on the same scale, providing a more accurate reflection of the actual impact of multi-order contacts on transmission risk.

The threshold normalization can be expressed as:where θ 0 is the underlying infection threshold, α ∈ [ 0,1 ] is the normalized index (uniformly set to one in this experiment), c is the coupling strength, I j , v ( t ) indicates whether node v is infected with another disease j (j ≠ i; one if infected, 0 otherwise). θ base = 0.1 is the baseline infection threshold. c ∈ { − 7,0,7 }: c = 7 (positive coupling) reduces the threshold by 40%; c = − 7 (negative coupling) increases the threshold by 40%; c = 0 (uncoupled) no interaction.

θ i , v t = 0.01 , θ 0 N v α ⋅ ∑ v ∈ N v 1 − c ⋅ I j , v t < 0.01 θ 0 N v α ⋅ ∑ v ∈ N v 1 − c ⋅ I j , v t , 0.01 ≤ ⋯ ≤ 0.5 0.5 , θ 0 N v α ⋅ ∑ v ∈ N v 1 − c ⋅ I j , v t > 0.5 (11)

2.2.6 Network dynamic over-edge generation

In social networks, when a topic or behaviour becomes popular, individuals who adopt it tend to interact more frequently, forming new, temporary, tightly knit groups. Traditional propagation models typically operate on static network structures. By dynamically creating new superedges based on infection status, CTPM can simulate such real-world phenomena more realistically. To more clearly represent these phenomena, this paper introduces discrete probability addition: G H ( t ) and G H ( t + 1 ) represent the hypergraph structures at time steps t and t + 1, respectively. e new = { v , u 1 , u 2 }, where v ∈ infected nodes, u 1 , u 2 randomly selected from its neighbours. p create = 0.01.

G H t + 1 = G H t ∪ e new , with probability p create G H t , with probability 1 − p create (12)

2.3 Mathematical derivation of formulas

From:we can derive:Hence:

I i , v t ∈ 0,1 (13)

S + M 1 + M 2 + I 3 = 1 , I 1 = M 1 + I 3 , I 2 = M 2 + I 3 (14)

S + I 1 − I 3 + I 2 − I 3 + I 3 = 1 ⇒ S + I 1 + I 2 − I 3 = 1 (15)

The new infections of M 1 are given by β 1 S I 1 + γ I 3, and the new recoveries from M 1 are β 2 M 1 I 2 + γ M 1. Thus the dynamic equation for M 1 is:

d M 1 d t = β 1 S I 1 + γ I 3 − β 2 M 1 I 2 − γ M 1 (16)

The dynamic equations for M 2 and I 3 can be derived analogously, yielding the system:

d M 1 d t = β 1 S I 1 + γ I 3 − β 2 M 1 I 2 − γ M 1 d M 2 d t = β 2 S I 2 + γ I 3 − β 1 M 2 I 1 − γ M 2 d I 3 d t = β 1 I 1 M 2 + β 2 I 2 M 1 − 2 γ − γ 2 I 3 (17)

To better reflect reality, we introduce a time-decay function for the baseline transmission rate β base:The actual transmission rate then becomes:

y t = e − λ t (18)

β i t = β base ⋅ y t = β base ⋅ e − λ t (19)

Define ϕ i , v ( t ) as the proportion of neighbours of node v that are infected with disease i at time t, i.e., the “local infection ratio”. It is obtained by scanning the infection states of all neighbours of v:In Equation 6, N v denotes the set of unique neighbours of node v after removing duplicates. It is formally defined as:where e ( v ) represents a hyperedge that contains v, ε ( v ) is the set of all hyperedges incident to e ( v ), and V ( e ( v ) ) is the set of nodes contained in hyperedge e. Note that the node v itself is excluded from the neighbour set.

ϕ i , v t = v ∈ N v : I i , v t = 1 | N v | (20)

N v = ⋃ e v ∈ ε v V e v − v (21)

Let ϕ i , v ( t ) be the local infection ratio defined in (Equation 6), and let θ i , v ( t ) be the dynamic infection threshold for node v and disease i. The threshold condition for infection is:

ϕ i , v t > θ i , v t (22)

If the threshold condition is satisfied, a stochastic transmission condition must also hold. The transmission is successful when the following random draw meets the condition:

random 0,1 ≥ β i t (23)

Combining (Equations 8, 9), the update rule for the infection state of node v with respect to disease i at time t is:For the dynamic threshold θ i , v ( t ), let θ 0 be the baseline threshold, c the coupling strength between diseases, and α the normalisation exponent. For each neighbour u of v, the coupling correction term is given by 1 − c ⋅ I j , u ( t ), where I j , u ( t ) indicates whether neighbour u is infected with the other disease j ( j ≠ i ). Summing this correction over all neighbours of v and normalising by ( N v ) α, we obtain the raw threshold value:To reflect realistic constraints, this value is clipped to the interval [ 0.01 , 0.5 ]. The final normalised threshold is thus defined as:

I i , v t = 1 , if random 0,1 ≥ β i t and ϕ i , v t > θ i , v t , 0 , otherwise . (24)

θ 0 N v α ⋅ ∑ u ∈ N v 1 − c ⋅ I j , u t .

θ i , v t = 0.01 , θ 0 N v α ⋅ ∑ v ∈ N v 1 − c ⋅ I j , v t < 0.01 θ 0 N v α ⋅ ∑ v ∈ N v 1 − c ⋅ I j , v t , 0.01 ≤ ⋯ ≤ 0.5 0.5 , θ 0 N v α ⋅ ∑ v ∈ N v 1 − c ⋅ I j , v t > 0.5 (25)

3 Experimental results and analysis

3.1 Simulation experiment

This section presents simulation experiments for the CTPM proposed in Section 2, primarily considering the following factors: network type, hypergraph structure, initial transmission rate, and disease coupling relationships. The pseudocode for experimental construction is presented in Table 2. In the transmission simulation, node states are updated according to time steps. Uninfected nodes determine whether to become infected by comparing the proportion of infected neighbours to a dynamic threshold. Infected nodes recover to the susceptible state with a fixed probability. Network structures evolve dynamically in response to infection status, reflecting adaptive changes in contact patterns. The CTPM system simulates how the dynamics of multi-disease transmission are affected by topology, structure, initial conditions and coupling relationships.

TABLE 2. Experimental pseudocode algorithm.

The experiment was conducted on a hypergraph consisting of 3,000 nodes with an average degree of six to ensure the presence of significant higher-order interactions. The network initialisation employs a dual-disease coexistence mechanism: each disease has an initial infection rate of 5%, and nodes that are infected with more than one disease are permitted. A single-disease transmission group was also set up as a control. The network architecture uses ER and BA networks to construct uniform and non-uniform hypergraphs with degrees ranging from two to 5. This approach evaluates the impact of network heterogeneity, higher-order interactions and hypergraph structure on propagation.

The parameters of the diseases D1 and D2 in the experiment are set in Table 3. The maximum transmission probability for both diseases is set to 0.9, the baseline infection threshold is 0.1 (indicating that an infection risk exists when the neighbour infection proportion exceeds 10%), and the recovery probability is set to 0.15. In the heterogeneous hypergraph, the proportion of hyperedges of orders two to five is as follows: 40% for two-order hyperedges, 30% for three-order hyperedges, 20% for four-order hyperedges, and 10% for five-order hyperedges.

TABLE 3. Experimental parameter settings.

Coupling relationships can be classified into three types: positive coupling (mutual promotion), negative coupling (mutual inhibition) and no coupling (independent transmission). To analyze the impact of coupling on the transmission dynamics of two diseases, the effects of coupling are incorporated into the model through adjustments to the infection thresholds. In reality, the effects of co-infection are multidimensional—for example, 2co-infection with rhinovirus and nffuenza may exacerbate symptoms and prolong disease uration [23] To simplify the analysis, these complex interactions are collectively abstracted in this study as modifications to the threshold parameters.

For different proportion settings considered in the experiments, several evaluation metrics in heterogeneous networks are calculated, including the average excess degree ⟨ K ⟩, the variance var, and the basic reproduction number R 0. The calculation methods are as follows:

Average Superiority: ⟨ K ⟩ = ∑ i q i ⋅ K i (26)

Variance: var = ∑ i q i ⋅ K i − ⟨ K ⟩ 2 (27)

Basic reproduction number: R 0 = β γ ⋅ 〈 K 2 〉 〈 K 〉 = β γ ⋅ 〈 K 〉 + var 〈 K 〉 (28)

Ki represents the degree of node i, and qi represents the proportion of nodes with degree Ki out of the total number of nodes in the network. A higher value of ⟨ K ⟩ implies a larger number of transmission pathways and increased network density, thereby accelerating disease spread. A higher variance ( var ) indicates the presence of more hub nodes, which can enhance transmission efficiency even when ⟨ K ⟩ is held constant. Both ⟨ K ⟩ and var are positively correlated with the basic reproduction number R 0, and a synergistic effect exists whereby the amplifying role of var depends on ⟨ K ⟩. When ⟨ K ⟩ is large, the relative influence of var diminishes, although the combined effect of ⟨ K ⟩ and var leads to a stronger absolute increase in R 0. Conversely, when ⟨ K ⟩ is small, the impact of var becomes more pronounced and may lead to epidemic outbreaks even in low-density networks.

BA and ER networks represent heterogeneous and homogeneous structures, respectively. Using the no-coupling case as a reference for independent transmission, it can be observed that the infection rate in the heterogeneous ER network is consistently higher than that of its k-3 order homogeneous network. In contrast, for BA networks, the infection rate of the k-3 order homogeneous hypernetwork is generally higher than or comparable to that of the heterogeneous case, indicating that under equivalent conditions, homogeneous ER networks exhibit stronger transmission intensity than BA networks. Comparative analysis shows that as the hyperedge order k increases, the peak transmission rate rises markedly, with the most pronounced increase occurring from k = 2 to k = 3, after which the growth rate gradually diminishes as k continues to increase. Because real-world networks more closely resemble BA-type structures, the infection rate of their heterogeneous networks typically lies between those of the k = 2 and k = 3 order homogeneous networks.

Considering the heterogeneity of real-world interactions, this study defines various non-uniform networks with different hyperedge proportions to investigate their impact. The control group includes proportions to as follows: (3, 4, 2, 1), (3, 2, 4, 1), (3, 2, 1, 4), (2, 5, 2, 1), (2, 2, 5, 1), (2, 2, 1, 5), (2, 6, 1, 1), (2, 1, 6, 1), (2, 1, 1, 6), (1, 7, 1, 1), (1, 1, 7, 1), and (1, 1, 1, 7).

In the first four groups, the valence in the ER non-uniform network is consistently higher than that in the k-3 uniform network, as shown in Figures 4, 5. In group, the BA non-uniform network initially rises more slowly than the k-3 uniform network, but reaches a higher peak, eventually approaching the ER non-uniform curve and continuously suppressing the k-3 uniform network in the later stages. In group, the BA non-uniform network initially matches the k-3 uniform network, surpasses it near the peak, and ultimately aligns with the ER non-uniform network. In group, the three curves distinctly separate: ER non-uniform > BA non-uniform > k-3 uniform, with the difference in prevalence ordered as > >. In group, the ER and BA non-uniform curves become closer to each other and more similar to the k-4 uniform network, with the difference showing as >.

It is noteworthy that even when the highest-order hyperedge proportion reaches 40%, the peak prevalence of the non-uniform network remains lower than that of the k max − 1 uniform network. To further explore the role of higher-order structures, the experiment continues to increase the highest-order proportion and analyzes its impact on transmission.

In comparison with the last nine experimental groups, the highest-order proportion is first increased by 10% (groups –). Compared to groups –, a 10% increase in the highest-order proportion leads to a significant rise in prevalence in the non-uniform network at the same highest-order level. However, the peak prevalence shows the result of >. This is because the proportion of k-4 in group is higher than in group, and higher-order structures have a stronger amplifying effect on transmission.

In groups –, the proportion of the highest-order structure is increased to 7, while the proportions of all other orders are fixed at 1, in order to eliminate interference from other structural orders. In group, the prevalence in heterogeneous networks is substantially higher than that of the k-3 homogeneous network, and the peak prevalence follows the ordering: heterogeneous > k max homogeneous. The non-uniform prevalence is significantly higher than that of the k-3 uniform network, with the peak showing that non-uniform > k max uniform. In group, only the ER heterogeneous curve under D2 barely reaches the k max homogeneous case, whereas all other heterogeneous curves remain below the k max homogeneous curve and are overall lower than those in group, indicating that a higher proportion of high-order structures exerts a more pronounced promoting effect on transmission. All heterogeneous curves are higher than the k max − 1 order homogeneous case, but still remain below the k max homogeneous network. These results indicate that as the structural order increases, its marginal contribution to the enhancement of the global infection rate gradually diminishes, which is consistent with the expectation that heterogeneous networks remain below the k max homogeneous curve.

The experiments show that increasing the proportion of higher-order structures leads to a higher transmission peak. When the proportions of k i and k j are equal, the transmission peak associated with max ( k i , k j ) is higher. This suggests that epidemic control strategies should prioritize interventions targeting high-order homogeneous structures ( k ≥ 3 ) by introducing heterogeneity, rather than excessively focusing on whether other high-order structures remain present within heterogeneous networks. If conditions permit, these structures can be further reduced to k = 2 order ordinary networks, and network connectivity can be lowered through measures such as isolation, thereby suppressing disease transmission.

From the perspective of the initial transmission rates, the largest disparity arises when one disease has a high initial transmission rate while the other has a low one, under which the mutual inhibitory effect is most pronounced. This further confirms that the greater the difference in transmission rates, the more evident the impact of coupling effects.

Based on the calculated results in Table 4 and the experimental observations shown in Figures 6, 7, it can be seen that when var is relatively low, as in cases and, the spreading process is slow. As ⟨ K ⟩ increases, R 0 exhibits a slight increase, accompanied by a corresponding but limited enhancement in the transmission rate. By contrast, when both var and ⟨ K ⟩ are at relatively high levels, as in group, R 0 increases significantly, the overall prevalence rises markedly, and the growth process becomes much more rapid. A further comparison between cases and reveals that if var and ⟨ K ⟩ are simply added, the resulting value for case (5.59) exceeds that for case (5.44). However, the resulting R 0 shows the opposite ordering. Experimental results shown in Figures 6, 7 further indicate that the overall prevalence in case is higher than that in case. It demonstrates that var and ⟨ K ⟩ do not contribute to R 0 through a simple linear superposition, but instead act jointly through a synergistic effect that mutually reinforces their influence.

TABLE 4. ⟨ K ⟩ and variance effects on R 0.

Another noteworthy phenomenon is reflected in the competitive dynamics between diseases: under no coupling conditions, in the ER network, disease D 2 suppresses D 1 in k-5 to k-3 uniform networks, while in k-2 and non-uniform networks, D 1 suppresses D 2. After introducing positive and negative coupling, all competitive patterns in ER networks are reversed, while those in BA networks remain unchanged. This can be attributed to the heterogeneity of BA networks, which disperses the uniformly propagating coupling effects and thus preserves the initial competitive hierarchy. In contrast, the homogeneity of ER networks allows coupling effects to act globally, thereby altering the overall system balance. When the initial transmission rates differ markedly, the disease with the higher transmission rate exhibits strong suppressive dominance, and the competitive pattern remains unchanged. The transition in competitive patterns originates from the system-wide shift induced by coupling effects in homogeneous networks. This raises a key question: how does the transition in competitive patterns occur?

In BA networks, the initial structure follows a typical power-law distribution, characterized by a small number of hub nodes and a large number of low-degree nodes. Once D 1 becomes dominant, positive coupling rapidly amplifies the super-spreading capacity of hub nodes through positive feedback, driving the transmission process toward monopolization and granting D 1 a sustained advantage. Negative coupling, by contrast, compresses transmission pathways through hierarchical suppression, thereby rigidifying the transmission structure. Despite their opposing mechanisms, both forms of coupling ultimately stabilize the competitive pattern.

In ER networks, the node degree distribution is homogeneous, and no prominent hub nodes are present. Although positive coupling temporarily accelerates transmission, the absence of super-spreaders causes the advantage of D 1 to be rapidly weakened by resource allocation mechanisms. Under negative coupling, mutual inhibition between the diseases is intensified, significantly suppressing the spread of D 1 and ultimately leading to a reversal of the competitive pattern.

In the experiments as shown in Figure 8, all black curves represent the proportion of co-infection between D 1 and D 2 in heterogeneous networks, and they exhibit clear differentiation. The peak prevalence in ER networks is consistently higher than that in BA networks. Negative coupling leads to a reduction in the peak prevalence for both types of networks. Under low-coupling conditions, the BA network surpasses the ER network in the later stage after the peak, because a small number of infected hub nodes are retained during the decay phase, thereby sustaining transmission, whereas ER networks lack such nodes. In addition, under negative coupling, the rising and declining rates of all curves are markedly slowed, which further enlarges the differences between different k orders. This indicates that the sensitivity of higher-order structures to coupling effects continuously increases with the structural order.

Phase diagrams in Figures 9–12 show that, in both BA and ER networks, the proportion of I 3 increases with the hyperedge order k, indicating that higher-order structures facilitate disease co-infection.

This increasing trend with respect to k remains consistent across different scenarios. One possible explanation is that the increase in co-infection is driven by rising single-disease infection rates. Under asymmetric negative coupling conditions, when the infection proportion of D 1 increases significantly, the proportion of I 3 rises accordingly. Under symmetric high negative coupling, the proportion of I 3 exhibits an even more pronounced increase, supporting this hypothesis. Extending the analysis to symmetric positive coupling scenarios, the I 3 proportion also increases markedly in all networks, indicating that positive coupling likewise promotes co-infection. Therefore, the growth of I 3 results from the combined effects of enhanced single-disease infection rates and coupling interactions.

When disease transmission continuously declines and remains at an extremely low level for an extended period, it is classified as extinction as shown in Figures 13 and 14. Extinction is observed only in the k = 2 order homogeneous networks. Analysis of the extinction time steps indicates that extinction is closely associated with the coupling relationships. In ER networks, the extinction time follows the ordering: negative coupling < positive coupling < no coupling. The I 3 curves in BA networks exhibit the same ordering, a result that appears counterintuitive at first glance.

This phenomenon can be attributed to differences in network topological properties: ER networks exhibit a homogeneous degree distribution with random connections and lack prominent transmission hotspots. Negative coupling directly weakens or even blocks transmission pathways, thereby accelerating extinction. Moreover, due to the low local clustering coefficient, transmission in ER networks relies strongly on global connectivity. Positive coupling accelerates the depletion of susceptible nodes, which likewise leads to earlier extinction.

Intuitively, in BA networks, super-spreaders and highly clustered communities can respectively control transmission through central nodes and form closed transmission subcommunities. Negative coupling weakens the capacity of central nodes and is therefore expected to accelerate extinction, whereas positive coupling strengthens their influence and prolongs transmission within communities, which would intuitively lead to the slowest extinction. However, the observed overall extinction ordering is: no coupling (fastest), positive coupling (slowest), with negative coupling in between. The delayed extinction under negative coupling arises because, after D 1 infects nodes in BA networks, negative coupling strongly suppresses the reinfection of D 2, locally blocking transmission and generating isolated small clusters, which require longer times to die out. Under asymmetric negative coupling conditions, where transmission abilities differ markedly, the effect of negative coupling dominates, causing D 2 to go extinct much faster than in other scenarios. The extinction time of I 3 in BA networks is consistent with that observed in ER networks, indicating a pattern of cooperative extinction.

Cooperative extinction refers to a process in which the elimination of one disease accelerates the elimination of the other. In theory, positive coupling should make diseases more difficult to eliminate by lowering the infection threshold; however, the experimental results show the opposite trend, with extinction occurring faster than in the no-coupling case. Further analysis indicates that this process can be divided into three stages. First, the reduced threshold enables the remaining infected individuals to deplete susceptible individuals more rapidly. Second, as the transmission rate continues to decline, the number of new infections becomes lower than the number of recoveries, leading to a contraction of the susceptible population. Finally, when the proportion of susceptible individuals falls below the threshold, the disease begins to go extinct. In addition, local communities in BA networks may sustain transmission chains, implying a potential risk of resurgence.

3.2 Empirical experiment on dengue fever

A real underlying network was constructed based on a dengue fever transmission dataset from Ecuador. Following the experimental procedure described in Section 3.1, diseases D 1 and D 2 were allowed to undergo coupled spreading on the corresponding homogeneous and heterogeneous hypergraphs of different orders generated from the dengue-based underlying network, and the experimental results were systematically analyzed. The dengue fever simulation differs somewhat from the previous section; the distinction lies in the construction of the underlying network. In the previous section, a synthetic network was used to generate a hypergraph, whereas in the dengue fever simulation, a network built from a real-world dengue fever dataset serves as the underlying network for generating the hypernetwork. In this section, there is no longer a distinction between network types, as networks constructed from actual datasets possess their own unique properties.

In the empirical dengue fever network, as shown in Figure 15 D 1 is only weakly affected by the coupling type in k = 3 to k = 5 order homogeneous networks. It was initially hypothesized that this phenomenon results from the constraints imposed on the infection rate by setting the upper bound of the transmission rate to 0.9 and the threshold to 0.5. However, in cases d, e, and f, the peak infection rate is only 0.5, and no noticeable differences are observed among higher-order k cases ( k ≥ 3 ), thereby refuting this hypothesis. A comparison across all nine scenarios further shows that the peak prevalence of all high-order homogeneous networks does not vary significantly with the coupling type.

Further control experiments that removed the upper bound on the threshold show that the upper limit of the infection rate does not increase accordingly. This is because the time-dependent function used to simulate real-world interventions constrains the transmission rate, making large-scale spreading difficult even when the threshold condition is satisfied. Under these conditions, the peak prevalence of k ≥ 3 homogeneous networks remains indistinguishable across the three coupling types. This indicates that high-order homogeneous networks exhibit strong robustness, and that coupling relationships have a limited influence on their transmission dynamics.

By contrast, coupling exerts a pronounced influence on k = 2 order homogeneous and heterogeneous networks. Under negative coupling, the peak prevalence is markedly reduced. Under positive coupling, the initially dominant disease can alter the competitive pattern, allowing the other disease to overtake it. In contrast, negative coupling locks in the dominance of the initially advantaged disease and suppresses the spread of the other disease. As illustrated by the comparison between cases g and h, in heterogeneous networks D 2 rapidly surpasses D 1 in case h, whereas in case g it remains dominant throughout, severely constraining the increase in D 1 infection prevalence and highlighting the suppressive effect of negative coupling.

Consistent with the experimental settings in Section 3.1, the proportions of 2-, 3-, 4-, and 5-hyperedges in heterogeneous networks range from case (4, 3, 2, 1) to case (1, 1, 1, 7). Figure 16 compares cases–and shows that, as the proportion of k max increases, the peak prevalence of the corresponding heterogeneous networks rises significantly. Under weak negative coupling, increasing the proportion of k max similarly induces a reversal of the competitive pattern, because an excessive number of hub nodes causes rapid network saturation, thereby amplifying hub effects. This trend is more pronounced in the empirical network than in the synthetic networks, indicating that real-world networks are more sensitive to variations in the proportion of higher-order structures.

From a longitudinal comparison, in Figure 17 it is observed that as the proportion of k max gradually increases, the prevalence in heterogeneous networks rises significantly. A cross-sectional comparison further reveals that as the order k max increases, the rate at which prevalence increases in heterogeneous networks becomes faster. Comparing cases and, where var remains unchanged, only ⟨ K ⟩ is increased, resulting in an increase of R 0 by 0.57 and a substantial rise in prevalence. By contrast, although R 0 increases by 0.69 in case relative to case, the corresponding enhancement in prevalence is weaker than that observed between cases and. This indicates that prevalence is closely related to R 0; however, due to diminishing marginal returns in the increase of R 0, a saturation effect emerges.

The results are consistent with those observed in the synthetic networks, except for one anomalous case: >. The underlying reason is consistent with the previous analysis and can be attributed to the more pronounced influence of higher-order k structures.

In the dengue empirical network, as shown in Figure 18 I 3 exhibits characteristics consistent with those observed in the synthetic networks. Negative coupling not only reduces the peak prevalence across different k orders, but also significantly slows the growth rate of infection, further enlarges the disparities among different k orders, and accelerates the decline of prevalence in the later stage. Although positive coupling slightly increases the peak prevalence, it shortens the duration of the peak and accelerates the extinction process of the disease.

Figure 19 shows that as the order k increases, the co-infected region I 3 (purple area) expands markedly, corroborating the simulation result that higher k orders lead to an increased proportion of co-infection. A comparison with Figure 20 reveals that under positive coupling, the proportion of I 3 is higher, which significantly amplifies the differences among different k orders. This effect is relatively weak for the k = 2 homogeneous network but becomes pronounced for higher-order k networks, further indicating that higher-order structures are more sensitive to coupling effects.

The empirical dengue network exhibits characteristics of cooperative extinction as shown in Figure 21; however, the differences among positive coupling, negative coupling, and no coupling are more pronounced than those observed in the synthetic ER network. Combined with the preceding analysis of competitive pattern transitions, the empirical network can be regarded as lying between an ideal BA network and an ideal ER network. This result indicates that the extinction behavior in the empirical network is not fully equivalent to that of ER networks. Instead, because the network heterogeneity does not reach the level of a fully developed BA network, some effects that would be evident in BA networks are attenuated or masked in the empirical network.

4 Conclusion

The construction of a threshold-triggered hypergraph multi-disease transmission model is presented in this paper, with a comprehensive examination of the effects of network structure, disease coupling, and environmental dynamics on disease spread. The hypergraph model provides a more accurate representation of real-world population interactions, while the threshold mechanism enhances the interpretability of the transmission process.

Experimental findings suggest that the prevalence peak in ER networks is marginally higher than in BA networks. Furthermore, competitive dynamics are more susceptible to reversal due to coupling in ER networks, remaining relatively stable in BA networks. Network heterogeneity has been demonstrated to influence the amplification of coupling effects. It has been shown that the structure of BA networks provides survival space for disadvantaged diseases and enhances the system stability.

The augmentation of higher-order super-edges has been demonstrated to engender a substantial escalation in peak prevalence rates and proportions of co-infection, concomitant with a delay in disease resolution and an enhancement in transmission persistence. Once higher-order structures attain a certain threshold, the networks become more sensitive to coupling but exhibit diminishing marginal returns due to network saturation effects. This phenomenon exemplifies the emergent property, whereby microscopic higher-order interactions initiate qualitative alterations in macroscopic patterns.

The coupling is quantified through the adjustment of dynamic infection thresholds, which reflects diseases interact through immune responses or resource competition in the real-world mechanism. The phenomenon of “synergistic regression” observed in experiments indicates that coupling can drive multiple diseases to stability or extinction, providing a mechanistic explanation for the coexistence of multiple diseases.

Furthermore, the study found that the proportion of higher-order hyperedges in a network significantly impacts the disease prevalence rates. This finding suggests that it is inadequate to characterise the real-world risks only relying on average contact intensity, emphasising the necessity for the integration of higher-order structural features into assessment metrics. By combining network structural metrics and propagation dynamics, we can further establish a basis for estimating critical thresholds. This facilitates the prediction of epidemic risks under known network structural conditions to develop intervention strategies that are better aligned with structural characteristics.

This study provides a unified explanation for disease transmission from three perspectives including structure, mechanism, and outcome. The hypergraph structure determines exposure status, while thresholds and coupling dictate triggering and interaction patterns. Together, they shape the macro-level dynamics of multi-disease co-occurrence. The framework’s parameters are of clear epidemiological significance, readily interface with real-world data, and can be extended to general multi-disease systems, providing a theoretical foundation for controlling multi-disease transmission. This model assumes that the coupling strength and recovery rates are constant. In reality, however, these rates may vary according to the severity of the disease. Although transmission rates are modelled linearly, actual rates may exhibit step changes or oscillations. Hyperedge structures are fixed, which fail to capture the complexity of real-world population scales, spatial and mobility factors.

Within this research framework, the following new approaches to public health prevention and control can be proposed: First, utilize hypergraph structures to identify key hyperedges as intervention targets—rather than focusing on individuals or average exposure—to implement precise structural disruption; second, formulate asynchronous quarantine sequencing strategies based on differences in the strength of multi-disease coupling and recovery rates, thereby disrupting the coupling trigger conditions by delaying the intervention window for a specific disease and reducing the risk of co-outbreaks; Third, by using critical thresholds estimated by the model and combining them with real-time data to calculate how close the current network is to an outbreak, dynamic alert thresholds can be set in advance, rather than relying on fixed thresholds.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Xiaopeng Li, Guizhou University of Finance and Economics, China

Reviewed by: Yijun Ran, Guizhou Normal University, China

Zhishuang Wang, Wuyi University, China