Yuqi Qiao, Haiying Wang
Business School, University of Shanghai for Science and Technology, Shanghai, China.
DOI: 10.4236/ojapps.2025.1512253   PDF    HTML   XML   78 Downloads   351 Views   Citations

Abstract

The study of spreading dynamics in two-layer network systems is of great importance for understanding various real-world propagation phenomena, such as disease transmission and information diffusion. In this paper, we investigate the impact of higher-order structures—specifically, triangular interactions—on contagion processes within a two-layer network, by integrating the SIS epidemic model with a simplicial complex framework. We derive mean-field equations to characterize the system dynamics, analyze the steady-state solutions and their stability, and validate the theoretical findings through numerical simulations. Our results demonstrate that the presence of higher-order structures significantly enhances the spreading capability and induces bistability phenomena. Furthermore, both the inter-layer coupling strength and the asymmetry in node population size are shown to influence the infection density in each layer. This study provides theoretical insights into higher-order contagion mechanisms in multilayer networks and offers novel perspectives for designing effective spread control strategies.

Keywords

Cooperative Games, Multilayer Networks, Higher-Order Interactions, Propagation Dynamics, Bistability, Mean-Field Theory

Share and Cite:

1. Introduction

Complex network theory has become a crucial framework for studying propagation processes, finding extensive applications in fields such as epidemiology, social contagion, and information diffusion [1]-[3]. Traditional contagion models predominantly rely on pairwise interactions between nodes, often overlooking the potential impact of group interactions (e.g., triangular structures) on propagation dynamics. In reality, contact transmission in the real world frequently involves interactions among three or more individuals, forming clustered transmission through multiple concurrent contacts. Examples include small-scale household transmissions, medium-scale transmissions in workplaces or schools, and large-scale transmissions in public venues. From a network perspective, the simultaneous contact of multiple individuals forms higher-order structures within a network, which can be characterized by simplices or hyperedges. Networks incorporating such structures are broadly termed higher-order networks [4], primarily categorized into two types: simplicial complexes [5] and hypergraphs [6]. Evidently, studying epidemic spreading on higher-order networks can transcend the limitations of traditional complex network frameworks. By focusing on higher-order clustering, it more accurately reflects the characteristics of real-world spreading.

In recent years, with the advancement of research on higher-order networks (such as simplicial complexes and hypergraphs), scholars have begun to investigate the role of higher-order structures in contagion processes [7]-[9]. For instance, Iacopini et al. [8] discovered that higher-order interactions can lead to discontinuous phase transitions and bistability in contagion dynamics. Subsequently, Zhao et al. proposed a simplicial contagion model, exploring its discontinuous phase transitions and complex bistable and periodic oscillatory dynamics [10]. Lin et al. established and studied a two-strain SIS competitive model on simplicial complexes [11]. Bukyoung Jhun et al. further considered an SIS model on scale-free simplicial complexes, revealing both continuous and hybrid phase transition behaviors of the disease [12].

On the other hand, real-world spreading processes often occur within multiple interacting networks, such as the coupling between transportation and social networks, or the interaction between online and offline social platforms [13] [14]. These multilayer network structures further increase the complexity of propagation dynamics. Currently, systematic research on higher-order contagion mechanisms in multilayer networks remains insufficient. Propagation behaviors under asymmetric conditions, such as asymmetric node sizes and varying inter-layer coupling strengths, require further in-depth investigation.

This paper aims to study the impact of higher-order structures on propagation dynamics within two-layer networks. We construct a contagion model that incorporates both intra-layer/inter-layer pairwise interactions and triangular higher-order interactions. Using the mean-field approach, we derive the system’s evolution equations, analyze its steady-state solutions and phase diagrams, and combine this with numerical simulations to explore the effects of parameters such as node population size, number of triangles, and inter-layer coupling strength on spreading behavior. Our results reveal the mechanisms through which higher-order structures induce bistability and enhance spreading capability in multilayer networks, providing a theoretical basis for predicting and controlling propagation in real-world multilayer systems.

2. Model

