Introduction

Simplicial complex is a set composed of simplexes, such as nodes, edges, triangles, and tetrahedrons1. In the real world, many systems can be effectively represented as complex networks consisting of nodes and edges2,3. Within these networks, higher-order structures like simplicial complexes are frequently observed4,5,6,7. For example, a feed-forward loop is a significant network structure in transcriptional regulation networks8,9. Additionally, closely connected structures, such as triangles, are commonly found in networks of scientific cooperation10 and social communication6,7, reflecting the dynamics of teamwork and intra-group communication. Simplicial complexes can also be identified in diverse domains, including brain networks11, email networks6, and human contact networks6. Spreading dynamics in complex networks is one of the most active research areas in network science. Diseases and information are usually assumed to spread through connections between nodes12,13,14,15,16,17. However, empirical evidence reveals that pairwise interactions between nodes cannot capture group interactions in real networks4. For instance, threshold models primarily consider the impact of neighbors without considering the group interaction within the neighborhood17,18,19,20. Therefore, network scientists have increasingly focused on the impact of higher-order interactions21,22,23,24,25. For example, the higher-order interactions of 2-simplices can lead to a phase transition from continuous to discontinuous in social contagion26. In competitive spreading processes involving two types of epidemics, an increase in the infection probability within higher-order structures can result in alternating dominance between the two epidemic types27. Furthermore, in temporal higher-order networks, the temporal properties of the underlying network structure constrain the impact of higher-order epidemic spreading28. In multi-layer networks, resource layers or social layers with 2-simplex structures can effectively suppress epidemic spreading and outbreaks29,30. Finally, the coupling between signals defined on nodes and edges results in explosive topological synchronization31.

In classical spreading dynamic models such as epidemic spreading models, the infection (recovery) process is typically assumed to be a Markovian process for simplicity. In other words, an infection (recovery) event occurs at a constant rate, leading to inter-event times that follow an exponential distribution32,33,34,35. However, studies have shown that many real-world processes are non-Markovian36,37,38,39,40,41,42,43, and the inter-event time does not follow an exponential distribution. For instance, empirical evidence indicates that the inter-event time of infection (recovery) events in epidemic spreading processes does not follow an exponential distribution but a curve distribution with obvious peaks44,45,46. In other words, individuals become infected (recovered) within a concentrated period of time. In innovation diffusion processes, the inter-event time for adopting a new product follows a bell-shaped curve distribution, which is a non-Markovian process47. The significance of non-Markovian processes has garnered increasing attention48,49,50,51,52, including efforts to comprehend their impact on epidemic spreading thresholds53, investigate equivalences between non-Markovian and Markovian epidemic spreading processes54,55,56, examine how non-Markovian recovery processes influence failure propagation57, and uncover the conditions under which Markovian and non-Markovian models are transient-state equivalent58.

However, our understanding of non-Markovian dynamic processes in networks with higher-order structures remains limited. In this paper, we address this gap through a comparative study of two types of dynamical processes: one with Markovian characteristics and the other exhibiting non-Markovian behavior. We extend the higher-order Markovian social contagion model to a higher-order non-Markovian model, considering scenarios where both the infection and recovery processes are non-Markovian. Subsequently, we develop a theoretical framework to describe the temporal evolution of higher-order non-Markovian propagation processes and find that there is an equivalence between higher-order non-Markovian and higher-order Markovian contagion dynamics. Finally, we investigate the impact of higher-order non-Markovian processes on propagation. Our results demonstrate that non-Markovian recovery can enhance network resilience against large-scale infections or small-scale infections under certain conditions.

Results

Simplicial complex

To improve the clarity and coherence of the article, we first define each parameter as shown in Supplementary Table 1. Figure 1a illustrates that a z-simplex σz is a convex polyhedron consisting of (z + 1) nodes, denoted as σz = {p0, p1,..., pz}, also known as z-dimensional simplex, where z ≥ 01. For instance, a node is a 0-simplex, an edge is a 1-simplex, a triangle is a 2-simplex, and a tetrahedron is a 3-simplex. A z-simplex is called as a higher-dimensional simplex if z ≥ 2, which is also referred to as a higher-order structure. A simplicial complex is a structure composed of simplices, such as points, edges, triangles, and tetrahedra, where any two simplices in the simplicial complex do not intersect unless they share common faces. The dimension of a simplicial complex is z if the largest dimension of its simplices is z. Figure 1b presents an example of a 3-simplicial complex.

Fig. 1: Schematic illustration of model and theory-related concepts.
figure 1

a z-dimensional simplices with z ≤ 3. b A simplicial complex, which consists of z-dimensional simplices with z ≤ 3. c Schematic representation of higher-order non-Markovian infection and recovery processes in simplicial complexes. d Schematic representation of virtual nodes, virtual infected nodes, and virtual edges. There are 3 virtual nodes and 3 virtual edges in a 2-simplex, represented by dashed circles and dashed lines, respectively. A virtual node j(2) in a 2-simplex is in the infected state if its 2 real nodes are both in the infected state. Meanwhile, the virtual edge connecting the virtual node is active. In (c, d), a green (red) circle indicates a node in S (I) state, and a white circle indicates a node in S or I state.

Full size image

Higher-order non-Markovian social contagion model

