Robustness of higher-order networks with synergistic protection

From chemical reactions to human communications, higher-order interactions are ubiquitous in real-world networks. Entities within higher-order interactions often exhibit collective behaviors that could create synergistic effects on robustness of the underlying system. Here we propose an analytical model to investigate the robustness of higher-order networks, in which potential higher-order synergistic protection is incorporated. In this model, higher-order networks are described with simplicial complexes, and robustness is studied under the proposed dynamics of extended bond percolation. We provide theoretical analysis for robustness quantities including the relative size of the giant component and percolation threshold. We discover that the percolation threshold could drop to zero, which is an indicator of notably strong robustness, with synergistic protective effects and dense higher-order simplices. We also find that higher-order interactions have strong impacts on the association between robustness and clustering. Specifically, a larger clustering coefficient could invariably indicate stronger robustness once the strength of protective effects exceeds a certain value. Our theoretical solutions are verified by simulation results in simplicial complexes with Poisson, exponential and power-law distributions.


Introduction
In recent years, many real-world systems have been successfully modeled as complex networks, which provide convenient abstract representations for systems in terms of components and their interactions [1][2][3][4].One property of special interest is the fault tolerance of a system, which is usually characterized by the robustness of its underlying network structure [5].When a fraction of nodes or edges are attacked and removed, the network could be torn apart into several isolated fragments, resulting in disastrous damage to the functionality of the entire system.There have been extensive studies on network robustness since the beginning of this century, the majority of which formalize the problem under the unified framework of percolation theory [6][7][8][9][10][11][12].The basic paradigms in network percolation include site percolation [9,10] and bond percolation [11,12].The site percolation randomly removes a fraction of nodes along with all the edges attached to them, while the bond percolation directly removes a given fraction of edges in the network.The changes in network connectivity are observed during the percolation process, and the functionality of the network system is characterized by the existence of a percolating component, named as giant component.Network robustness can then be analyzed by studying the giant component size or the percolation threshold at which the percolating component ceases to exist [6,13,14].
Previous works usually define networks as a collection of nodes and edges, where interactions are naturally limited to dyadic ones.However, in many empirical systems such as neuroscience [15,16], biology [17], ecology [18,19] and social systems [20] to name a few, interactions may take place not only between pairs of nodes, but at the level of groups [21].All these systems are better characterized as higher-order networks such as simplicial complexes and hypergraphs [22], which intrinsically take into account the presence of higher-order interactions.Various studies have investigated dynamic processes that may occur in higher-order systems [23], such as spreading [24], diffusion [25] and synchronization [26].Percolation, in particular, has also been extended to higher-order networks, leading to many interesting results.Hu et al [27] proposed a theoretical framework for group percolation.They discovered that after assembling randomly selected nodes into groups, the robustness of the interdependent network is improved while the phase transition remains discontinuous.Sun and Bianconi [28] generalized percolation to multiplex hypergraphs and observed a wide range of critical behaviors and multiple percolation transitions.Ghoshal et al [29] further studied percolation on hypergraphs with hyperedges of fixed cardinality.K-core percolation was also extended and investigated on hypergraphs [30,31].Liu et al [32] presented a threshold model on hypergraphs and revealed the dual effect of hyperedges on network robustness.For simplicial complexes, Sun and Bianconi investigated percolation on several special types of geometric structures, such as hyperbolic simplicial complexes [33] and branching simplicial and cell complexes [34].Zhao et al [35] presented an analytical model and discovered that simplicial complexes with dense high-order simplices become highly fragile if nodes in simplicial components survive or fail altogether.
When concerning the relationship between robustness and higher-order structures, many previous studies emphasized the cascading effects of interdependent survival arising from higher-order interactions.For instance, one common assumption is that the survival of one node strongly relies on its neighbors' surviving states within one higher-order group, which could account for many cascading processes in real-world systems [27,32,35].There is, however, another kind of collective phenomenon that reveals an essential aspect distinct from the aforementioned interdependent survival: entities unite into groups to reinforce internal interactions and enhance resistance to external risks or attacks.These phenomena, aptly referred to as collaborative risk mitigation [36], can also be commonly observed in various real-life network systems [17,[36][37][38][39][40][41][42].For instance, in the financial market, each company has trading and sales connections with other companies that can be regarded as pairwise interactions within the network [43].In addition, it is natural to see that several companies share a common chairman or a major shareholder, as illustrated in figure 1(a).Companies within such higher-order relationships usually cooperate more closely, and the trading chains between them are thus less prone to be damaged due to market fluctuations or financial shocks.Take online social networks as another example.Individual users often communicate with their friends in private chat windows, forming an online social network with dyadic relationships [44].It is also common for some individuals to form a group chat or an online community together to facilitate communication owing to common hobbies or shared organizations, illustrated in figure 1(b).Contact or information exchange between members is then less possible to be cut down within such higher-order groups.Inspired by these collaborative risk-mitigating behaviors, we present the novel concept of synergistic protection, in order to better explore the potential impacts of higher-order interactions on network robustness.To be specific, we assume that higher-order interactions, once they exist, may exert potential synergistic protective effects on their covered elements such as edges and nodes.To the best of our knowledge, incorporating this kind of synergistic effects into network robustness analysis is still an understudied aspect.Bridging this gap by designing a more generalized model is thus of theoretical and practical significance for further detailed and in-depth study of network dynamics and robustness in the presence of higher-order interactions.
In this work, we develop a mathematical framework to study the robustness of higher-order networks with special attention to the synergistic protective effects brought by higher-order interactions.We utilize simplicial complexes, involving both dyadic and triadic interactions, to describe higher-order networks, and propose a triplex configuration method for network generation.Aiming at formalizing synergistic protection, we extend the classic bond percolation process by introducing a carefully designed 'attack and protect' mechanism as well as a generalized synergistic function.Extended bond percolation is then performed on simplicial complexes, and both analytical and numerical solutions to the relative size of the giant component as well as percolation threshold are given.Through extensive experiments, we discover that with moderate synergistic protective effects, simplicial complexes become more robust and the percolation threshold decreases to zero when the density of 2-simplices exceeds a certain value.We also demonstrate that the existence of higher-order interactions could impact the correlation between robustness and clustering of the model: a larger clustering coefficient could invariably indicate stronger robustness after the strength of protective effects reaches a particular threshold.This paper is organized as follows.In section 2, we first introduce the concepts of simplicial complex as well as the triplex configuration method for network construction.The extended bond percolation model is then presented, along with several robustness measures adopted in this work.Section 3 demonstrates the detailed analytical solutions to the relative size of the giant component as well as percolation threshold in the proposed extended bond percolation model.In section 4, we ground the analytical framework with numerical simulations on simplicial complexes that follow homogeneous and heterogeneous joint degree distributions.Quantitative impacts of higher-order The presence of higher-order interactions is represented by yellow dotted lines, while pairwise edges are drawn in solid lines.Edges are reinforced with synergistic protection if they are covered by higher-order interactions, represented by the change in color from blue to yellow.synergistic protection are studied, with the main findings illustrated in detail.Section 5 concludes the whole work with the outlook for future research.