This study deepens the understanding of higher-order propagation dynamics in multilayer complex systems, providing new theoretical perspectives and analytical tools for predicting and controlling real-world spreading processes.

We adopt the mean-field approach to investigate disease spread within the proposed network framework. Based on the SIS contagion model, a node $i$ can reside in one of the following states: $S^A$ , $I^A$ , $S^B$ , or $I^B$ , representing a susceptible or infected node in layer A or layer B, respectively. The schematic diagram of possible infection mechanisms for a node $i$ in a given layer is illustrated in Figure 1, which primarily includes the following pathways:

• Intra-Layer Pairwise Infection: Node $i$ can be infected by an infected neighbor within the same layer through pairwise interactions.

• Inter-Layer Pairwise Infection: Node $i$ can be infected by its counterpart node (or neighboring nodes) in the opposite layer through inter-layer pairwise interactions.

• Intra-Layer Higher-order Infection: Node $i$ can be infected through a group interaction, specifically within a triangular simplex (a 2-simplex) in its own layer. In this scenario, node $i$ experiences an amplified infection pressure when it is connected to two infected neighbors within the same triangle.

These mechanisms collectively define the complex contagion dynamics within our two-layer network model with higher-order interactions.

Figure 1. Simplicial contagion model (SCM). (a) A randomly generated two-layer interconnected network. (b) In the simplicial contagion model (SCM) of order $D=1$ , susceptible and infected nodes are represented by blue and orange colors, respectively. A susceptible node contacts an infected node through links (1-simplex) and gets infected with probability ${\lambda }_a$ (${\lambda }_b$ ) and ${\lambda }_{ba}$ (${\lambda }_{ab}$ ) per time step via each link, and recovers with probability ${\mu }_a$ (${\mu }_b$ ). (c) (d) In the simplicial contagion model (SCM) of order $D=2$ , a susceptible node contacts one or more infected nodes through links (1-simplex) and gets infected with probability ${\lambda }_a$ (${\lambda }_b$ ) per time step via each link. Additionally, it can acquire infection from 2-faces with probability ${\lambda }_a^{\Delta }$ (${\lambda }_b^{\Delta }$ ), and recovers with probability ${\mu }_a$ (${\mu }_b$ ).

Consider two homogeneous network layers A and B with no degree-degree correlations between them. Nodes in layer A (B) are characterized by average degrees $\langle k_a\rangle$ ($\langle k_b\rangle$ ), while $\langle k_{ab}\rangle$ and $\langle k_{ba}\rangle$ represent the average inter-layer degrees for layers A and B, respectively. $\langle k_{a,\Delta }\rangle$ and $\langle k_{b,\Delta }\rangle$ represent the intra-layer average simplicial degrees (of order 2) for layers A and B, respectively. A susceptible node contacts an infected node through links (1-simplex) and gets infected with probability ${\lambda }_a$ (${\lambda }_b$ ) and ${\lambda }_{ba}$ (${\lambda }_{ab}$ ) per time step via each link, and recovers with probability ${\mu }_a$ (${\mu }_b$ ). In the simplicial contagion model (SCM) of order $D=2$ , a susceptible node contacts one or more infected nodes through links (1-simplex) and gets infected with probability ${\lambda }_a$ (${\lambda }_b$ ) per time step via each link. Additionally, it can acquire infection from 2-faces with probability ${\lambda }_a^{\Delta }$ (${\lambda }_b^{\Delta }$ ), and recovers with probability ${\mu }_a$ (${\mu }_b$ ).

Let ${\rho }_A^t$ and ${\rho }_B^t$ denote the fractions of infected nodes in layers A and B at time $t$ . The evolutionary dynamics can be described by the following system of differential equations:

$\begin{array}{l}\frac{\text{d}{\rho }_A^t}{\text{d}t}=-{\mu }_a{\rho }_A^t+{\lambda }_a\langle k_a\rangle {\rho }_A^t\left(1-{\rho }_A^t\right)+{\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle {\left({\rho }_A^t\right)}^2\left(1-{\rho }_A^t\right)\\ +{\lambda }_{ba}\langle k_{ba}\rangle {\rho }_B^t\left(1-{\rho }_A^t\right)\\ \frac{\text{d}{\rho }_B^t}{\text{d}t}=-{\mu }_b{\rho }_B^t+{\lambda }_b\langle k_b\rangle {\rho }_B^t\left(1-{\rho }_B^t\right)+{\lambda }_b^{\Delta }\langle k_{b,\Delta }\rangle {\left({\rho }_B^t\right)}^2\left(1-{\rho }_B^t\right)\\ +{\lambda }_{ab}\langle k_{ab}\rangle {\rho }_A^t\left(1-{\rho }_B^t\right)\end{array}$ (1)

In layer A, the first term on the right-hand side represents the recovery probability of infected nodes at time $t$ ; the second term corresponds to the infection probability of susceptible nodes through pairwise interactions with infected neighbors; the third term denotes infection through triangular (higher-order) interactions; and the fourth term represents cross-layer infection from infected neighbors in the opposite layer through pairwise interactions. A similar interpretation applies to layer B.

After a transient period, the dynamical system evolves toward a stationary state. To obtain non-trivial stationary solutions, we set $\frac{\text{d}{\rho }_A^t}{\text{d}t}=0$ and $\frac{\text{d}{\rho }_B^t}{\text{d}t}=0$ , yielding:

$\begin{array}{l}0=-{\mu }_a{\rho }_A^t+{\lambda }_a\langle k_a\rangle {\rho }_A^t\left(1-{\rho }_A^t\right)+{\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle {\left({\rho }_A^t\right)}^2\left(1-{\rho }_A^t\right)\\ +{\lambda }_{ba}\langle k_{ba}\rangle {\rho }_B^t\left(1-{\rho }_A^t\right)\\ 0=-{\mu }_b{\rho }_B^t+{\lambda }_b\langle k_b\rangle {\rho }_B^t\left(1-{\rho }_B^t\right)+{\lambda }_b^{\Delta }\langle k_{b,\Delta }\rangle {\left({\rho }_B^t\right)}^2\left(1-{\rho }_B^t\right)\\ +{\lambda }_{ab}\langle k_{ab}\rangle {\rho }_A^t\left(1-{\rho }_B^t\right)\end{array}$ (2)

To facilitate the analysis of non-trivial stationary solutions, we impose symmetry conditions by setting identical parameters for both layers and their interconnections: $\langle k_a\rangle =\langle k_b\rangle$ , $\langle k_{a,\Delta }\rangle =\langle k_{b,\Delta }\rangle$ , $\langle k_{ab}\rangle =\langle k_{ba}\rangle$ , ${\lambda }_a={\lambda }_b$ , ${\lambda }_{ab}={\lambda }_{ba}$ , ${\lambda }_a^{\text{Δ}}={\lambda }_b^{\text{Δ}}$ , along with identical initial infection fractions and recovery rates. This simplification reduces the two-layer network dynamics to an effective single-layer formulation.

Under these symmetric conditions, the system becomes:

$\begin{array}{l}0=-{\mu }_a{\rho }_A^t+{\lambda }_a\langle k_a\rangle {\rho }_A^t\left(1-{\rho }_A^t\right)+{\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle {\left({\rho }_A^t\right)}^2\left(1-{\rho }_A^t\right)\\ +{\lambda }_{ab}\langle k_{ab}\rangle {\rho }_B^t\left(1-{\rho }_A^t\right)\\ 0=-{\mu }_a{\rho }_B^t+{\lambda }_a\langle k_a\rangle {\rho }_B^t\left(1-{\rho }_B^t\right)+{\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle {\left({\rho }_B^t\right)}^2\left(1-{\rho }_B^t\right)\\ +{\lambda }_{ab}\langle k_{ab}\rangle {\rho }_A^t\left(1-{\rho }_B^t\right)\end{array}$ (3)

At equilibrium, we have ${\rho }_A={\rho }_B={\rho }^{\text{*}}$ , leading to:

$\begin{array}{l}0=\left({\lambda }_a\langle k_a\rangle +{\lambda }_{ab}\langle k_{ab}\rangle -{\mu }_a\right){\rho }^{*}+\left({\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle -{\lambda }_a\langle k_a\rangle -{\lambda }_{ab}\langle k_{ab}\rangle \right){\left({\rho }^{*}\right)}^2\\ -{\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle {\left({\rho }^{*}\right)}^3\end{array}$ (4)

This can be factorized as:

$\begin{array}{l}0={\rho }^{*}[\left({\lambda }_a\langle k_a\rangle +{\lambda }_{ab}\langle k_{ab}\rangle -{\mu }_a\right)+\left({\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle -{\lambda }_a\langle k_a\rangle -{\lambda }_{ab}\langle k_{ab}\rangle \right){\rho }^{*}\\ -{\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle {\left({\rho }^{*}\right)}^2]\end{array}$ (5)

Thus, we obtain the trivial solution ${\rho }^{\text{*}}=0$ , along with two additional solutions:

${\rho }_{\pm }^{\text{*}}=\frac{-b\pm \sqrt{b^2-4ac}}{2a}$ (6)

where

$a={\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle$ (7)

$b={\lambda }_a\langle k_a\rangle +{\lambda }_{ab}\langle k_{ab}\rangle -{\lambda }_a^{\Delta }\langle k_{a,\Delta }\rangle$ (8)

$c={\mu }_a-{\lambda }_a\langle k_a\rangle -{\lambda }_{ab}\langle k_{ab}\rangle$ (9)

Consequently, the steady-state equation $\frac{\text{d}{\rho }_A^t}{\text{d}t}=0$ admits up to three solutions within the physically meaningful range $\rho \in \left[0,1\right]$ . The solution ${\rho }^{\text{*}}=0$ corresponds to the absorbing disease-free state where all individuals recover, and the epidemic vanishes. However, a careful stability analysis of this state and the two additional solutions ${\rho }_{+}^{\text{*}}$ and ${\rho }_{-}^{\text{*}}$ is required to fully characterize the system’s phase diagram.

Therefore, through theoretical analysis of symmetric scenarios, we establish a theoretical baseline for the potential bistable mechanisms in the system; subsequent numerical simulations build upon this foundation to further investigate the specific effects of asymmetric conditions—such as node sizes, triangle distribution, and coupling strengths—on propagation dynamics, thereby bridging the gap between theoretical benchmarks and real-world complexity.

3. Simulation

In this section, we present numerical simulations to validate our theoretical predictions.

In Figure 2, we simulate the temporal evolution of the infected node density under symmetric node populations (equal sizes in both layers) and varying initial conditions. When the initial infection density is set to 0.1, the system eventually converges to a disease-free state with zero infection density. In contrast, when the initial infection density exceeds 0.2, the system reaches a stable endemic state. Under these symmetric conditions, the system exhibits a clear bistability: depending on the initial condition, it converges to either the disease-free state or an endemic state with identical steady-state infection densities in both layers.

Figure 2. A randomly synthesized simplicial complex (SC) with dimension $D=2$ (RSC) was generated with the following parameters: both the upper and lower layers contain 412 triangles, with network sizes $N_1=N_2=1000$ . The average intra-layer degrees are $\langle k_a\rangle =\langle k_b\rangle =3.4$ , and the average inter-layer degrees are $\langle k_{ab}\rangle =\langle k_{ba}\rangle =1$ . The average degrees related to triangular interactions are $\langle k_{a,\Delta }\rangle =\langle k_{b,\Delta }\rangle =1.2$ . The transmission rates are set as ${\lambda }_a={\lambda }_b=0.05$ for intra-layer pairwise contacts, ${\lambda }_{ab}={\lambda }_{ba}=0.05$ for inter-layer pairwise contacts, and ${\lambda }_a^{\Delta }={\lambda }_b^{\Delta }=0.8$ for higher-order (2-simplex) interactions. The recovery rate is ${\mu }_a={\mu }_b=0.3$ for both layers. The left panel shows the temporal evolution of the infected density $I_a$ under different initial conditions, while the right panel shows the corresponding dynamics for $I_b$ . The system exhibits bistability, as the final state depends on the initial infection density.

In Figure 3, we investigate the case with asymmetric node populations between the layers while keeping the number of triangles identical to that in Figure 2. The system still reaches a steady state after a transient period, but the bistability observed in the symmetric case disappears. The infection density in layer A (the larger layer) remains relatively unchanged compared to the symmetric scenario. However, layer B (the smaller layer) exhibits a significantly higher infection density than layer A. Furthermore, the total system-wide infection density in Figure 3 is higher than that in Figure 2.

These results demonstrate that, under identical network parameters, the layer

Figure 3. A randomly synthesized simplicial complex model (SCM) with dimension $D=2$ (RSC) was constructed with the following parameters: both the upper and lower layers contain 412 triangles, with node populations set to $N_1=1000$ and $N_2=500$ . The average intra-layer degrees are $\langle k_a\rangle =\langle k_b\rangle =3.4$ , while the average inter-layer degrees are $\langle k_{ab}\rangle =\langle k_{ba}\rangle =1$ . The average degrees associated with triangular (2-simplex) interactions are $\langle k_{a,\Delta }\rangle =\langle k_{b,\Delta }\rangle =1.2$ . The transmission rates are defined as ${\lambda }_a={\lambda }_b=0.05$ for intra-layer pairwise contacts, ${\lambda }_{ab}={\lambda }_{ba}=0.05$ for inter-layer pairwise contacts, and ${\lambda }_a^{\Delta }={\lambda }_b^{\Delta }=0.8$ for higher-order simplex interactions. The recovery rate is ${\mu }_a={\mu }_b=0.3$ for both layers. The left panel displays the temporal evolution of the infected density $I_A$ under different initial conditions, and the right panel shows the corresponding dynamics for $I_B$ . The system ultimately converges to a steady state.

with the smaller node population exhibits stronger disease transmissibility and a higher endemic infection level than the larger layer. Under identical network and transmission parameters, as shown in Figure 4, the infection densities in both upper and lower layers increase with the number of triangles in the network. In the same system with different node populations in the two layers, when the number of triangles is identical, their proportions in the respective networks differ ($\langle k_{a,\Delta }\rangle \ne \langle k_{b,\Delta }\rangle$ ). Specifically, the network with fewer nodes exhibits a higher triangle density, and the infection density $I_B$ is greater than $I_A$ . Therefore, a larger number of triangles facilitates disease transmission and diffusion, and the network layer with a higher triangle density experiences a greater infection density.

In Figure 5, under identical network and transmission parameters, panels (a)-(c) illustrate the system’s spreading behavior with triangles placed in different network locations.

In case (a), where both upper and lower layers contain triangles, the system achieves bistability at a steady state.

Figure 4. A randomly synthesized simplicial complex model (SCM) with dimension $D=2$ (RSC) was generated with the following fixed parameters: $N_1=1000$ , $N_2=500$ , ${\lambda }_a={\lambda }_b=0.06$ , ${\lambda }_{ab}={\lambda }_{ba}=0.06$ , ${\lambda }_a^{\text{Δ}}={\lambda }_b^{\text{Δ}}=0.8$ , and ${\mu }_a={\mu }_b=0.225$ . The number of triangles in the complex was systematically varied across simulations. The left panel shows the steady-state infection density $I_A$ as a function of the number of triangles, while the right panel shows the corresponding steady-state infection density $I_B$ .

In case (b), where only the lower layer contains triangles (with no triangles in the upper layer), the system also reaches bistability. However, compared to case (a), the steady-state infection densities in both layers A and B are reduced.

In case (c), where only the upper layer contains triangles (with no triangles in the lower layer), the infection density drops to zero at steady state.

In case (d), where the total number of triangles is distributed equally between layers (52 triangles in each layer, totaling 104 triangles), the infection density also converges to zero.

These results demonstrate that: 1) when triangles are present in both layers, networks with more triangles exhibit faster disease transmission and higher infection densities than those with fewer triangles; 2) when triangles are concentrated in specific layers, placing triangles exclusively in layer B leads to stronger transmission and higher steady-state infection densities than placing them exclusively in layer A. Consequently, triangles located in layers with smaller node populations enhance disease transmission more effectively.

In Figure 6, it is observed that the infection densities in both the upper and lower layers increase with the rise of ${\lambda }_{a,\text{Δ}}$ and ${\lambda }_{b,\text{Δ}}$ . When ${\lambda }_{b,\text{Δ}}$ is fixed and ${\lambda }_{a,\text{Δ}}$ increases, or when ${\lambda }_{a,\text{Δ}}$ is fixed and ${\lambda }_{b,\text{Δ}}$ increases, a distinct asymmetric effect is revealed. In the upper layer, the increase in infection density is more pronounced when ${\lambda }_{b,\text{Δ}}$ increases with fixed ${\lambda }_{a,\text{Δ}}$ , compared to the case where ${\lambda }_{a,\text{Δ}}$ increases with fixed ${\lambda }_{b,\text{Δ}}$ . Conversely, in the lower layer, the infection density increases more significantly when ${\lambda }_{a,\text{Δ}}$ rises with fixed ${\lambda }_{b,\text{Δ}}$ , compared to increasing ${\lambda }_{b,\text{Δ}}$ with fixed ${\lambda }_{a,\text{Δ}}$ .

(a)

(b)

(c)

(d)

Figure 5. A randomly synthesized simplicial complex model (SCM) with dimension $D=2$ (RSC) was generated using the following parameters: $N_1=1000$ , $N_2=500$ , ${\lambda }_a={\lambda }_b=0.06$ , ${\lambda }_{ab}={\lambda }_{ba}=0.06$ , ${\lambda }_{a,\Delta }={\lambda }_{b,\Delta }=0.8$ , and ${\mu }_a={\mu }_b=0.225$ . The distribution of triangles across layers was configured as follows: (a) 104 triangles in both the upper and lower layers; (b) No triangles in the upper layer and 104 triangles in the lower layer; (c) 104 triangles in the upper layer and no triangles in the lower layer; (d) A total of 104 triangles evenly split, with 52 triangles in each layer. The left panel illustrates the temporal evolution of the infected density $I_A$ under different initial conditions, while the right panel shows the corresponding dynamics for $I_B$ . In all configurations, the system eventually converges to a steady state.

In Figure 7, two asymmetric higher-order configurations are examined:

• In case (a), where higher-order interactions are absent in the upper layer but present in the lower layer, both $I_A$ and $I_B$ increase with ${\lambda }_{b,\text{Δ}}$ , and the increase in $I_B$ is substantially greater than that in $I_A$ .

Figure 6. A randomly synthesized simplicial complex model (SCM) with dimension $D=2$ (RSC) was generated with the following parameters: $N_1=1000$ , $N_2=500$ , intra-layer degrees $\langle k_a\rangle =\langle k_b\rangle =8$ , inter-layer degrees $\langle k_{ab}\rangle =\langle k_{ba}\rangle =5$ , higher-order degrees $\langle k_{a,\Delta }\rangle =3$ , $\langle k_{b,\Delta }\rangle =5$ , infection rates ${\lambda }_a={\lambda }_b=0.03$ (intra-layer), ${\lambda }_{ab}={\lambda }_{ba}=0.01$ (inter-layer), and recovery rates ${\mu }_a={\mu }_b=0.2$ . The left panel shows the phase diagram of the steady-state infection density $I_A$ in the upper layer under varying combinations of higher-order infection rates ${\lambda }_a^{\text{Δ}}$ and ${\lambda }_b^{\text{Δ}}$ , while the right panel displays the corresponding phase diagram for the steady-state infection density $I_B$ in the lower layer.

(a)

(b)

Figure 7. A randomly synthesized simplicial complex model (SCM) with dimension $D=2$ (RSC) was generated with the following parameters: $N_1=1000$ , $N_2=500$ , ${\lambda }_a={\lambda }_b=0.08$ , ${\mu }_a={\mu }_b=0.3$ . Two configurations with higher-order structures are examined: (a): 410 triangles in the upper layer, showing the variation of steady-state infection densities $I_A$ and $I_B$ with respect to ${\lambda }_{b,\Delta }$ (with ${\lambda }_{a,\Delta }=0$ ). (b): 410 triangles in the upper layer, showing the variation of steady-state infection densities $I_A$ and $I_B$ with respect to ${\lambda }_{a,\Delta }$ (with ${\lambda }_{b,\Delta }=0$ ).

Figure 8. A randomly synthesized simplicial complex model (SCM) with dimension $D=2$ (RSC) was generated with the following parameters: both layers contain 104 triangles, with network sizes $N_1=1000$ and $N_2=500$ . The average degrees are $\langle k_a\rangle =1.5$ , $\langle k_b\rangle =2.3$ for intra-layer connections, while the inter-layer coupling strength $\langle k_{ab}\rangle =\langle k_{ba}\rangle$ is varied systematically. The higher-order degrees are $\langle k_{a,\Delta }\rangle =0.31$ and $\langle k_{b,\Delta }\rangle =0.62$ . The transmission rates are set as ${\lambda }_a={\lambda }_b=0.06$ for intra-layer spreading, ${\lambda }_{ab}={\lambda }_{ba}=0.06$ for inter-layer spreading, and ${\lambda }_{a,\Delta }={\lambda }_{b,\Delta }=0.5$ for higher-order interactions. The recovery rate is ${\mu }_a={\mu }_b=0.3$ . The figure shows how the steady-state infection densities $I_a$ and $I_b$ vary with different values of the inter-layer coupling strength $\langle k_{ab}\rangle$ .

• In case (b), where higher-order interactions exist only in the upper layer, both infection densities grow with ${\lambda }_{a,\text{Δ}}$ , but the increase in $I_A$ is markedly larger than that in $I_B$ .

These results collectively demonstrate that variation in ${\lambda }_{a,\text{Δ}}$ exerts a stronger influence on the infection density of the upper layer, while variation in ${\lambda }_{b,\text{Δ}}$ more dominantly affects the lower layer.

In Figure 8, we examine the influence of inter-layer coupling strength on the infection densities $I_A$ and $I_B$ . The results demonstrate that as the inter-layer coupling degree $\langle k_{ab}\rangle$ (equivalent to $\langle k_{ba}\rangle$ ) increases, both infection densities $I_A$ and $I_B$ rise. However, the effect of inter-layer coupling is more pronounced on $I_B$ , indicating that the layer with the smaller node population (the lower layer) is more significantly influenced by changes in inter-layer connectivity.

4. Conclusions

Based on theoretical analysis and numerical simulations, this study systematically investigates the impact of higher-order structures on spreading dynamics in two-layer networks. The main conclusions are as follows:

• Higher-order structures enhance spreading capability: The presence of triangular interactions significantly increases the infection density within the system and can induce bistability phenomena.

• Node population size influences spreading intensity: Under identical parameters, the network layer with a smaller node population is more susceptible to infection, exhibiting greater sensitivity to the spread.

• The distribution of triangles affects spreading pathways: When triangles are concentrated in the layer with fewer nodes, that layer experiences a higher infection density and more intense propagation.

• Inter-layer coupling strength exhibits an asymmetric impact on spreading: Stronger inter-layer connections lead to higher infection densities, with a more pronounced effect on the layer possessing the smaller node population.

• Higher-order infection rates more directly affect their respective layers: The parameter ${\lambda }_{a,\text{Δ}}$ has a greater influence on the infection density of the upper layer, while ${\lambda }_{b,\text{Δ}}$ predominantly affects the lower layer.

These findings reveal the synergistic effects between higher-order structures and multilayer coupling in contagion processes, providing a new theoretical basis for understanding and controlling complex spreading phenomena in the real world.

5. Discussion

While this study has advanced the understanding of higher-order spreading dynamics in two-layer networks, several aspects warrant further investigation. First, our model assumes homogeneous network structures, leaving the impact of degree heterogeneity, such as that in scale-free networks, on higher-order contagion unexplored. Future work could integrate heterogeneous network topologies to deepen the mechanistic understanding.

Second, the current framework does not incorporate dynamic feedback mechanisms, such as behavioral adaptations or immunization strategies, which play critical roles in real-world spreading processes. Incorporating such adaptive dynamics would enhance the model’s practical relevance.

Furthermore, while our study focuses on the SIS model, extending it to more realistic frameworks like SIR or SEIR could provide broader insights into different types of contagion phenomena.

On the practical front, this research offers implications for areas such as epidemic control and public opinion management. For instance, in multi-platform social systems, identifying and intervening in higher-order structures could serve as an effective strategy for curbing the spread of misinformation or disease. Similarly, regulating inter-layer coupling strength may offer a viable approach for managing cross-platform propagation.

Promising future directions include exploring the effects of even higher-order structures, investigating higher-order contagion in temporal networks, and developing control algorithms specifically designed for systems with higher-order interactions. We believe that with the continued development of higher-order network theory, the understanding of spreading behaviors in multilayer interconnected systems will become increasingly thorough and comprehensive.

Funding

This work was supported by National Natural Science Foundation of China (Grants No. 12305047).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Pastor-Satorras, R., Castellano, C., Van Mieghem, P. and Vespignani, A. (2015) Epidemic Processes in Complex Networks. Reviews of Modern Physics, 87, 925-979.[CrossRef]
[2]	Keeling, M.J. and Eames, K.T.D. (2005) Networks and Epidemic Models. Journal of The Royal Society Interface, 2, 295-307.[CrossRef] [PubMed]
[3]	Centola, D. and Macy, M. (2007) Complex Contagions and the Weakness of Long Ties. American Journal of Sociology, 113, 702-734.[CrossRef]
[4]	Bick, C., Gross, E., Harrington, H.A. and Schaub, M.T. (2023) What Are Higher-Order Networks? SIAM Review, 65, 686-731.[CrossRef]
[5]	Bianconi, G. (2021) Higher-Order Networks. Cambridge University Press.[CrossRef]
[6]	Feng, Y., You, H., Zhang, Z., Ji, R. and Gao, Y. (2019) Hypergraph Neural Networks. Proceedings of the AAAI Conference on Artificial Intelligence, 33, 3558-3565.[CrossRef]
[7]	Battiston, F., Cencetti, G., Iacopini, I., Latora, V., Lucas, M., Patania, A., et al. (2020) Networks Beyond Pairwise Interactions: Structure and Dynamics. Physics Reports, 874, 1-92.[CrossRef]
[8]	Iacopini, I., Petri, G., Barrat, A. and Latora, V. (2019) Simplicial Models of Social Contagion. Nature Communications, 10, Article No. 2485.[CrossRef] [PubMed]
[9]	Wang, D., Zhao, Y., Luo, J. and Leng, H. (2021) Simplicial SIRS Epidemic Models with Nonlinear Incidence Rates. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31, Article ID: 053112.[CrossRef] [PubMed]
[10]	Li, W., Xue, X., Pan, L., Lin, T. and Wang, W. (2022) Competing Spreading Dynamics in Simplicial Complex. Applied Mathematics and Computation, 412, Article ID: 126595.[CrossRef]
[11]	Jhun, B., Jo, M. and Kahng, B. (2019) Simplicial SIS Model in Scale-Free Uniform Hypergraph. Journal of Statistical Mechanics: Theory and Experiment, 2019, Article ID: 123207.[CrossRef]
[12]	de Arruda, G.F., Petri, G. and Moreno, Y. (2020) Social Contagion Models on Hypergraphs. Physical Review Research, 2, Article ID: 023032.[CrossRef]
[13]	Matamalas, J.T., Gómez, S. and Arenas, A. (2020) Abrupt Phase Transition of Epidemic Spreading in Simplicial Complexes. Physical Review Research, 2, Article ID: 012049.[CrossRef]
[14]	Kivelä, M., Arenas, A., Barthelemy, M., Gleeson, J.P., Moreno, Y. and Porter, M.A. (2014) Multilayer Networks. Journal of Complex Networks, 2, 203-271.[CrossRef]