We propose a higher-order non-Markovian social contagion model in simplicial complexes, where the higher-order effect refers to group interactions through higher-dimensional simplices such as triangles. We consider the susceptible-infected-susceptible (SIS) spreading process. In other words, a node can be in one of two states: susceptible (S) state and infected (I) state. Then, we introduce the definitions of z-dimensional virtual nodes j(z), virtual infected nodes, virtual edges, virtual active edges, and ages of virtual infected nodes. A z-dimensional virtual node j(z) refers to a set consisting of z real nodes in a z-dimensional simplex. For example, in a z-simplex σz = {p0, p1,…, pz}, there are z + 1 z-dimensional virtual node {j(z) ~ σz⧹{pi}∣i = 0, 1, 2...z}. For instance, the first row in Fig. 1d depicts that a 2-simplex contains three virtual nodes consisting of 2 real nodes. A z-dimensional virtual node j(z) is in the infected state if all of its z real nodes are in the infected state, as shown in the second row in Fig. 1d. The age τ of a virtual/real node is the elapsed time since it turned into the current state. For a z-dimensional virtual infected node j(z) such as {j(z) = σz⧹{p0}}, its infection age τ depends on the smallest infection age of its z real infected nodes(i. e., \(\tau =\min ({\tau }_{1},{\tau }_{2},\ldots ,{\tau }_{z})\)). A virtual edge connects the z-dimensional virtual node {j(z) = σz⧹{p0}} and the remaining node p0. When the virtual node becomes infected, the connected virtual edge becomes active as well and its active age is initially set as τ = 0. As shown in Fig. 1c, behavioral information is propagated through simplices such as edges (1-simplices), triangles (2-simplices), and so on. For the infection process, a z-dimensional virtual infected node j(z) propagates behavior information at a z-dimensional infection rate η(z)(τ) to the remaining node in the simplex, where τ is the infection age of j(z). The z-dimensional infection rate η(z)(τ) varies with τ, rendering the infection processes non-Markovian. In the case of z = 1, if a node at the end of an edge is in the I state and has been in that state for a duration of τ, it will infect the remaining node at the other end at the rate η(1)(τ), as illustrated in the first row of Fig. 1c. In the case of z = 2 if two nodes within a triangle are in the I state and their smallest duration in that state is τ, they will infect the remaining node in the triangle at the rate η(2)(τ), as shown in the second row of Fig. 1c. If the remaining node is in the S state, it can be infected with the probability η(2)(τ)dt; Otherwise, nothing happens to it. Information propagation through z-simplices follows a similar way. For the recovery process, an infected node will spontaneously recover at the rate of ωrec(τ), where τ is its duration in the I state at time t, as shown in the last row of Fig. 1c. Note that the propagation through edges (i.e., pairwise interactions) is a simple propagation process54, while propagation through higher-dimensional simplices such as triangles (i.e., group interactions) constitutes a more complex propagation process. Our model takes into account the presence of both simple and complex propagation in the networks. Then, we define the specific details of the rules governing the infection process and the recovery process54. We assume a stochastic infection process between a virtual infected node j(z) and a susceptible node i, where the infection dynamics solely depend on the state of node j(z). Upon becoming infected, node j(z) initiates independent infection attempts towards each of its neighbors, irrespective of their states. This process resembles a renewal process characterized by inter-event times drawn from the distribution \({\psi }_{\inf }^{(z)}(\tau )\). Conceptually, this process can be envisaged as a sequence of activation events occurring at random intervals τ, commencing upon node j(z)’s infection. If such an activation event coincides with the susceptible state of neighbor node i, node i becomes infected (refer to Fig. 2). Similarly, upon infection, a node triggers a recovery process, with the duration of recovery, denoted by τ, governed by the distribution ψrec(τ).

Fig. 2: The mechanism of infection propagation from a virtual infected node, denoted as j(z), to a neighboring node i is as follows.
figure 2

A virtual infected node j(z) initiates infection events towards its neighbor i at random time intervals denoted by τ. If node i is already in an infected state, the infection event triggered by virtual infected node j(z) is deemed invalid. In the context of the representation, the red circle symbolizes the occurrence of an infection event triggered by the virtual infected node, with state 0 indicating the susceptible state of the node and state 1 indicating the infected state.

Full size image

Simulation verification

To validate our theory, we carry out Monte Carlo simulations54,59 for the higher-order non-Markovian SIS social contagion dynamics in simplicial complexes. The networks used in the experiments are random simplicial complexes26, where both the 1-simplicial and the 2-simplicial degree distributions of real nodes in the network follow the Poisson distribution. The random simplicial complex (RSC) model produces simplicial complexes of dimension z = 2 as follows. Given a set \({{{{{{{\mathcal{V}}}}}}}}\) of N nodes we connect any two nodes \(i,j\in {{{{{{{\mathcal{V}}}}}}}}\) with probability ζ1 ∈ [0, 1]. Then,for any \(i,j,l\in {{{{{{{\mathcal{V}}}}}}}}\), we add a 2-simplex {i, j, l} with probability ζ2 ∈ [0, 1]. By adjusting ζ1 and ζ2, we can obtain the desired value of 〈k(1)〉 and 〈k(2)〉. 〈k(z)〉 is the z-dimensional degree.

Since the Weibull distribution is a convenient and widely applicable choice for data analysis44,45,46,47, we adopt it as the form for \({\psi }_{\inf }^{(z)}(\tau )\) and ψrec(τ):

$${\psi }_{\inf }^{(z)}(\tau )=\frac{{\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}}{{\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}}{\left[\frac{\tau }{{\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}}\right]}^{{\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}-1}{e}^{-{\left[\frac{\tau }{{\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}}\right]}^{{\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}}}$$
(1)

and

$${\psi }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )=\frac{{\alpha }_{{{{{{{{\rm{R}}}}}}}}}}{{\beta }_{{{{{{{{\rm{R}}}}}}}}}}{\left(\frac{\tau }{{\beta }_{{{{{{{{\rm{R}}}}}}}}}}\right)}^{{\alpha }_{{{{{{{{\rm{R}}}}}}}}}-1}{e}^{-{\left(\frac{\tau }{{\beta }_{{{{{{{{\rm{R}}}}}}}}}}\right)}^{{\alpha }_{{{{{{{{\rm{R}}}}}}}}}}},$$
(2)

where \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}\) (\({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}\)) and αR (βR) are the shape (scale) parameters of the distribution. When βR is smaller (larger), the distribution is more (less) concentrated, indicating a more homogeneous (heterogeneous) recovery time required for each real infected node in the network. When αR ≥ 1, the larger (smaller) αR is, the more (less) concentrated the function distribution is, indicating that the more homogeneous (heterogeneous) the recovery time needed for each real infected node is. When αR ≤ 1, although the distribution is typically ”fat-tailed,” it decays exponentially around τ ≈ 0. Compared with the case of αR = 1, The function distribution decays faster when αR < 1, indicating that most real infected nodes take less time to recover. For simplicity, we only consider simplices of dimension up to z = 2. We simulate the higher-order non-Markovian and higher-order Markovian social contagions in random simplicial complexes with N = 2000, 〈k(1)〉 ≈ 14, and 〈k(2)〉 ≈ 6, where the distributions of \({k}_{i}^{(1)}\) and \({k}_{i}^{(2)}\) follow the Poisson distribution26.