Simplicial complex and triplex configuration method
We first give a brief introduction to simplicial complexes, a natural tool for characterizing higher-order networks [21].Just like a graph is a collection of edges, a simplicial complex is a collection of simplices.A simplex is a complete subgraph of a certain order in a network, also called a clique in graph theory.Specifically, a p-simplex is composed of p + 1 completely connected nodes, where all subsets of the p + 1 nodes are also contained in the p-simplex.In Euclidean space, we can simply consider a 0-simplex as a point, a 1-simplex as an edge, and a 2-simplex as a filled triangle, illustrated in figures 1(a)-(c).2-simplices are the most simple yet fundamental higher-order units in the network.What calls for special note is that a 2-simplex is different from the combination of three connected 1-simplices, illustrated in figure 2(d), since the former involves an extra three-body interaction and ensures the presence of the latter, yet the converse is not true.Based on the definition of simplices, a simplicial complex, denoted as G n = (C 0 , C 1 , . .., C n ), is a network of simplices, where n is the highest order of its simplices and C 0 , C 1 ,. .., C n are the set of 0-simplices, 1-simplices, . .., n-simplices in the network, respectively.In this work, we focus on simplicial complexes G 2 = (C 0 , C 1 , C 2 ), which includes both pairwise and three-body interactions.In keeping with convention, we hereinafter refer to 0-simplices as nodes and 1-simplices as edges.To distinguish between 2-simplices and three connected 1-simplices, which are both in a triangular shape, we refer to the latter structure as '3-1-simplex' , and the former as 2-simplex without change.
We next propose a triplex configuration method, a direct extended version of the classic configuration method [45,46] and the simplicial configuration method [47], for the construction of simplicial complexes.The major improvement of our method lies in the ability to explicitly distinguish between 2-simplices and 3-1-simplices, providing convenience for the study of higher-order synergistic effects in simplicial complexes.
In the triplex configuration model, we specify the number of single edges (1-simplices that are not part of 3-1-simplices or 2-simplices), the number of 3-1-simplices and the number of 2-simplices attached to each node respectively.Specifically, we define s i to be the number of single edges node i belongs to, h i and t i to be the number of 3-1-simplices and 2-simplices that node i participates, leading to the conventional degree k i = s i + 2h i + 2t i .The complete joint degree sequence (s i , h i , t i ) specifies the connecting information for node i, and is generated by a given joint degree distribution p s,h,t .For simplicity, we suppose the value distributions of s i , h i and t i are uncorrelated.Simplicial complexes could then be constructed with triplex configuration method following algorithm 1. Figure 2(e) illustrates an example of a constructed simplicial complex, with a joint degree sequence (2, 1, 1) assigned to the centered red node of the network.(named as a 3-1-simplex).Note that a 2-simplex is distinguished from a 3-1-simplex since the former includes an additional three-body interaction.(e) Example of a simplicial complex constructed by triplex configuration method.The centered red node is assigned a joint degree sequence (2, 1, 1), meaning that the node has two single edges, one 3-1-simplex and one 2-simplex attached to it.