Figure 3 shows the time evolution of the fraction of infected nodes. We start with 2% of nodes randomly chosen as the infection seeds, and their infection ages are initially set to zero. Figure 3a presents the results for z = 1, where information is propagated through edges, representing simple propagation. It can be observed that information propagation is more inhibited when \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}\) is larger. The reason for this is that with a larger \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}\), the distribution of \({\psi }_{\inf }^{(z)}(\tau )\) becomes more concentrated, indicating not only more homogeneity but also longer time needed for information propagation. Figure 3b displays the results for z ≤ 2 with parameters taken in left triangles of Fig. 3a. Information is propagated through edges and triangles, representing both simple and complex propagation processes coexisting. Similarly, it can be observed that the larger \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}\) is, the more inhibited information propagation becomes. However, compared to simple contagion, complex interplay promotes spreading. The dashed lines represent the results predicted by the higher-order non-Markovian homogeneous mean-field theory, as given by Eqs. (13)–(18). It can be seen from Fig. 3a that the theoretical results are in good agreement with the simulation results. However, they begin to deviate when \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=4\), primarily due to nodal-pair dynamic correlations. As shown in Fig. 3b, dynamic correlations are stronger when considering the propagation through 2-simplices, which increases the deviation between theoretical results and simulation results. When the system reaches to the steady state, the theoretical results align well with the simulated ones. We do not go further to study the theory due to the extremely high computational complexity associated with considering nodal-pair dynamical correlations.

Fig. 3: Comparison of time evolution between simulation and theoretical results.
figure 3

a Time evolution of the fraction of infected nodes when only the simplices of dimension z = 1 are considered, namely propagation through edges (simple contagion process only). The blue circles, orange lower triangles, green upper triangles, and red left triangles are the simulation results when \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=1,2,3,4\), respectively. b Time evolution of the fraction of infected nodes when the simplices of dimension z ≤ 2 are considered, namely propagation through edges and triangles (simple and complex contagion processes coexisted). The different symbols represent the results with different \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}\) while keeping \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=4\) fixed. The blue circles, orange lower triangles, green upper triangles, and red left triangles correspond to \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}=1,2,3,4\), respectively. The simulation results are obtained by averaging over 100 realizations, and the dashed lines are the results predicted by the higher-order non-Markovian homogeneous mean-field theory given by Eqs. (13)–(18). Error bars are smaller than the size of the data point and are omitted for clarity. Error bars are defined as standard deviation. Other parameters are N = 2000, 〈k(1)〉≈14, 〈k(2)〉 ≈ 6, I(0) = 0.02, αR = 2.0, βR = 0.45, and \({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=0.9\).

Full size image

To better compare non-Markovian processes with Markovian ones, we firstly define the average infection rate \(\langle {\lambda }^{(z)}\rangle =\frac{1/E({\tau }_{\inf }^{(z)})}{1/E({\tau }_{{{{{{{{\rm{rec}}}}}}}}})}=\frac{E({\tau }_{{{{{{{{\rm{rec}}}}}}}}})}{E({\tau }_{\inf }^{(z)})}\), where \({\tau }_{\inf }^{(z)}\) and τrec are the infection and recovery time needed, respectively. The expected value and variance of \({\psi }_{\inf }^{(z)}(\tau )\) with the form of Weibull distribution are \(E({\tau }_{\inf }^{(z)})={\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}\Gamma (1+\frac{1}{{\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}})\) and \({{{{{{{\rm{Var}}}}}}}}({\tau }_{\inf }^{(z)})={({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)})}^{2}\Gamma (1+\frac{2}{{\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}})-{[E({\tau }_{\inf }^{(z)})]}^{2}\), respectively. Similarly, the expected value and variance of ψrec(τ) with the form of Weibull distribution are \(E({\tau }_{{{{{{{{\rm{rec}}}}}}}}})={\beta }_{{{{{{{{\rm{R}}}}}}}}}\Gamma (1+\frac{1}{{\alpha }_{{{{{{{{\rm{R}}}}}}}}}})\) and \({{{{{{{\rm{Var}}}}}}}}({\tau }_{{{{{{{{\rm{rec}}}}}}}}})={({\beta }_{{{{{{{{\rm{R}}}}}}}}})}^{2}\Gamma (1+\frac{2}{{\alpha }_{{{{{{{{\rm{R}}}}}}}}}})-{[E({\tau }_{{{{{{{{\rm{rec}}}}}}}}})]}^{2}\), respectively. It is worth noting that when the shape parameter equals to 1, the Weibull distribution degenerates into the exponential distribution, indicating that the corresponding process degenerates into a Markovian process. Therefore, the average infection time of \({\psi }_{\inf }^{(z)}(\tau )\) and the average recovery time of ψrec(τ) degenerate into \(E({\tau }_{\inf }^{(z)})=1/{\beta }^{(z)}\) and E(τrec) = 1/μ, respectively. Here, β(z) and μ are the Markovian z-dimensional infection and the Markovian recovery rates, respectively. Both are constants. Hence, the average infection rate corresponds to 〈λ(z)〉 = β(z)/μ, which is the same with the Markovian effective infection rate. For each pair of shape parameter \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}\) and average infection time \(E({\tau }_{\inf }^{(z)})\) which are known, the scale parameter \({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}\) can be obtained as \({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}={(\frac{1}{E({\tau }_{\inf }^{(z)})}\Gamma (1+\frac{1}{{\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(z)}}))}^{-1}\). The scale parameter βR can be obtained in the same way. Figure 4a shows the dependence of stationary \(\tilde{I}\) on the average infection rate 〈λ(1)〉 in the case of z = 1, where E(τrec) = 2.0 is fixed. The dot-dashed lines in Fig. 4a are the results of the higher-order non-Markovian homogeneous mean-field theory given by Eqs. (13)–(18). It can be seen that the theoretical result and the simulation result are in good agreement in the steady state. We can also see that the larger \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}\) is, the more inhibited the information propagation is.

Fig. 4: Equivalence between higher-order non-Markovian and higher-order Markovian propagation.
figure 4