Extended bond percolation
We propose the extended bond percolation model, incorporating the synergistic protection into classic bond percolation process.We split the edge removal step into a two-stage 'attack and protect' process, with the aim of formalizing the assumption that edges within 2-simplices could obtain additional synergistic protection during attacks: (1) we first select a fraction of edges to attack according to a given attack ratio γ ∈ [0, 1]. ( 2) we next consider whether each attacked edge is protected: if the attacked edge is within a 2-simplex, then it obtains extra synergistic protection that may prevent the edge from failing; attacked edges that are not protected fail directly, which is the same case in conventional bond percolation.After the 'attack and protect' process, we find the percolating cluster with classical graph searching algorithms such as breadth-first searching and compute the relative size of giant component S. We note that the process after 'attack and protect' is the same as the classic bond percolation if we simply project the network into its 'simplex-reduced' graph structure.
We next quantify the synergistic protection during the 'protect' stage.We define ϕ ∈ [0, 1] as synergistic protective strength, representing the probability that a protected edge within 2-simplices may survive after being attacked.ϕ → 1 implies the protected edge becomes completely immune to attacks while ϕ → 0 indicates that three-body interactions make no difference during bond percolation.The whole process of extended bond percolation is exemplified in figure 3, and can be simulated following algorithm 2.

Robustness evaluation
Here we introduce the different measures of network robustness utilized in this work.

(1) Relative Size of the Giant Component
The percolation theory mainly studies the emergence of the giant component under different occupation probabilities or attack ratios [6].In this work, the giant component refers to the largest connected component in which nodes are connected via any type of simplices.The relative size of the giant component is adopted to evaluate network structure integrity, defined as where N is the number of nodes in the complete network and M is the number of nodes in the giant component after percolation. (

2) Percolation Threshold
With the decreasing attack ratio γ, there exists a critical threshold γ c above which a non-zero S can be found.The percolation threshold is simply p c = 1 − γ c .A large percolation threshold p c indicates that just a small amount of failed nodes or edges can dismantle the network system, meaning that the robustness is poor, and vice versa.
For numerical experiments, it is often impractical to calculate p c by directly locating the first value that reduces S to zero owing to limited network scales and stochasticity.An alternative approach is to record the The four protected edges survive with a probability that equals the synergistic protective strength ϕ = 0.5, leading to two survival edges (edges 13,15).(e) The giant component, colored in red, is then found in the remaining simplicial complex with size M = 10, resulting in the relative size of the giant component S = 0.625.(f) For comparison, in classic bond percolation where no synergistic protection is involved (i.e.ϕ = 0), the remaining network would have two giant components with equal size M = 5, leading to S = 0.3125 which is much smaller than 0.625.fluctuations of the order parameter S with a proper normalization, named as susceptibility χ [48,49].The location where susceptibility reaches its peak provides an accurate estimate of the percolation threshold value p c .In this work, we follow Colomer-de-Simón and Boguñá [49].and define susceptibility χ as (3) R-index Despite the fact that percolation threshold p c well describes the robustness of a network from the view of phase transition, it fails to characterize robustness in situations where the system suffers partial damage and undergoes an incomplete collapse [50].R-index [50,51], in contrast, is a measure that considers the relative size of the giant component in the face of all possible attacks.In the standard version, R-index is designed for site percolation and is defined as where n = γN is the number of removed nodes during attacks and S(n) is the corresponding fraction of nodes in the giant component.
In this work, we modify this measure for bond percolation and define R-index as where L is the total number of pairwise edges in the network and l = γL is the number of edges being attacked during bond percolation.