a Dependence of \(\tilde{I}\) on 〈λ(1)〉 in the case of z = 1. b Dependence of \(\tilde{I}\) on \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) in the case of z = 1. c Dependence of \(\tilde{I}\) on 〈λ(1)〉 in the case of z ≤ 2. d Dependence of \(\tilde{I}\) on \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) in the case of z ≤ 2. The blue circles, orange lower triangles and green upper triangles are the simulation results when \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=1,5,10\), respectively. The red left triangles are the simulation results for Markovian propagation, where μ = 0.5 in (a) and (b), and β(2) = 0.315 and μ = 0.5 in (c, d). The dot-dashed lines in (a, c) are the theoretical results obtained by Eqs. (13)–(18). The dashed lines and the dot-dashed lines in (b, d) are the theoretical results obtained by the Markovian mean-field and the Markovian pair approximation methods, respectively. The simulated and theoretical results in the figure reveal an equivalence between higher-order non-Markovian and higher-order Markovian propagation. The simulation results are obtained by averaging over 100 realizations. Error bars are smaller than the size of the data point and are omitted for clarity. Error bars are defined as standard deviation. In (c) and (d), \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(2)}\,\approx \,0.63\) by letting \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}=15.0\) and \({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}\approx 1.425815\). Other parameters are N = 2000, 〈k(1)〉 ≈ 14, 〈k(2)〉 ≈ 6, I(0) = 0.02, and αR = 15.0. The average recovery time E(τrec) is set as 2.0 by letting βR≈2.071118.

Full size image

Obviously, the results for the non-Markovian processes can differ significantly from those for the Markovian case, as demonstrated in Fig. 4a. In the preceding section, we theoretically derived that there is an equivalence between higher-order non-Markovian and higher-order Markovian processes, when the non-Markovian effective infection rate takes the form of Eq. (42). Therefore, we replace the abscissa 〈λ(1)〉 in Fig. 4a with the effective infection rate \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\). Figure 4b shows the dependence of \(\tilde{I}\) on \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) with different \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}\) for z = 1, where \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(2)}\approx 0.63\) is fixed. It can be seen that the symbols and the lines overlap each other, which implies equivalence. Our simulation results agree well with the theoretical results. Figure 4c, d are the results for z ≤ 2. In comparison to Fig. 4a, the results in Fig. 4c reveal that when the group interactions are considered, there is a phase transition from the continuous one to the discontinuous one. Analogous results have been previously observed for the Markovian model26. Furthermore, the equivalence between higher-order non-Markovian propagation and higher-order Markovian propagation is confirmed by the simulation and theoretical results depicted in Fig. 4d. As shown in Fig. 4b, the Markovian mean-field theoretical results are consistent with the simulation results when z = 1. However, the power of its prediction is not as good as that of pair approximation when z ≤ 2, due to the strong nodal-pair dynamic correlations, as shown in Fig. 4d. Figure 5 displays the averaged results depicting how the steady-state behavior changes with variations in 〈k(2)〉. It can be observed that as 〈k(2)〉 increases, the steady-state infection density, and the spreading threshold also increase. Furthermore, an observed transition from continuous to discontinuous phase transition occurs as the average 2-simplicial degree 〈k(2)〉 increases from 0 to 7.

Fig. 5: The effect of the average 2-simplicial degree 〈k(2)〉 on the the steady-state infection density.
figure 5

a Dependence of \(\tilde{I}\) on 〈λ(1)〉 for \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=1.0\). b Dependence of \(\tilde{I}\) on 〈λ(1)〉 for \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=5.0\). c Dependence of \(\tilde{I}\) on 〈λ(1)〉 for \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=10.0\). Other parameters are 〈k(1)〉≈14, E(τrec) = 2.0, I(0) = 0.02, and αR = 15.0. The blue circles, orange lower triangles, green upper triangles, red left triangles, and purple right triangles are the simulation results when 〈k(2)〉 = 0, 2, 4, 6, 7, respectively. Error bars are smaller than the size of the data point and are omitted for clarity. Error bars are defined as standard deviation.

Full size image

The effect of non-Markovian processes on system resilience

System function, such as system resilience is a hot topic in the network science60,61,62,63. It’s very meaningful to investigate the impact of the non-Markovian process on the system resilience since the real process is non-Markovian. In order to simplify simulations, we consider the case of z ≤ 2, where both simple and complex propagation processes exist. We fix \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(2)}\approx 0.63\) to observe the dependence of the fraction of stationary infected nodes \(\tilde{I}\) on \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\), as shown in Fig. 6a. The simulation results are obtained by swabbing \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) up and down in small step, starting with I(0) = 0.02 for \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}=0\). The final state for a value of \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) is used as the initial state of the simulations for the next value of \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\)–which is an adiabatic process. The green upper (orange lower) triangles are the simulation results where both z-dimensional infection and recovery processes are Markovian (non-Markovian) processes. The blue circles are the simulation results where z-dimensional infection processes are Markovian processes, while the recovery processes are non-Markovian. These results reveal that a discontinuous phase transition occurs as \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) adiabatic varies either from 0.00 to 0.15 or from 0.15 to 0.00–the signature of a hysteresis. When moving \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) up (down) to a specific value, \(\tilde{I}\) varies from zero (a large value) to a large value (zero), indicating that the state of the system jumps from the low-infection (high-infection) state to the high-infection (low-infection) state. Within the region of hysteresis, there exist two stable solutions for each \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\).

Fig. 6: The effect of non-Markovian processes on system resilience in the case of z ≤ 2, where \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(2)}=0.63\).
figure 6

a Dependence of \(\tilde{I}\) on \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\). The green upper triangles are the simulation results that both the z-dimensional infection and recovery processes are Markovian processes, where β(2) = 0.315, μ = 0.5. The blue circles are the simulation results that z-dimensional infection processes are Markovian processes while the recovery one are non-Markovian processes, where \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=1.0\), \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}=1.0\), \({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}\approx 3.1927\), αR = 15.0, and βR ≈ 2.071118. The orange lower triangles are the simulation results that both the z-dimensional infection and recovery processes are non-Markovian processes, where \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=10.0\), \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}=15.0\), \({\beta }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}\approx 1.425815\), αR = 15.0, and βR �� 2.071118. The dashed lines and the dot-dashed lines are the theoretical results of Markovian mean-field and pair approximation methods, respectively. b, d The dependence of \({\lambda }_{{{{{{{{\rm{c}}}}}}}}}^{(1)}\) for reaching a high-failure state on the initial value I(0). The blue circles are the simulation results that both the z-dimensional infection and recovery processes are Markovian processes, where β(2) = 0.315 and μ = 0.5. The remaining symbols in (b) are the results that z-dimensional infection processes are Markovian processes while the recovery one are non-Markovian processes, where the lower, upper, left and right triangles correspond to αR = 0.5, 5.0, 10.0, 15.0. The remaining symbols in (c, d) are the results that both z-dimensional infection and recovery processes are non-Markovian processes when αR = 0.5 and αR = 15.0, respectively. Besides, \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(2)}=15.0\) in (c, d). The upper, left, and right triangles in (c) and (d) are the results when \({\alpha }_{{{{{{{{\rm{I}}}}}}}}}^{(1)}=1.0,5.0,10.0\), respectively. The simulation results in the figure are averaged. Error bars are shown in the figure. Error bars are defined as standard deviation. Other parameters are N = 2000, 〈k(1)〉≈14, 〈k(2)〉≈6, and E(τrec) = 2.0.