Analytical framework
We provide analytical solutions to the relative size of giant component S as well as percolation threshold p c in the extended bond percolation model, with the help of generating functions.We first define the generating function of the joint degree distribution p s,h,t as A further quantity involved in this framework is the excess degree distribution [46].Conventionally, the excess degree is a property of a node reached by traversing a randomly chosen edge in a network, and is defined to be the number of edges connected to that node excluding the one it traverses.Here we define three types of excess degrees, that are the excess degree of a node reached by traversing a single edge, a 3-1-simplex and a 2-simplex respectively.Distributions of these excess degrees are calculated where ⟨s⟩, ⟨h⟩ and ⟨t⟩ are the average number of single edges, 3-1-simplices and 2-simplices a node belongs to within the network.The generating functions for the three excess degree distributions are written as We now present the procedure of calculating the giant component size with the generating functions.In simplicial complexes constructed by triplex configuration method, a node belonging to the giant component at least connects to one single edge, 3-1-simplex or 2-simplex that leads to the percolating cluster.We specify three vital probabilities to be calculated, denoted as u, v and w, for each type of connection.For any randomly chosen node i, we first define u to be the mean probability that one random neighbor of node i is not a member of the giant component when this neighbor is reached by traversing a random single edge of node i.We then let v to be the mean probability that one neighboring node is not in the giant component by traversing a 3-1-simplex of node i.Since one 3-1-simplex involves two neighbors of node i, we can naturally use v 2 to denote the mean probability that a 3-1-simplex of node i fails to lead node i to the giant component.Similarly, w and w 2 are defined as the corresponding mean probabilities for 2-simplices.
We next demonstrate how the values of u, v and w can be calculated in detail.If a given single edge l of node i fails to lead to the giant component, one of the following two cases must be true: either the edge l is attacked and failed, or the node at the other end of edge l, denoted as node j, is not a member of the percolating component via any of its other single edges, 3-1-simplices or 2-simplices.The first case happens with probjability of the attack ratio γ.The probability of the second case happening depends on the excess degrees of node i, i.e. the degrees of node j.Suppose node j owns s + 1 single edges (including edge l), h 3-1-simplices and t 2-simplices, then the probability of the second case equals (1 − γ)u s v 2h w 2t .Recalling that s, h and t follow the excess degree distribution m s,h,t , we can simply average over this distribution to reach the self-consistent equation for u The equation for v is more complicated since each of the three edges in the 3-1-simplex can be functional or damaged, bringing about a few possible cases after edge removal.Summing up all cases, we derive the following equation for v 2 The derivation for the equation of w is slightly different since it involves the additional synergistic function f.The two-step 'attack and protect' process for edges covered by 2-simplices can be equivalently considered as reducing the attack ratio to (1 − ϕ)γ, and we can then reuse the form of the previously derived equation.The equation for w 2 is written as Equations ( 12)-( 14) together constitute the core self-consistent equation system.The equations are often in complex nonlinear forms in real scenarios, but can be solved by numerical methods.
The probability of a node itself not being in the giant component is now able to be solved.A given node i with a joint degree sequence (s i , h i , t i ) does not connect to the giant component with probability u s i v 2h i w 2t i .Averaging over the joint degree distribution, we obtain the probability s,h,t P sht u s v 2h w 2t = G(u, v 2 , w 2 ), which is the mean probability that a randomly chosen node not belong to the giant component.Correspondingly, the expected relative size of the giant component in the entire simplicial complex is The critical value of attack ratio γ c that completely dismantles the giant component, along with the percolation threshold p c = 1 − γ c , is further investigated.Given the fact that there exists a trivial solution u = v = w = 1 that always satisfies equations ( 12)-( 14), the existence of a non-trivial solution, i.e. a solution of (u, v, w) with values unequal to (1, 1, 1), would indicate the presence of the giant component.The problem is then converted into finding the critical point where u, v and w first have values unequal to one.To this end, a straightforward way is to take the derivatives of the equations and impose equivalence when u, v and w take the trivial value We note that in practice, equation ( 16) could be complicated and in general has multiple solutions which can be puzzling for selecting the right γ c .Hence, performing variable elimination on the self-consistent equations, as in the case of triplex-Poisson distribution (see section 4.1), is preferred in order to simplify the computation of γ c and p c .