Full size image

Similar to Fig. 4, the overlap between simulated and theoretical results in Fig. 6 confirms the equivalence between higher-order non-Markovian and higher-order Markovian propagation. It can be seen that when z ≤ 2, the nodal-pair dynamic correlations are so strong that there is a deviation between the higher-order Markovian mean-field theory and the simulation results (green upper triangles), while the higher-order Markovian pair approximation theory tends to provide a more accurate prediction of the simulation results.

We next investigate network resilience, which assesses how networks respond when subjected to attacks or changing external conditions. Resilience is a key metric for characterizing the critical functionality of networks when perturbations occur61,62,63,64. As shown in Fig. 6a, when the effective infection rate \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) increases (decreases) adiabatically to a critical value, the systems jump from the low-infection (high-infection) state to the high-infection (low-infection) state. We can define the critical value \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) as \({\lambda }_{{{{{{{{\rm{c}}}}}}}}}^{(1)}\) to characterize the resilience, representing the point at which the system changes state. The resilience performance shown in Fig. 6a is similar for both the non-Markovian and Markovian models. We then examine resilience for both the non-Markovian and Markovian models, by starting with different fractions of initially infected nodes I(0). Each \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) owns a stationary result \(\tilde{I}\) by starting the system from the same I(0). We search for the critical value of \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\), beyond which the stationary state of system will experience from the low-infection one to the high-infection one. Obviously, a larger value of \({\lambda }_{{{{{{{{\rm{c}}}}}}}}}^{(1)}\) indicates that the network is more resilient to the large-scale infection as it requires a larger value of \({\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}\) to drive the system into a high-infection state. Therefore, a critical value \({\lambda }_{{{{{{{{\rm{c}}}}}}}}}^{(1)}\) is calculated when I(0) is given. Figure 6b, d display the dependence of \({\lambda }_{{{{{{{{\rm{c}}}}}}}}}^{(1)}\) on I(0), where E(τrec) = 2.0 is fixed. Figure 6b demonstrates the simulated results of the case where the z-dimensional infection processes are Markovian while the recovery processes are non-Markovian. Here, z ≤ 2. Comparing to the case where both the infection and recovery processes are Markovian processes (blue circles), we find that non-Markovian recovery tends to enhance the resilience of system to withstand a large-scale infection of I(0) ≥ 0.8, when αR is larger than 1.0 and is increased to the large enough value such as αR = 10.0 (red left triangles) or αR = 15.0 (purple right triangles). This is because when αR > 1.0, the larger the value is, the more concentrated the function distribution of ψrec(τ), which means that the recovery time needed for the infected nodes in the network tend to be more homogeneous and longer. When there are many infected nodes initially, many of them will have to wait for long and close time period to recover. At that time, the system becomes one with only a few infected nodes, which is effectively equivalent to one with small value of I(0) and requires a larger value of \({\lambda }_{{{{{{{{\rm{c}}}}}}}}}^{(1)}\) to evolve into the high-infection state. When αR → + ∞, the distribution of ψrec(τ) is approaching to δ distribution, and analogous results have been observed57. In addition, when αR < 1.0 and is decreased to a small enough value such as αR = 0.5 (orange lower triangles), we find that non-Markovian recovery can improve system resilience to withstand a small-scale infection of I(0) ��� [0, 0.7], as shown in Fig. 6b. The reason is that when αR < 1.0, most real infected nodes recover more quickly compared to when αR = 1.0. Most of nodes who are infected initially will recover so soon that the state of system goes to a low-infection state, which enhances the resilience. The results in Fig. 6b is the case of Markovian infection processes. We also do simulations in the case of that both the infection and recovery processes are non-Markovian processes, as shown in Fig. 6c with αR = 0.5 fixed and Fig. 6d with αR = 15.0 fixed. Similar results to Fig. 6b can be observed.

Conclusions

In this study, we considered the impact of higher-order interactions between nodes on spreading dynamics, and propose a general formalism to study non-Markovian propagation dynamics on higher-order networks. We proposed a higher-order non-Markovian social propagation model in simplicial complexes and developed the theoretical framework. Then, we theoretically deduced that there is an equivalence between the higher-order non-Markovian and higher-order Markovian social propagation processes, which is verified by the simulation results. Results indicated that the higher-order interactions can significantly promote the spreading dynamics. Furthermore, we investigated the effect of higher-order non-Markovian propagation on the system resilience. We found that non-Markovian recovery can boost the system resilience to withstand a large-scale infection or a small-scale infection under different conditions. There are shortcomings in this study as well. For example, due to the strong dynamic correlations between nodes in higher-dimensional simplices, there is a deviation between the simulation and theoretic results. The pursuit of more precise theoretical models poses an intriguing challenge for future research. In summary, we studied the higher-order non-Markovian social contagion in simplicial complexes and investigated the impact of non-Markovian processes on system resilience. These results may provide some guidance for the understanding and control of social propagation in the real world, such as public opinion and rumor propagation.

Methods

Mean-field theory for non-Markovian social contagion in simplicial complexes

Inspired by the recently proposed theoretical framework for non-Markovian epidemic spreading processes55, we develop a mean-field theory for higher-order non-Markovian social contagion in simplicial complexes. In the non-Markovian propagation process, \({\psi }_{\inf }^{(z)}(\tau )\) (ψrec(τ)) is the probability density function of the time interval τ between two consecutive z-dimensional infection events (recovery events). The recovery rate ωrec(τ) and the z-dimensional infection rate η(z)(τ) can be computed as55,65

$${\omega }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )=\frac{{\psi }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )}{{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )}$$
(3)

and

$${\eta }^{(z)}(\tau )=\int\nolimits_{0}^{\tau }{\eta }^{(z)}({\tau }^{{\prime} }){\psi }_{\inf }^{(z)}(\tau -{\tau }^{{\prime} })d{\tau }^{{\prime} }+{\psi }_{\inf }^{(z)}(\tau )$$
(4)