Results
We apply the extended bond percolation model to simplicial complexes with a variety of joint degree distributions, including both homogeneous and heterogeneous cases.Quantitative impacts of higher-order synergistic protection are studied through extensive numerical and analytical experiments.

Simplicial complex with homogeneous distribution
Simplicial complexes that follow homogeneous degree distributions are first investigated.Here we focus on Poisson distribution, a well-studied case in homogeneous random network models.In particular, simplicial complexes with triplex-Poisson degree distributions are analyzed, in which the numbers of single edges, 3-1-simplices and 2-simplices are all following an independent Poisson distribution where µ, ν and ξ are the mean values of Poisson random variables s, h and t, respectively.The average total degree can then be calculated as ⟨k⟩ = µ + 2ν + 2ξ.In this case, the generating functions of the joint degree distribution and excess degree distributions are easy to calculate and are coincidentally in the same form Substitute equation ( 18) into equations ( 12)-( 14), we reach the self-consistent equation system for triplex-Poisson distribution The calculation of S follows naturally after equation ( 19) is solved One property of interest in this case is that the values of S, γ c and p c can be solved in a simplified approach.Utilizing the equivalence between generating functions in equation ( 18), performing variable elimination becomes feasible.After eliminating variables u, v and w, equations ( 19) and (20) together are reduced to a single equation For simplicity, we denote the right-hand side of this equation as 1 − F(S; γ).This equation can be solved graphically.To this end, we construct an auxiliary function The solution of equation ( 21) corresponds to the crossing point of w(S) and S-axis.Since w(S) is continuous with w(0) = 0 and w(1) > 0, we can draw a qualitative curve of w(S) to observe whether a non-trivial solution exists.The critical values γ c and p c can then be easily found as it corresponds to the tangency of w(S) and S-axis with S = 0, written as dw (S) Compared with equation ( 16), solving equation (23) numerically is easier to implement.Figure 4(a) illustrates an example of the graphical solutions.We can see that when p ≈ p c , the curve is tangent to the S-axis, and a non-trivial solution exists when p > p c .We fix µ = 2.22, ν + ξ = 0.88 and change ρ, ϕ, respectively.The dotted lines are analytical solutions and symbols are simulation results.For each (ϕ, ρ) configuration, we run simulations on networks with 10 000 nodes and all estimates are averaged under 100 realizations.We can observe that high ϕ and ρ push the S curve upward and lead to smaller pc.For the case of (ϕ, ρ) = (0.5, 0.9), we observe that pc reduces to zero.(c) and (d) record the fluctuations for each group of 100 simulations as susceptibility χ.The peaks of χ match well with the analytical values of pc.
Quantitative impacts of synergistic protection on network robustness are further investigated by varying two model attributes: (1) the strength of synergistic protection ϕ, and (2) the density of higher-order interactions, defined as ρ = ξ ν+ξ , representing the proportion of 2-simplices among all triangle structures.With the increase of synergistic protection strength, the edges within 2-simplices obtain a bigger chance to survive in the face of attacks, and would decrease the overall speed of network collapse during percolation.A larger density of higher-order interactions implies that more edges in the network could be protected, leading to stronger resilience against attacks as well.Figures 5(a) and (b) display the numerical and analytical results of the giant component sizes in simplicial complexes versus attack ratio γ, from which one can find the analytical results are in good agreement with numerical simulation.Figures 5(c) and (d) show that the peaks of susceptibility χ also match the analytical values of percolation threshold p c , verifying the correctness of our framework.Meanwhile, it can be observed that as the density ρ or synergistic protection strength ϕ increases, the whole curve of S is shifted to the right and the percolation threshold is reduced, indicating the enhancement of the network robustness.One can also find that for the case of (ϕ, ρ) = (0.9, 0.5) , the percolation threshold has dropped to zero, which is an indicator of notably strong robustness.In figure 4(b) and (c), we provide direct results for percolation thresholds under different configurations.Again, results show that once the synergistic protection strength ϕ exceeds 0.45, the percolation threshold could drop to zero when density ρ is large enough.
For triplex-Poisson distributed simplicial complexes, we also investigate the relationship between synergistic protection and clustering coefficient, a structural property that is closely related to robustness [51][52][53][54].The reason why we did not directly select clustering coefficient as a robustness metric is that the answer to the correlation between clustering and robustness remains controversial in the literature.Intuitively, a high clustering coefficient could imply high robustness because the number of alternative paths increases along with the number of triangles [52].However, the opposite conclusion, i.e. increasing clustering coefficient results in decreasing robustness, is also reached in several works [51,53].Studying the correlation between these two properties is thus of great interest especially after synergistic protection is involved.By simply ignoring the existence of three-body interactions, we project the simplicial complexes onto their underlying graph structures, and calculate the global clustering coefficient c, which by convention is defined as where N △ and N 3 represent the number of triangles and the number of connected triples in the network, respectively [45].For the projected graph in our model, we can calculate N △ and N 3 as follows where t ′ = h + t is the total triangle degree, k = s + 2h + 2t is the conventional total degree, G t ′ (x, y) = G(x, y, y) is the generating function for s and t ′ , and G k (x) = G(x, x 2 , x 2 ) is the generating function for k.Substituting equation (18) into equations ( 24)-( 26), we obtain the global clustering coefficient c for triplex-Poisson joint distribution Experiments are conducted with fixed average total degree ⟨k⟩ and higher-order interaction density ρ. Figure 6(a) shows the corresponding curves of S under varied synergistic protection strength ρ as well as clustering coefficient c.For each fixed ϕ, one interesting observation is that the curves with different c always intersect at the same point, and the location of the intersection point moves right as ϕ increases.On the left side of the intersection point, higher clustering brings a larger giant component size and lower percolation threshold, indicating the reinforcement of robustness.On the right side, however, the giant component size becomes smaller as clustering coefficient grows, leading to an opposite conclusion.To avoid such conflict and reach a straightforward insight, we further refer to R-index as the robustness measure.Results for R-index in different configurations are illustrated in figure 6(b).One can observe that as the synergistic protection strength ϕ enhances, the network with c = 0.2 changes from the most vulnerable system to the most robust one in terms of R-index.This phenomenon reveals that the potential higher-order protective effects could strongly influence the relationship between clustering and robustness.There exists a certain threshold for ϕ, below which clustering and robustness are negatively correlated.As ϕ exceeds this threshold, however, a large clustering coefficient becomes invariably indicative of stronger robustness.