Performing Laplace transform, we have

$${\eta }^{(z)}(\tau )=\frac{1}{2\pi i}\int\nolimits_{\sigma -i\infty }^{\sigma +i\infty }\frac{{\hat{\psi }}_{\inf }^{(z)}(s)}{1-{\hat{\psi }}_{\inf }^{(z)}(s)}{e}^{s\tau }ds,$$
(5)

The survival probability Ψrec(τ) represents the probability that a recovery event will not occur before age τ, which are given by

$${\Psi }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )=\int\nolimits_{\tau }^{+\infty }{\psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{{\prime} })d{\tau }^{{\prime} }$$
(6)

In order to describe the non-Markovian propagation dynamics, we use Ii(τ; t) and Si(τ; t) to represent the probability density function that a node i is in the I state and the S state with age τ at time t, respectively. When the time becomes t + dt after the infinitesimal time interval dt, the ages of Ii and Si increase by dτ, we thus have dt = dτ. The evolution equations of a system can be written as

$$\left(\frac{\partial }{\partial \tau }+\frac{\partial }{\partial t}\right){I}_{i}(\tau ;t)=-{I}_{i}(\tau ;t){\omega }_{{{{{{{{\rm{rec}}}}}}}}}(\tau ),$$
(7)
$$\left(\frac{\partial }{\partial \tau }+\frac{\partial }{\partial t}\right){S}_{i}(\tau ;t)=-{S}_{i}(\tau ;t){\Phi }_{i}(t),$$
(8)
$${I}_{i}(0;t+dt)=\int\nolimits_{0}^{+\infty }{S}_{i}(\tau ;t){\Phi }_{i}(t)d\tau ,$$
(9)
$${S}_{i}(0;t+dt)=\int\nolimits_{0}^{+\infty }{I}_{i}(\tau ;t){\omega }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )d\tau .$$
(10)

Eq.(7) [Eq. (8)] represents the decrease in the probability density of infected (susceptible) node i with age τ at time t, due to the recovery (infection) of the node. Eq. (9) [Eq. (10)] represents the increase in the probability density of node i being in the I (S) state with age τ = 0 at time t + dt, due to the infection (recovery) of the other susceptible (infected) nodes with different ages at time t. There are relations \({I}_{i}(t)=\int\nolimits_{0}^{+\infty }{I}_{i}(\tau ;t)d\tau\) and \({S}_{i}(t)=\int\nolimits_{0}^{+\infty }{S}_{i}(\tau ;t)d\tau\). We use Φi(t) to represent the rate that a susceptible node i is infected at time t, which can be written as

$${\Phi }_{i}(t)=\sum\limits_{z=1}^{N-1}\sum\limits_{{j}^{(z)}\in {Z}^{(z)}(i)}\int\nolimits_{0}^{+\infty }{\eta }^{(z)}({\tau }^{{\prime} }){B}_{{j}^{(z)}}({\tau }^{{\prime} };t)d{\tau }^{{\prime} },$$
(11)

where z is the dimension of simplices and z ∈ [1, N − 1]. Z(z)(i) is the set of virtual infected nodes, all of which consist of z real nodes and are the virtual neighbors of node i. j(z) is an element in the set Z(z)(i). \({B}_{{j}^{(z)}}({\tau }^{{\prime} };t)\) is the probability density function that the virtual node j(z) is in the infected state and its infection age is \({\tau }^{{\prime} }\) at time t. The age \({\tau }^{{\prime} }\) of \({B}_{{j}^{(z)}}({\tau }^{{\prime} };t)\) is determined by the smallest infection age of the real infected nodes of the virtual node j(z). \({B}_{{j}^{(z)}}({\tau }^{{\prime} };t)\) can thus be written as

$${B}_{{j}^{(z)}}({\tau }^{{\prime} };t)= \sum\limits_{m\in {j}^{(z)}}{I}_{m}({\tau }^{{\prime} };t)\prod\limits_{n\in {j}^{(z)}\setminus m}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{I}_{n}({\tau }^{\prime\prime};t)d{\tau }^{\prime\prime}\\ = -\frac{d}{d{\tau }^{{\prime} }}\prod\limits_{n\in {j}^{(z)}}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{I}_{n}({\tau }^{\prime\prime};t)d{\tau }^{\prime\prime}d{\tau }^{{\prime} },$$
(12)

Eqs. (7)–(12) depict the temporal evolution of the higher-order non-Markovian propagation dynamics in simplicial complexes. When the number of z-simplices on a real node is homogeneously distributed, the equations can be simplified as

$$\left(\frac{\partial }{\partial \tau }+\frac{\partial }{\partial t}\right)I(\tau ;t)=-I(\tau ;t){\omega }_{{{{{{{{\rm{rec}}}}}}}}}(\tau ),$$
(13)
$$\left(\frac{\partial }{\partial \tau }+\frac{\partial }{\partial t}\right)S(\tau ;t)=-S(\tau ;t)\,\Phi \,(t),$$
(14)
$$I(0;t+dt)=\int\nolimits_{0}^{+\infty }S(\tau ;t)\,\Phi \,(t)d\tau$$
(15)

and

$$S\,(0;t+dt)=\int\nolimits_{0}^{+\infty }I(\tau ;t){\omega }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )d\tau ,$$
(16)

where

$$\Phi (t)=\sum\limits_{z=1}^{N-1}\langle {k}^{(z)}\rangle \int\nolimits_{0}^{+\infty }{\eta }^{(z)}({\tau }^{{\prime} }){B}^{(z)}({\tau }^{{\prime} };t)d{\tau }^{{\prime} }$$
(17)

and

$${B}^{(z)}({\tau }^{{\prime} };t)=zI({\tau }^{{\prime} };t){\left(\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }I({\tau }^{\prime\prime};t)d{\tau }^{\prime\prime}\right)}^{z-1}.$$
(18)

Here, \({\eta }^{(z)}({\tau }^{{\prime} })\) is depicted by Eq. (5). 〈k(z)〉 represents the average number of z-simplices incident on a real node in the network, namely the average z-simplicial degree of a real node. In the case of z = 1, 〈k(1)〉 is the average degree. The relationship among B(z)(t), \({B}^{(z)}({\tau }^{{\prime} };t)\), and I(t) is \({B}^{(z)}(t)=\int\nolimits_{0}^{+\infty }{B}^{(z)}({\tau }^{{\prime} };t)d{\tau }^{{\prime} }={(I(t))}^{z}\), where B(z)(t) is the probability of virtual node j(z) being in the infected state. We have \({B}^{(1)}({\tau }^{{\prime} };t)=I({\tau }^{{\prime} };t)\) for the reason that a virtual node j(1) is actually a real node in the network.