Simplicial complex with heterogeneous distribution
We next study simplicial complexes with heterogeneous distributions, in which the degrees of nodes display a high variance.The exponential distribution is first considered.In this case, simplicial complexes follow exponential edge degree distribution p s , given by where κ is a constant.For continuity of illustration, we suppose that the 3-1-simplices and 2-simplices of simplicial complexes remain Poisson distributed.The generating functions for this joint degree distribution are as follows Substituting equation ( 29) into equations ( 12)-( 16), we can calculate the percolation quantities to measure robustness.
Power-law distribution is studied as the second case.Specifically, we consider power-law edge degree distribution with an exponential cutoff, which is defined as where τ, κ are constants and Li n (x) = ∞ k=1 x k k n is the nth polylogarithm of x.The reasons for including exponential cutoff are twofold: firstly this cutoff is shown in many real-world networks [46,55]; secondly it makes the normalization of the distribution possible for all power exponent τ , not just τ ⩾ 2 [46].The corresponding generating functions can then be conveniently written down in closed form, which are Similarly, with these generating functions and equations ( 12)-( 16), the model can be solved theoretically.Analytical solutions of the giant component sizes for the aforementioned two kinds of simplicial complexes are illustrated in figures 7(a) and (b), compared with simulation results.We note that the phenomena of percolation thresholds diminishing to zero could be observed in both distributions when protective strength ϕ is large enough, similar to the previous case of the triplex-Poisson distribution.We also calculate the R-index for both distributions in figure 7(c) for comparison.Results show that simplicial complexes with power-law edge degree distribution are more robust than the exponential counterpart of the same average degree.This finding is consistent with the well-established conclusion that scale-free networks exhibit remarkable robustness against random attacks.It can also be observed that the advantage of power-law simplicial complexes in terms of R-index is maintained with the increase of synergistic protective strength ϕ in figure 7(c).

Discussion
In this work, we have proposed a theoretical framework to study the robustness of simplicial complexes.Inspired by the collective risk-mitigating behaviors, our model takes into account that individuals often form into groups to resist external attacks and we establish a hypothesis that these higher-order groups may exert synergistic protection on their group members and their interconnections.We describe higher-order networks with simplicial complexes constructed by triplex configuration method, and propose the extended bond percolation model utilizing an 'attack and protect' process to incorporate synergistic protective effects.We present exact analytical solutions for various robustness measures including the relative size of the giant component and percolation threshold.Simulation results closely match our analytical solutions in simplicial complexes with Poisson, exponential and power-law distributions.
We show that after synergistic protection is incorporated, two interesting observations can be made.First, as the density of 2-simplices increases, we find that the percolation threshold could gradually diminish to zero in both homogeneous and heterogeneous cases.This suggests that higher-order networks could obtain remarkably strong robustness with the aid of synergistic protection, providing inspiration for designing robust network structures when higher-order interactions are included.Second, we discover that the strength of synergistic protection could alter the correlation between clustering and robustness.To be specific, a larger clustering coefficient could indicate stronger robustness in terms of R-index when protective strength exceeds a certain threshold, below which an opposite conclusion would be reached.Since clustering is a common property in real-world networks and is intimately related to higher-order structures as well, our finding provides an insightful angle to further inspect the robustness of clustered networks in terms of higher-order synergistic protection.
Despite the valuable insights gained from this study, there are two major limitations.On the one hand, our framework only considers simplicial complexes of dimension two for simplified modeling, but in reality, higher-order interactions could manifest in more complex forms beyond three-body relationships.Extending our model to the general simplicial complex structure is thus of great significance.One feasible and natural approach to achieve this is to extend the proposed triplex configuration method to the multiplex version and adopt the n-dimensional generating function method for theoretical analysis, which shares similarities with the works by Karrer and Newman [56] and Zhao et al [35].On the other hand, the synergistic protective strength ϕ in the model is simplified as a constant, but the actual form could be complicated.In this framework, it is applicable to set ϕ as a function of the model parameters, such as ϕ (ξ), ϕ(⟨k⟩) or ϕ(ρ).Exact solutions for robustness quantities could be obtained by simply replacing the symbol ϕ with the specific functional form ϕ(x) in equation (14).Since the exact mechanism of higher-order synergistic protection may vary in different real-world situations, further research could focus on exploring the concrete forms of synergistic protection in the context of specified scenarios, such as social networks [20] and ecological systems [17].Meanwhile, defining a series of synergistic functions for higher-order interactions of different dimensions, such as defining ϕ 2 (x), ϕ 3 (x), . .., ϕ n (x) for 2-simplices, 3-simplices, . .., n-simplices, respectively, is also worthy of investigating.We expect that by incorporating multi-dimensional synergistic protection rules, higher-order networks could exhibit remarkable robustness against a number of specific attack strategies, which would provide insights for the design of robust higher-order network systems.
In summary, we have introduced the novel concept of synergistic protection in the study of higher-order network robustness, and have provided an analytical framework to precisely solve several robustness quantities in simplicial complexes.Our work contributes to understanding the dynamics of higher-order networks from a long-neglected perspective and sheds light on the mechanisms underlying the relationship between collective risk-mitigating behaviors and network robustness.We hope that our findings will inspire further investigations in the field of higher-order network analysis, especially from the perspective of synergistic protection.

Figure 1 .
Figure 1.Schematic illustration of higher-order synergistic protection in (a) a financial network and (b) an online social network.The presence of higher-order interactions is represented by yellow dotted lines, while pairwise edges are drawn in solid lines.Edges are reinforced with synergistic protection if they are covered by higher-order interactions, represented by the change in color from blue to yellow.

Figure 2 .
Figure 2. Geometrical representations of a (a) 0-simplex, (b) 1-simplex, (c) 2-simplex and (d) a set of three connected 1-simplices(named as a 3-1-simplex).Note that a 2-simplex is distinguished from a 3-1-simplex since the former includes an additional three-body interaction.(e) Example of a simplicial complex constructed by triplex configuration method.The centered red node is assigned a joint degree sequence (2, 1, 1), meaning that the node has two single edges, one 3-1-simplex and one 2-simplex attached to it.

Figure 3 .
Figure 3. Schematic diagram of the extended bond percolation dynamics.(a) A simplicial complex is initially constructed using triplex configuration method, with N = 16 nodes, eight single edges, two 3-1-simplices and two 2-simplices, leading to L = 20 pairwise edges in total.(b) In the attack stage, given the attack ratio γ = 0.4, eight randomly selected edges are attacked.(c) In the protection stage, four protected edges that are within 2-simplices (edges 2, 3, 13, 15) obtain synergistic protection, while the other four unprotected edges (edges 5, 8, 16, 17) fail directly.(d)The four protected edges survive with a probability that equals the synergistic protective strength ϕ = 0.5, leading to two survival edges (edges13,15).(e) The giant component, colored in red, is then found in the remaining simplicial complex with size M = 10, resulting in the relative size of the giant component S = 0.625.(f) For comparison, in classic bond percolation where no synergistic protection is involved (i.e.ϕ = 0), the remaining network would have two giant components with equal size M = 5, leading to S = 0.3125 which is much smaller than 0.625.