When considering a homogeneous network, and assuming both the infection and the recovery processes are Markovian processes, the above equations can be simplified into the following homogeneous mean-field theory for the higher-order Markovian propagation dynamics as

$$\frac{dI(t)}{dt}=\sum\limits_{z=1}^{N-1}\langle {k}^{(z)}\rangle {\beta }^{(z)}S(t){\left(I(t)\right)}^{z}-\mu I(t),$$
(19)

where S(t) = 1 − I(t). μ and β(z) are the recovery rate and z-dimensional infection rate respectively, which are constant.

Furthermore, we take into account two-node spatial correlation based on the pairwise approximation (PA). Firstly, we define [UV](t) as the fraction of UV-type edges at time t, where U, V ∈ {S, B(z)}. U and V cannot be the virtual nodes at the same time. We have B(1) = I. If either U or V is a virtual infected node, [UV](t) is the fraction of virtual active edge. For instance, [SB(1)](t) = [SI](t), [SB(2)](t) = [SII](t) and [SB(3)](t) = [SIII](t). Based on the pair approximation66, which considers the evolution of different types of edges, the above equations can be written as

$$\frac{dI(t)}{dt}=\sum\limits_{z=1}^{N-1}\langle {k}^{(z)}\rangle {\beta }^{(z)}S(t)\frac{[S{B}^{(z)}](t)}{S(t)}-\mu I(t),$$
(20)
$$\frac{d[II](t)}{dt}=\, 2\sum\limits_{z=1}^{N-1}\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(z)}\rangle {\beta }^{(z)}[SI](t)\frac{[S{B}^{(z)}](t)}{[SI](t)}\\ +2\sum\limits_{z=1}^{N-1}\left(\langle {k}^{(z)}\rangle -\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(z)}\rangle \right){\beta }^{(z)}[SI](t)\frac{[S{B}^{(z)}](t)}{S(t)}\\ -2\mu [II](t).$$
(21)

We then have

$$\frac{dI(t)}{dt}=\sum\limits_{z=1}^{N-1}\langle {k}^{(z)}\rangle {\beta }^{(z)}[S{B}^{(z)}](t)-\mu I(t),$$
(22)
$$\frac{d[II](t)}{dt}=\, 2\sum\limits_{z=1}^{N-1}\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(z)}\rangle {\beta }^{(z)}[S{B}^{(z)}](t)\\ +2\sum\limits_{z=1}^{N-1}\left(\langle {k}^{(z)}\rangle -\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(z)}\rangle \right){\beta }^{(z)}[SI](t)\frac{\left[S{B}^{(z)}\right](t)}{S(t)}\\ -2\mu [II](t).$$
(23)

where \(\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(z)}\rangle\) represents the average number of z-simplices incident on a real edge in the network, namely the average z-simplicial degree of a real edge. For example, we have \(\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(1)}\rangle =1\) if z = 1 and \(\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(2)}\rangle =(N\langle {k}^{(2)}\rangle )/m\) if z = 2, where N is the network size, and m is the total number of real edges in the network. The first (second) term in Eq. (23) represents the transition that a susceptible node at the end of a SI-type edge is infected by the z-simplex (not) incident to the edge, which increases the number of II-type edges. The third term in Eq. (23) represents the transition that an infected node at one end of II-type edge recovers, which decreases the number of II-type edges. By combining the pair approximation method with the characteristics of the triangular structure67, we obtain the following closure approximation for z ≤ 2:

$$[UVW](t) =\frac{[UV](t)}{U(t)}\frac{[VW](t)}{V(t)}\frac{[WU](t)}{W(t)}\\ =\frac{[UV](t)[VW](t)[WU](t)}{U(t)V(t)W(t)}.$$
(24)

The above evolution equations can then be written as

$$\frac{dI(t)}{dt}=\langle {k}^{(1)}\rangle {\beta }^{(1)}[SI](t)+\langle {k}^{(2)}\rangle {\beta }^{(2)}[SII](t)-\mu I(t),$$
(25)
$$\frac{d[II](t)}{dt}\,=\, 2{\beta }^{(1)}[SI](t)+2\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(2)}\rangle {\beta }^{(2)}[SII](t)\\ +2(\langle {k}^{(1)}\rangle -1){\beta }^{(1)}[SI](t)\frac{[SI](t)}{S(t)}\\ +2(\langle {k}^{(2)}\rangle -\langle {k}_{{{{{{{{\rm{e}}}}}}}}}^{(2)}\rangle ){\beta }^{(2)}[SI](t)\frac{[SII](t)}{S(t)}\\ -2\mu [II](t),$$
(26)

where

$$[SII](t)=\frac{[SI](t)[II](t)[IS](t)}{S(t)I(t)I(t)}.$$
(27)

Equivalence between higher-order non-Markovian and higher-order Markovian spreading dynamics

Next, we investigate the stationary relationship between higher-order non-Markovian and higher-order Markovian spreading dynamics in simplicial complexes. We use \(\tilde{I}\) and \(\tilde{S}\) to represent the probabilities of nodes in the I state and in the S state at the steady state, respectively. When the system reaches a steady state, the higher-order non-Markovian propagation equations can then be written as

$$\frac{d}{d\tau }{\tilde{I}}_{i}(\tau )=-{\tilde{I}}_{i}(\tau ){\omega }_{{{{{{{{\rm{rec}}}}}}}}}(\tau ),$$
(28)
$$\frac{d}{d\tau }{\tilde{S}}_{i}(\tau )=-{\tilde{S}}_{i}(\tau ;t){\tilde{\Phi }}_{i},$$
(29)
$${\tilde{I}}_{i}(0)={\tilde{S}}_{i}(0).$$
(30)

Integrating the both sides of Eqs. (28) and (29), we have

$${\tilde{I}}_{i}(\tau )={\tilde{I}}_{i}(0){\Psi }_{{{{{{{{\rm{rec}}}}}}}}}(\tau ),$$
(31)
$${\tilde{S}}_{i}(\tau )={\tilde{S}}_{i}(0){e}^{-{\tilde{\Phi }}_{i}{\tau }^{{\prime} }}.$$
(32)