Figure 4 .
Figure 4. (a) Example of graphical solutions for triplex-Possion distributed simplicial complexes.The curves are drawn for different values of occupation rate p = 1 − γ under configuration µ = 2.22, ν = ξ = 0.44, ϕ = 0.2.For p = 0.1 < pc, there is only one trivial intersection S = 0.As we increase the value of p, the curve w(s) gradually moves downward.When p = 0.1844 ≈ pc, the curve is tangent to the S-axis, corresponding to the critical scenario.When p = 0.3 > pc there exists a nontrivial solution indicated by the second intersction point on the S−axis.(b) Analytical and numerical results of percolation threshold pc in simplicial complexes with triplex-Poisson distribution under configuration µ = 2.22, ν + ξ = 0.88 and different ϕ, ρ.We can observe that when synergistic protective strength ϕ exceeds 0.45, the percolation threshold could reach zero with large density ρ.(c) Detailed illustration of percolation thresholds under ρ = 0.3, 0.6 and 0.9.It can be observed that for larger ρ, pc drops faster along with increasing synergistic protective strength ϕ.

Figure 5 .
Figure 5. (a) and (b) show the relative sizes of the giant component in simplicial complexes with triplex-Poisson degree distribution.We fix µ = 2.22, ν + ξ = 0.88 and change ρ, ϕ, respectively.The dotted lines are analytical solutions and symbols are simulation results.For each (ϕ, ρ) configuration, we run simulations on networks with 10 000 nodes and all estimates are averaged under 100 realizations.We can observe that high ϕ and ρ push the S curve upward and lead to smaller pc.For the case of (ϕ, ρ) = (0.5, 0.9), we observe that pc reduces to zero.(c) and (d) record the fluctuations for each group of 100 simulations as susceptibility χ.The peaks of χ match well with the analytical values of pc.

Figure 6 .
Figure 6.(a) The relative size of the giant component in triplex-Poisson distributed simplicial complexes with different clustering coefficients.We set the average total degree ⟨k⟩ = 4 and clustering coefficient c = 0.0, 0.1 and 0.2 for each (ϕ, ρ) configuration.We observe that for each fixed clustering coefficient c, the three curves always intersect at the same point, which gradually moves up and left as ϕ increases.(b) Values of R-index for each clustering coefficient c in (a), with changing synergistic strength ϕ.Results show that from ϕ = 0.0 to ϕ = 0.5, the R-index value of networks with c = 0.2 rises from the lowest to the highest, representing an inversion of the correlation between robustness and clustering.

Figure 7 .
Figure 7. (a) and (b) show the relative size of the giant component S in simplicial complexes with exponential and power-law edge degree distribution, respectively, as a function of synergistic protective strength ϕ and attack ratio γ.The degree distributions of 3-1-simplices and 2-simplices remain double-Poisson for continuity of illustration, with ν = ξ = 0.5 (ρ = 0.5) in both cases.The average edge degree in both cases is fixed to ⟨s⟩ ≈ 3 by setting the constants κ = 3.2 in (a) and τ = 1.2, κ = 10 in (b), respectively.(c) shows the R-index values of both kinds of simplicial complexes under different ϕ.It can be observed that simplicial complexes with power-law distribution exhibit stronger robustness compared with the exponential counterpart in terms of R-index, and this advantage is maintained with the increase of synergistic protective strength ϕ.

Algorithm 1 .
Construction of Simplicial Complexes using Triplex Configuration Method.Input: joint degree distribution p s,h,t , number of nodes N Output: simplicial complex SC Create arrays S, H, T that each store N integers and array SC to store N empty node structures; Generate N joint degree sequences that follow P sht , and store values in arrays S, H, T in order; while elements in array S are not all 0 do Input: simplicial complex SC, attack ratio γ, synergistic protective strength ϕ Output: relative size of giant component S Enumerate pairwise edges in SC and generate edge list L; Randomly select ceil(γ • L.length) edges from L to generate attacked edge list L ′ ; foreach edge {i, j } in L ′ do if ∃k, {i, j } in SC[k].two_simplices then if rand() < ϕ then continue; Remove 2-simplex {i, j, k} from SC; else if ∃k, {i, j } in SC[k].three_one_simplices then Remove 3-1-simplex {i, j, k} from SC; else if i in SC[ j].one_simplices then Remove 1-simplex {i, j} from SC;