Substituting Eq. (31) into the following equation of \({\tilde{B}}_{{j}^{(z)}}({\tau }^{{\prime} })\) in the steady state

$${\tilde{B}}_{{j}^{(z)}}({\tau }^{{\prime} }) =-\frac{d}{d{\tau }^{{\prime} }}\prod\limits_{n\in {j}^{(z)}}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{\tilde{I}}_{n}({\tau }^{\prime\prime})d{\tau }^{\prime\prime}d{\tau }^{{\prime} }\\ =\sum\limits_{m\in {j}^{(z)}}{\tilde{I}}_{m}(\tau )\prod\limits_{n\in {j}^{(z)}\setminus m}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{\tilde{I}}_{n}({\tau }^{\prime\prime})d{\tau }^{\prime\prime},$$
(33)

we can get

$${\tilde{B}}_{{j}^{(z)}}({\tau }^{{\prime} })=z{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{{\prime} }){\left(\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{\prime\prime})d{\tau }^{\prime\prime}\right)}^{z-1}\prod\limits_{n\in {j}^{(z)}}{\tilde{I}}_{n}(0).$$
(34)

When \({\tau }^{{\prime} }=0\), we have

$${\tilde{B}}_{{j}^{(z)}}(0)=z{\left(\frac{1}{{\delta }_{{{{{{{{\rm{eff}}}}}}}}}}\right)}^{z-1}\prod\limits_{n\in {j}^{(z)}}{\tilde{I}}_{n}(0)=z{\delta }_{{{{{{{{\rm{eff}}}}}}}}}\prod\limits_{n\in {j}^{(z)}}{\tilde{I}}_{n},$$
(35)

where we define \({\delta }_{{{{{{{{\rm{eff}}}}}}}}}=1/\int\nolimits_{0}^{+\infty }{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )d\tau\). Substituting Eq. (35) into Eq. (34), we have

$${\tilde{B}}_{{j}^{(z)}}({\tau }^{{\prime} })={\tilde{B}}_{{j}^{(z)}}(0){\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{{\prime} }){\left({\delta }_{{{{{{{{\rm{eff}}}}}}}}}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{\prime\prime})d{\tau }^{\prime\prime}\right)}^{z-1}.$$
(36)

Since \({\tilde{\Phi }}_{i}\) is a constant in the steady state, Eqs. (31) and (32) are integrated to obtain

$${\tilde{I}}_{i}=\int\nolimits_{0}^{\infty }{\tilde{I}}_{i}(0){\Psi }_{{{{{{{{\rm{rec}}}}}}}}}(\tau )d\tau =\frac{{\tilde{I}}_{i}(0)}{{\delta }_{{{{{{{{\rm{eff}}}}}}}}}},$$
(37)
$${\tilde{S}}_{i}=\int\nolimits_{0}^{\infty }{\tilde{S}}_{i}(0){e}^{-{\tilde{\Phi }}_{i}\tau }\tau =\frac{{\tilde{S}}_{i}(0)}{{\tilde{\Phi }}_{i}}.$$
(38)

Combining Eqs. (37) and (38), we get

$${\tilde{I}}_{i}={\tilde{S}}_{i}{\tilde{\Phi }}_{i}/{\delta }_{{{{{{{{\rm{eff}}}}}}}}}.$$
(39)

Combining Eqs. (39), (36), and (11), we get

$${\tilde{I}}_{i}= {\tilde{S}}_{i}\sum\limits_{z=1}^{N-1}\sum\limits_{{j}^{(z)}\in {Z}^{(z)}(i)}\int\nolimits_{0}^{+\infty }{\eta }^{(z)}({\tau }^{{\prime} })\left(z{\delta }_{{{{{{{{\rm{eff}}}}}}}}}\prod\limits_{n\in {j}^{(z)}}{\tilde{I}}_{n}\right){\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{{\prime} })\\ {\left[{\delta }_{{{{{{{{\rm{eff}}}}}}}}}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{\prime\prime})d{\tau }^{\prime\prime}\right]}^{z-1}d{\tau }^{{\prime} }/{\delta }_{{{{{{{{\rm{eff}}}}}}}}}\\ = {\tilde{S}}_{i}\sum\limits_{z=1}^{N-1}\sum\limits_{{j}^{(z)}\in {Z}^{(z)}(i)}\left(\prod\limits_{n\in {j}^{(z)}}{\tilde{I}}_{n}\right)\int\nolimits_{0}^{+\infty }z{\eta }^{(z)}({\tau }^{{\prime} }){\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{{\prime} })\\ {\left[{\delta }_{{{{{{{{\rm{eff}}}}}}}}}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{\prime\prime})d{\tau }^{\prime\prime}\right]}^{z-1}d{\tau }^{{\prime} }.$$
(40)

Then Eq. (40) reduces to

$${\tilde{I}}_{i}=\sum\limits_{z=1}^{N-1}{\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(z)}{\tilde{S}}_{i}\sum\limits_{{j}^{(z)}\in {Z}^{(z)}(i)}\,\prod\limits_{n\in {j}^{(z)}}{\tilde{I}}_{n},$$
(41)

where we define

$${\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(z)}=\int\nolimits_{0}^{+\infty }z{\eta }^{(z)}({\tau }^{{\prime} }){\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{{\prime} }){\left[{\delta }_{{{{{{{{\rm{eff}}}}}}}}}\int\nolimits_{{\tau }^{{\prime} }}^{+\infty }{\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{\prime\prime})d{\tau }^{\prime\prime}\right]}^{z-1}d{\tau }^{{\prime} }$$
(42)

as the effective infection rate. It can be seen that the form of Eq. (40) is the same as that of the Markovian equation Eq. (41), which indicates that there is an equivalence between the higher-order non-Markovian and higher-order Markovian propagation dynamics at the steady state in simplicial complexes.

In the case of z = 1, the above equations are simplified as Eq. (43), which are the effective infection rate of the classic SIS non-Markovian spreading processes55:

$${\lambda }_{{{{{{{{\rm{eff}}}}}}}}}^{(1)}=\int\nolimits_{0}^{+\infty }{\eta }^{(1)}({\tau }^{{\prime} }){\Psi }_{{{{{{{{\rm{rec}}}}}}}}}({\tau }^{{\prime} })d{\tau }^{{\prime} }.$$
(43)