Percolation transitions in partially edge-coupled interdependent networks with different group size distributions

1. Introduction

Complex networks provide a powerful paradigm for characterizing complex systems because they can be used to simulate and study a wide variety of real-world network systems, including the Internet, social networks and food chains.^1,2,3,4 In most real-world scenarios, systems are not isolated but coupled; a typical scenario is that of cyber-physical systems (CPSs).^{5,6,7,8,9,10,11} Although the interactions between networks enable such systems to operate more efficiently, they also introduce many risks. The interdependency of interdependent networks can significantly increase their vulnerability to cascading failures; i.e. if either external or internal perturbations disrupt one system, the resulting failures could spread not only within the disrupted system itself but also within the other systems with which it interacts, leading to abrupt failures and catastrophic collapses, like the 2003 Italian blackouts and the 2020 global COVID-19 breakout.^12,13

Since the introduction of the framework proposed by Buldyrev et al. in 2010 for analyzing the robustness of interdependent networks based on percolation theory,¹² numerous studies have been conducted within this framework. To date, the models of interdependent systems have primarily concentrated on the pairwise dependencies between individual nodes. Considering the structural and functional patterns, the models above vary from structural and functional patterns, including but are not limited to coupling strength,¹⁴ coupling schemes,¹⁵ and coupling preferences.^16,17 Additionally, these models have explored various dependency relationships, including multiple support^18,19 and multiple dependencies,²⁰ as well as various network properties, such as directed,^21,22 weighted,^23,24 and capacity load.^11,25

Most studies above, however, have focused on node-coupled interdependent networks (NINs), where an interdependency relationship exists between nodes from different layers. However, interactions usually happen on the edges connecting nodes rather than at the nodes themselves in many real-world networks. For example, congestion in the backbone of a transportation network can also cause failures to propagate to other routes between cities. To enhance the performance of designated backbone networks, such as an intermodal bus and rail transport network, bridging services may be needed to temporarily replace some parts of a metro network.^26,27 Therefore, considering the relationships between edges from different layers has become an emerging and promising topic in studying interdependent networks.^{28,29,30,31,32,33}

For studying the edge-coupled interdependent networks (EINs), in which the functioning of an edge depends on the functioning of edges in other layers, was proposed by Gao et al.²⁸ and a mathematical framework based on percolation theory was established. They found that EINs have a lower threshold than NINs and that a broader degree distribution increases the robustness of EINs, a finding that contrasts with the results observed for NINs. Later, Gao et al.²⁹ proposed a model for partially edge-coupled interdependent networks (partially EINs) to facilitate the study of arbitrarily coupled EINs, aiming to represent the case in which only some edges in a network depend on edges in another network better. This work revealed how changes in the dependency strength q can trigger alterations in the types of phase transitions occurring in networks. Then, Zhou et al.³¹ analyzed the mechanism that makes EIN generally more robust than NIN based on a quenched network framework and showed that this property is rooted in the fact that in a network, the excess degree of an edge is on average larger than the degree of a node. Subsequently, Xie et al.,³⁰ and Gao et al.³² developed a weighted and a directed EIN model, respectively, and also established comprehensive theoretical frameworks to enrich the research on percolation transitions of EINs. In real-world scenarios, relationships between systems may not only contain pairwise interactions but also exist in high-order interactions.^34,35,36 specific nodes may sometimes band together into cliques or groups to improve the robustness of a system against risk.^{37,38,39,40,41,42,43,44,45,46,47} For example, in service function chaining (SFC), multiple virtual network functions (VNFs) are deployed on the same physical node and tend to fail or survive together.^37,38,39 In networks with groups, there exist two kinds of connections. One can be regarded as a hyperedge containing several regular nodes that survive or fail together; another is the link that connects two nodes. So pairwise and high-order interactions exist together in the network with groups.

The theory of group percolation, which studies the above phenomenon, has thus been widely developed to improve further the robustness of interdependent networks with groups (NGINs) and prevent their abrupt collapse. Wang et al.⁴⁰ proposed a general model of “group percolation” in interdependent networks and found that the formation of groups can significantly improve the network robustness. From the perspective of coupling schemes, Wang et al.⁴² investigated the effects of different coupling schemes and found that the network robustness improves as the assortativity of interdependent networks with cliques increases. Similarly, Su et al.⁴³ formulated a cascading failure model for interdependent networks with cliques and multiple dependency relationships and found that the system’s robustness improves as the number of multiple dependency relationships increases. Notably, the type of phase transition remains unchanged for homogeneous degree distributions and multiple dependency relations but can be altered by heterogeneous degree distributions. From the perspective of dependency strength, Zang et al.⁴⁴ explored the impact of weak interdependence on interdependent networks with cliques and found that reducing the dependency strength can cause the type of the phase transitions in both layers to change from discontinuous phase transitions to continuous phase transitions. From the attack strategies, Peng et al.⁴⁸ proposed a series of edge attack and removal strategies and investigated the impact under different strategies in studying the robustness of NGINs.

In addition, realistic interdependent networks are not as fragile as most realistic networks because the network can be heterogeneous, and this heterogeneity is reflected in the probability of failure of the connected edges or the size of the supernode.^45,46,47 Lu et al.⁴⁵ and Zang et al.⁴⁶ analyzed the effects of different clique size distributions on network robustness under random and targeted attacks, respectively, and found that increasing the heterogeneity of the clique size distribution can improve network robustness. More recently, Li and colleagues⁴⁷ examined the combined effects of reinforcement of critical nodes and the dependency group distribution and concluded that suitable modifications to the density of reinforced critical nodes and the heterogeneity of the dependency group could significantly enhance the robustness of networks and cause them to exhibit multiple phase transitions.

However, we do not know whether the strength of dependence and heterogeneity of supernodes size affects the robustness of the EIN network, as it did before in the NINs and NGINs. To address such scenarios, based on the original EIN model, we develop a model of group percolation in partially edge-coupled interdependent networks with groups (partially EGINs) with different group size distributions to fill this gap in this paper. Our study reveals that dependency strength and heterogeneity in group size distribution can greatly affect the robustness of interdependent networks. These findings provide valuable theoretical insights for improving the robustness of such networks in real-life applications.

This work is structured as follows. In Sec. 2, we introduce the two-layer partially EGINs model. In Sec. 3, we present theoretical derivations to investigate the sizes of the giant components (GCs) and the phase transition properties of the proposed model. In Sec. 4, we report the results for the general case of EGINs with groups of equal size (HomEGINs) as well as for the case of EGINs with groups of two different group size distributions (HetEGINs). Finally, we offer conclusions and further discussion in Sec. 5.

2. Model Description

In this section, we develop a model of group percolation in a partially EGIN that consists of two layers, namely, layer A, with node set $N_A$ and edge set $E_A$, whose degree distribution is $P_A(k_A)$, and layer B, with node set $N_B$ and edge set $E_B$, whose degree distribution is $P_B(k_B)$. As described in Ref. 40, we randomly divide the nodes in layers A and B into small, none-overlapping groups (i.e. groups with no shared nodes) with sizes of $m_1$ with equal size and $m_2$ with a given group size distribution $g(m_2)$, respectively, which can be regarded as supernodes. After this transformation, layers A and B are changed into new networks $Ã$ and $\tilde{B}$, respectively, each with a number of supernodes equal to $N_{Ã}=N_{\tilde{B}}=N$. In the constructed networks $Ã$ and $\tilde{B}$, each supernode corresponds to a group, and the size of each supernode is equal to the number of regular nodes in the corresponding group. If edges are present between the original regular nodes, corresponding edges are also added between the supernodes to which these nodes belong after the network is divided; note that in this process, edges are not added repeatedly between the same two supernodes. In this model, a fraction $q_1$ of the edges in network $Ã$ are randomly dependent on edges in network $\tilde{B}$ (through dependency links), and a fraction $q_2$ of the edges in network $\tilde{B}$ is similarly defined, as described in Ref. 49. It is assumed that each network edge has at most one dependent edge. If an edge in subnet A and an edge in subnet B exhibit a one-to-one correspondence, this is referred to as a one-to-one no-feedback dependency. For explicit illustration, we first consider transformed networks $Ã$ and $\tilde{B}$, where $Ã=A$ with a supernode size of $m_1=1$, and assume that the nodes of layer B cooperate and form groups with $m_2=2$. The illustration of EGIN and its bipartite network representation are shown in Fig. 1.

3. Theoretical Analysis

In this section, we theoretically derive the giant components (GCs) and giant grouped components (GGCs) of the proposed partially EGINs. Based on group percolation theory established in Ref. 40, the generation functions of the transformed networks $Ã$ and $\tilde{B}$ can be expressed as $G_{Ã0}(x)={\sum }_{k_{Ã}}{P}_{Ã}(k_{Ã})x^{k_{Ã}}$ and $G_{\tilde{B}0}(x)={\sum }_{k_{\tilde{B}}}{P}_{\tilde{B}}(k_{\tilde{B}})x^{k_{\tilde{B}}}$, with degree distributions of $P_{Ã}(k_{Ã})$ and $P_{\tilde{B}}(k_{\tilde{B}})$, respectively. Here, $P_{Ã}(k_{Ã})$ or $P_{\tilde{B}}(k_{\tilde{B}})$ denotes the probability that an arbitrarily selected supernode in network $Ã$ or $\tilde{B}$ will have a degree of $k_{Ã}$ or $k_{\tilde{B}}$, respectively. $〈k_{Ã}〉$ and $〈k_{\tilde{B}}〉$ represent the average degrees.

In addition, we define the probability of randomly selecting a pair of supernodes in networks $Ã$ and $\tilde{B}$ with sizes $m_1$ and $m_2$, respectively, as $Q(m_1,m_2)$. In this model, we assume that the size of each supernode in network $Ã$ is set to $m_1=1$, while the sizes of the supernodes in network $\tilde{B}$ follow a distribution $g(m_2)$. Then, the expression for $Q(m_1,m_2)$ can be simplified to $Q(1,m_2)$, where $Q(1,m_2)=g(m_2)$. The generating function corresponding to $Q(m_1,m_2)$ can be expressed as

where M is the largest size among the supernodes in layer B. When the sizes of both layers approach infinity, i.e. $N_A\to \infty$ and $N_B\to \infty$, and the ordinary nodes in layers A and B are randomly divided into none-overlapping groups, the probability of internal links existing between nodes belonging to the same node group approaches zero in the thermodynamic limit and becomes negligible. We can derive the relationship between $D(x,y)$ and $G(x,y)$ as follows :where $G_A(x)={\sum }_{k_A}{P}_A(k_A)x^{k_A}$ and $G_B(y)={\sum }_{k_B}{P}_B(k_B)y^{k_B}$ are the generating functions corresponding to the degree distributions of the original networks $Ã$ and $\tilde{B}$, respectively.

For a partially EGIN with different group size distributions in the two layers, we define $f_{Ã}$ and $f_{\tilde{B}}$ as the probabilities of reaching the GC in one direction along a randomly selected edge in networks $Ã$ and $\tilde{B}$, respectively. For the constructed network $Ã$, which has a degree distribution of $P_{Ã}(k_{Ã})$, the probability that a supernode with a degree of $k_{Ã}$ does not lead to the GC is ${(1-f_{Ã})}^{k_{Ã}}$. By averaging the probability ${(1-f_{Ã})}^{k_{Ã}}$ overall possible degrees $k_{Ã}$, the probability that a supernode belongs to the GC can be calculated as ${\sum }_{k_{Ã}}^{\infty }P(k_{Ã})[1-{(1-f_{Ã})}^{k_{Ã}}]$, which is also the size of the GC, denote as ${\mu }_{\infty }$. In addition, if a supernode along a selected edge leads to the GC, i.e. at least one of the remaining $k_{Ã}-1$ edges of that supernode leads to the GC, then the selected edge is also considered to lead to the GC. Accordingly, the probability that a randomly selected supernode belongs to the GC is ${\sum }_{k_{Ã}}\frac{P_{Ã}(k_{Ã})k_{Ã}}{〈k_{Ã}〉}[1-{(1-f_{Ã})}^{k_{Ã}-1}]$.

Different from traditional node-coupled interdependent networks (NGINs),⁴⁰ our study focuses on the partially EGINs, which involve coupled relationships between the edges connecting supernodes. Therefore, it makes more sense to analyze this network based on the bond percolation theory.⁸ Gao et al.²⁸ using the self-consistent probabilities method to derive the percolation equations of the fully EINs, where the coupled relationships exist between the edges that connect regular nodes. Following this approach, we establish an analysis framework of the EGINs model proposed in this paper. To calculate the value of ${\mu }_{\infty }$, we need to obtain the expressions for $f_{Ã}$ and $f_{\tilde{B}}$. For a randomly selected edge “$l_a$” in network $Ã$, suppose that not only does it belong to the GC, but its dependent partner edge “$l_b$” in network $\tilde{B}$ also exists and belongs to the GC. These two edges are thus part of the mutually connected giant component (MCGC). For network $\tilde{B}$, it is easy to find the probability of the opposite case, that is, the probability that neither end of the edge is in the GC. From the above analysis, we know that the probability that a randomly selected supernode “b” in network $\tilde{B}$ is not in the GC is ${\sum }_{k_{\tilde{B}}}\frac{P_{\tilde{B}}(k_{\tilde{B}})k_{\tilde{B}}}{〈k_{\tilde{B}}〉}{(1-f_{\tilde{B}})}^{k_{\tilde{B}}-1}$. Furthermore, because each edge is connected to two supernodes, one at each end, we know that the probability that a randomly selected edge “$l_b$” in network $\tilde{B}$ does not belong to the GC is ${[{\sum }_{k_{\tilde{B}}}\frac{P_{\tilde{B}}(k_{\tilde{B}})k_{\tilde{B}}}{〈k_{\tilde{B}}〉}{(1-f_{\tilde{B}})}^{k_{\tilde{B}}-1}]}^2$; then, the probability of the opposite case, namely, the probability that edge “$l_b$” belongs to the GC, is $1-{[{\sum }_{k_{\tilde{B}}}\frac{P_{\tilde{B}}(k_{\tilde{B}})k_{\tilde{B}}}{〈k_{\tilde{B}}〉}{(1-f_{\tilde{B}})}^{k_{\tilde{B}}-1}]}^2$.

Since this model is not fully but partially EGINs, we assume that there is a probability of $q_1$ that an edge in network $Ã$ is dependent on an edge in network $\tilde{B}$ and a probability of $1-q_1$ that an edge in network $Ã$ is an autonomous edge. First, we randomly remove a proportion $1-p$ of the edges in network $\tilde{B}$ to trigger a cascading failure process. Therefore, for a randomly chosen edge “$l_a$” in network $Ã$, the probability that this edge is in the GC is equal to the probability that one of the following conditions holds: (i) it is an autonomous edge and leads to the GC or (ii) it is a nonautonomous edge, but its dependent partner in network $\tilde{B}$ was not removed after the initial attack and also leads to the GC. Thus, the self-consistent probability equation for subnet $Ã$ in the partially EGINs can be transformed as follows using the generating function $G(x,y)$ described in Eq. (3)

where $G_x(x,y)=\frac{\partial G(x,y)}{\partial x}$ is the partial derivative of $G(x,y)$ with respect to x and $G_{xy}(x,y)=\frac{\partial G(x,y)}{\partial x\partial y}$ is the simultaneous partial derivative of $G(x,y)$ with respect to both x and y.

Similarly, we can obtain the self-consistent equation for $f_{\tilde{B}}$ as follows :

Because there is no interaction between the supernodes in the two layers, the expressions for ${\mu }_{\infty }^{Ã}$ and ${\mu }_{\infty }^{\tilde{B}}$ are the same as in a single network, as illustrated in what follows :

The probability that a randomly chosen regular node in layer A is connected to the GGC is denoted by ${\mu }_A$. Since layer A is isomorphic to network $Ã$, it holds that ${\mu }_A={\mu }_{\infty }^{Ã}$. In layer B, where the supernodes have sizes of $m_2$ with a distribution of $g(m_2)$, it also follows that ${\mu }_B={\mu }_{\infty }^{\tilde{B}}$.

4. Percolation Transition Properties and Critical Tipping Points

In this section, we focus on analyzing the effect of dependency strength and heterogeneity on the robustness of the EGIN model. First, we present an analysis and discussion of the percolation thresholds and critical tipping points in HomEGINs, where supernodes of equal size exist within subnetworks $Ã$ and $\tilde{B}$. We then explore the HetEGINs model in terms of heterogeneity by investigating the effect of two different group size distributions on the robustness of EGINs. The simulation results are in excellent agreement with our theoretical predictions, highlighting the validity of our model.

Numerical simulations are performed on ErdsRnyi (ER) and scale-free (SF) networks, commonly used to represent real-world network systems. In an ER-ER partially EGINs, the distributions of the network degrees satisfy $P_A(k_A)=e^{-〈k_A〉}{〈k_A〉}^{k_A}∕k_A!$ and $P_B(k_B)=e^{-〈k_B〉}{〈k_B〉}^{k_B}∕k_B!$, and the generating functions are $G_{A0}(x)=G_{A1}(x)=e^{〈k_A〉(x-1)}$ and $G_{B0}(x)=G_{B1}(x)=e^{〈k_B〉(y-1)}$. In an SF-SF partially EGINs, the degree distribution satisfies $P(k)=c{k}^{-\lambda }$ for both networks, where $\lambda$ is the degree exponent and $k_{min}\le k\le K$. In this latter case, $P(k)={(k_{min}∕k)}^{\lambda -1}-{(k_{min}∕k+1)}^{\lambda -1}$ is a reasonable approximation of the degree distribution. The generating functions for SF networks are as follows^50,51 :

4.1. HomEGINs

In this part, we first analyze how to determine the critical dependency strength $q_2$ and the percolation threshold $p_c$ for a given $q_1$, and show the characterization of different types of phase transitions and finally give the phase diagrams with respect to p and $q_2$ in the following sections. Without loss of generality, in the case of partially HomEGINs, the number of nodes in layer A is set to be $N_A=10^6$ and the number of nodes in layer B is set to be $N_B=2\times 10^6$. The sizes of supernodes in network $Ã$ are fixed to be $m_1=1$ and the sizes of supernodes in network $\tilde{B}$ are fixed to be $m_2=2$ in the simulations.

4.1.1. Determination of critical tipping points

Theoretically, Eqs. (4) and (5) can be reformulated as $f_{Ã}=F_1(p,f_{\tilde{B}},f_{Ã})$ and $f_{\tilde{B}}=F_2(p,f_{Ã},f_{\tilde{B}})$, respectively. At the point of a discontinuous phase transition for the system, the following condition applies :

For a continuous phase transition, the size of the GC experiences a continuous increase at the transition point $p_c^{\mathrm{\text{II}}}$. Because the system exhibits a critical tipping point from a discontinuous phase transition to a continuous phase transition, Eqs. (4) and (5) yield two distinct types of solutions, as follows:

(i)

When $p\to p_{cÃ}^{\mathrm{\text{II}}}$, the GC of network $Ã$ is zero, and the GC of network $\tilde{B}$ is nonzero ($f_{Ã}=0$ and $f_{\tilde{B}}>0$). In this case, a Taylor expansion of Eq. (4) at $f_{Ã}\to 0$ is performed as follows : $$f_{Ã}={F_1}^{\prime }(p_c^{\mathrm{\text{II}}},f_{\tilde{B}},f_{Ã})f_{Ã}+\frac{1}{2}{F_1}^{\prime \prime }(p_{cÃ}^{\mathrm{\text{II}}},f_{\tilde{B}},f_{Ã}){f_{Ã}}^2+o({f_{Ã}}^3),$$(9) which leads to $${F_1}^{\prime }(p_{cÃ}^{\mathrm{\text{II}}},f_{\tilde{B}},0)=\frac{[1-q_1\cdot (1-p_{cÃ}^{\mathrm{\text{II}}})]G_x^{\prime }(1,1)}{G_x(1,1)}-\frac{q_1\cdot p_{cÃ}^{\mathrm{\text{II}}}\cdot G_{xy}^{\prime }(1,1-f_{\tilde{B}})}{G_{xy}(1,1)},$$(10) where $f_{\tilde{B}}$ can be regarded as a constant $f_{c\tilde{B}}$. In addition, substituting $f_{Ã}=0$ into Eq. (5) yields $$f_{c\tilde{B}}=p_{cÃ}^{II}\cdot \left[1-\frac{G_y(1,1-f_{\tilde{B}})}{G_y(1,1)}\right]-q_2\cdot p_{cÃ}^{II}\cdot \left[1-\frac{G_{xy}(1,1-f_{\tilde{B}})}{G_{xy}(1,1)}\right].$$(11) By combining Eqs. (10) and (11), the solutions for $p_{cÃ}^{II}$ and $f_{c\tilde{B}}$ can be derived. At the phase transition point, the network’s phase transition threshold satisfies $p_{cÃ}^{\mathrm{\text{I}}}=p_{cÃ}^{\mathrm{\text{II}}}$. Additionally, by concurrently solving Eq. (8), the junction point $q_{2cÃ}^{\mathrm{\text{I-II}}}$ can be determined. This is the point at which the phase transition type of network $Ã$ shifts from discontinuous to continuous and is also the point at which network $\tilde{B}$ transitions from a discontinuous phase transition to a hybrid phase transition, i.e. $q_{2cÃ}^{\mathrm{\text{I-II}}}=q_{2c\tilde{B}}^{\mathrm{\text{I-H}}}$. (ii) When $p\to p_{c\tilde{B}}^{\mathrm{\text{II}}}$, the GCs of networks $Ã$ and $\tilde{B}$ are both zero ($f_{Ã}=0$ and $f_{\tilde{B}}=0$). In this case, a Taylor expansion of Eq. (5) at $f_{\tilde{B}}\to 0$ is performed as follows : $$f_{\tilde{B}}={F_2}^{\prime }(p_c,f_{Ã},f_{\tilde{B}})f_{\tilde{B}}+\frac{1}{2}{F_2}^{\prime \prime }(p_c,f_{Ã},f_{\tilde{B}}){f_{\tilde{B}}}^2+o({f_{\tilde{B}}}^3).$$(12) When $f_{Ã}\to 0$, $f_{\tilde{B}}\to 0$, which leads to $${F_2}^{\prime }(p_{c\tilde{B}}^{\mathrm{\text{II}}},0,0)=p_{c\tilde{B}}^{\mathrm{\text{II}}}\cdot \left\{\frac{G_y^{\prime }(1,1)}{G_y(1,1)}-q_2\cdot \frac{G_{xy}^{\prime }(1,1)}{G_{xy}(1,1)}\right\}.$$(13) Then, we have $$p_{c\tilde{B}}^{II}=\frac{1}{\frac{G_y^{\prime }(1,1)}{G_y(1,1)}-\frac{q_2\cdot G_{xy}^{\prime }(1,1)}{G_{xy}(1,1)}}.$$(14)

The threshold $p_{c\tilde{B}}^{\mathrm{\text{II}}}$ can be found by solving Eq. (14). At this phase transition point, because the network phase transition threshold satisfies $p_{c\tilde{B}}^{\mathrm{\text{I}}}=p_{c\tilde{B}}^{\mathrm{\text{II}}}$, by simultaneously considering Eqs. (10) and (11), we can determine the solution $q_{2c\tilde{B}}^{H-\mathrm{\text{II}}}$ where network $\tilde{B}$ transitions from a hybrid phase transition to a continuous phase transition.

Graphical illustrations provide further clarification of this process. Taking ER-ER partially EGINs as an example, the curves represented by Eqs. (4) and (5) are plotted, and the values of $f_{Ã}$ and $f_{\tilde{B}}$ at the intersection are labeled, as shown in Fig. 2. Specifically, for $q_2=0.8$, which is greater than $q_{2c}$ in both networks $Ã$ and $\tilde{B}$, an examination of Fig. 2(a) shows that as the proportion p of preserved edges increases, $f_{Ã}$ and $f_{\tilde{B}}$ abruptly shift from zero to nonzero values, representing discontinuous phase transitions in both networks $Ã$ and $\tilde{B}$. For $q_2=0.4$ in Fig. 2(b), where $q_2$ decreases to $q_{cÃ}^{\mathrm{\text{I-II}}}<q_2<q_{2c\tilde{B}}^{H-\mathrm{\text{I}}}$, the value of $f_{Ã}$ starts at zero and undergoes a sudden change when the two curves become tangent, indicating a discontinuous phase transition in network $Ã$, while the value of $f_{\tilde{B}}$ varies continuously from zero to nonzero and then abruptly shifts to a new value when the two curves become tangent, indicating a hybrid phase transition in network $\tilde{B}$. For $q_2=0.1$, where $q_2<q_{cÃ}^{\mathrm{\text{I-II}}}$ in Fig. 2(c), as the proportion p of preserved edges increases, the two curves continuously intersect without tangency, resulting in values of $f_{Ã}$ and $f_{\tilde{B}}$ that do not show abrupt changes. The networks $Ã$ and $\tilde{B}$ undergo continuous phase transitions.

4.1.2. Characterization of different types of phase transitions

Figure 3 gives the behaviors of ${\mu }_A({\mu }_B)$ as a function of p for ER–ER and SF–SF networks. Essential properties between NINs, EINs, NGINs and EGINs are revealed by comparing their performance, as can be seen in Fig. 3. In both ER–ER and SF–SF networks, the EGINs showed better robustness than the remaining three models, NINs, EINs and NGINs, which can be reflected by $p_c(\mathrm{\text{EGIN}})<p_c(\mathrm{\text{EIN}})<p_c(\mathrm{\text{NGIN}})<p_c(\mathrm{\text{NIN}})$. Our model is based on EGIN and then further investigates how the dependency strength affects the robustness of the EGINs.

Figures 4 and 5 depict the phase transition behaviors in both networks $Ã$ and $\tilde{B}$ with respect to the proportion of preserved edges p when $q_2$ takes values of 0.1, 0.3, 0.5, 0.7 and 0.9. As depicted in Figs. 4(a) and 4(b), for a fixed dependency strength $q_1=1$, the fraction of the giant grouped component of layer A and the fraction of the giant grouped component of layer B increase with the increase of the proportion of preserved edges p. For another aspect, when p is fixed, ${\mu }_A$ and ${\mu }_B$ increase as the decrease of $q_2$, and the system exhibits different types of phase transitions. As can be seen in these figures, when $q_2$ is equal to 0.9 or 0.7, both ${\mu }_A$ and ${\mu }_B$ exhibit a discontinuous phase transition. When $q_2$ decreases to 0.5 or 0.3, the percolation diagrams of ${\mu }_A$ and ${\mu }_B$ exhibit distinct differences, with ${\mu }_A$ still showing a discontinuous phase transition, while the behavior of ${\mu }_B$ indicates a hybrid phase transition. Specifically, in network $\tilde{B}$, a continuous phase transition occurs first, followed by a discontinuous phase transition. From Figs. 4(a) and 4(b) and Figs. 5(a) and 5(b), we can draw the preliminary conclusion that as $q_2$ decreases for both ER–ER and SF–SF EGINs, network $Ã$ gradually changes from a discontinuous phase transition to a continuous phase transition, while network $\tilde{B}$ changes first from a discontinuous phase transition to a hybrid phase transition and then from a hybrid phase transition to a continuous phase transition.

4.1.3. Phase diagrams with respect to p and $q_2$

Next, in Fig. 6, we intuitively present the variations in the sizes of ${\mu }_A$ and ${\mu }_B$ in the $(p,q_2)$ plane for both networks $Ã$ and $\tilde{B}$ in the cases of ER–ER and SF–SF HomEGINs. It is observed that ${\mu }_A$ and ${\mu }_B$ increase with increasing p and decrease with increasing $q_2$. A red curve divides each panel into two regions, namely, the functional and nonfunctional regions, where a GC is present in the functional region but absent in the nonfunctional region. Specifically, as shown in Figs. 6(a) and 6(c), for network $Ã$, when $q_2>q_{2cÃ}^{\mathrm{\text{I-II}}}$, network $Ã$ undergoes a discontinuous phase transition, and the corresponding threshold $p_{cÃ}^{\mathrm{\text{I}}}$ can be determined by solving Eq. (8). On the other hand, when $q_2\le q_{2cÃ}^{\mathrm{\text{I-II}}}$, network $Ã$ undergoes a continuous phase transition, and the corresponding threshold $p_{cÃ}^{\mathrm{\text{II}}}$ can be determined by combining Eqs. (10) and (11). For network $\tilde{B}$, as shown in Figs. 6(b) and 6(d), when $q_2>q_{2c\tilde{B}}^{\mathrm{\text{I}}-H}$, a discontinuous phase transition occurs in $\tilde{B}$, with the corresponding threshold $p_{c\tilde{B}}^{\mathrm{\text{I}}}$ being the same as that for network $Ã$. When $q_{2c\tilde{B}}^{H-\mathrm{\text{II}}}<q_2<q_{2c\tilde{B}}^{H-\mathrm{\text{I}}}$, $\tilde{B}$ undergoes a hybrid phase transition. In this transition, as p decreases, ${\mu }_B$ first jumps discontinuously to a nonzero value, which occurs when the curves of $f_{Ã}$ and $f_{\tilde{B}}$ become tangent to each other, and then continuously decreases from nonzero to zero at $p_{c\tilde{B}}^{H-\mathrm{\text{II}}}$, the point where Eq. (12) intersects with the point $(0,0)$. Furthermore, when $q_2<q_{2c\tilde{B}}^{H-\mathrm{\text{II}}}$, network $\tilde{B}$ undergoes a continuous phase transition, and the corresponding threshold $p_{c\tilde{B}}^{\mathrm{\text{II}}}$ can be determined from Eq. (14).

In Fig. 7, we further show the phase diagrams of the systems for different average degrees $〈k〉$ in ER–ER partially HomEGINs (solid lines) and NGINs (dashed lines) and for different degree exponents $\lambda$ in SF–SF partially HomEGINs (solid lines) and NGINs (dashed lines). In general, the robustness of the proposed EGINs model is better than the traditional NGINs, which is reflected in the percolation threshold $p_c^{\mathrm{\text{I}}}$ and $p_c^{\mathrm{\text{II}}}$ of EGINs are lower than that of NGINs. The meanings of $q_{2cÃ}^{\mathrm{\text{I-II}}}$ ($q_{2c\tilde{B}}^{H-\mathrm{\text{II}}}$), and $q_{2c\tilde{B}}^{\mathrm{\text{I}}-H}$ remain consistent with those depicted in Fig. 6. It is noteworthy that the threshold $q_{2c\tilde{B}}^{\mathrm{\text{I}}-H\prime }$ within NGINs exceeds that of EGIN (hereby referred to as $q_{2c\tilde{B}}^{\mathrm{\text{I}}-H}$). The boundary $q_{2c\tilde{B}}^{\mathrm{\text{I}}-H\prime }$ divided the phase transition characteristics of layer $\tilde{B}$ into two classes: discontinuous and hybrid. It suggests that layer $\tilde{B}$ in NGIN commences its hybrid phase transition at an earlier stage compared to its counterpart in EGIN. These diagrams explore the relationship between $q_2$ and the phase transition threshold $p_c$. It is observed that for both ER–ER partially HomEGINs and NGINs, a higher $〈k〉$ results in greater robustness. Conversely, for both SF–SF partially HomEGINs and NGINs, an increase in $\lambda$ has a detrimental impact on network robustness. Additionally, both Figs. 7(a) and 7(b) illustrate that a decrease in the interdependency strength, $q_2$, leads to a decrease in the phase transition threshold and can cause the types of percolation transitions observed in these systems to change.

4.2. HetEGINs

In this section, we present an extensive series of numerical simulations to investigate the robustness of interdependent networks with heterogeneous group sizes (HetEGINs) to assess the validity and effectiveness of the theoretical framework established in Sec. 3. Without loss of generality, we fix the number of nodes in network B to $N_B=2\times 10^6$. Accordingly, when the group distribution $g(m_2)$ changes, the value of $N_A=N$ also changes. In our investigation, we consider two different types of $g(m_2)$ distributions, namely, a Gaussian heterogeneous distribution and a power-law heterogeneous distribution.

4.2.1. Gaussian distribution

In this part, we explore the effect of Gaussian-distributed group sizes of the supernodes on the robustness and phase transition behaviour of the system. The Gaussian distribution is expressed as follows :

where ${\bar{m}}_2$ is the mean and $\sigma$ is the variance. When $\sigma =0$, the heterogeneous network transforms into a homogeneous network, with a $g(m_2)$ distribution that is a delta function, indicating that all supernodes are of the same size. Figure 8 presents the simulation results for ER–ER HetEGINs in which the size distributions of the supernodes in network $\tilde{B}$ follow the aforementioned Gaussian distribution. As depicted in this figure, the phase transition point, denoted by $p_c$, monotonically decreases, while the network robustness increases as the variance parameter of the Gaussian distribution increases. A larger variance indicates a more heterogeneous network. Consequently, there is a positive relationship between network robustness and network heterogeneity. This finding helps to explain why real complex networks exhibit superior robustness against catastrophic failure compared to theoretical expectations. Here, our theory (solid lines) can predict the numerical simulations (empty symbols) well.

4.2.2. Power-law distribution

In this part, we further investigate another case of HetEGINs in Fig. 9. Specifically, we assume that the group sizes of $m_2$ in network $\tilde{B}$ follow a power-law distribution and explore its impact on the robustness and phase transition behavior of the system. The power-law distribution is expressed as follows :

where the exponent $\beta$ is the power-law parameter of group size distribution in network $\tilde{B}$. Specifically, decreasing $\beta$ corresponds to increasing heterogeneity of group sizes. The values of ${\mu }_A$ and ${\mu }_B$ increase as $\beta$ decreases, implying that ${\mu }_A$ and ${\mu }_B$ positively correlate with heterogeneity. Once again, our theoretical predictions agree well with the numerical results.

5. Conclusions

This paper investigates a cascading failure model for partially edge-coupled interdependent networks with heterogeneous group size distributions. Instead of interdependencies between supernodes, this model considers interdependencies between the edges connecting supernodes in each network layer. To analyze the phase transition behaviors of such systems, we develop a mathematical framework for the proposed HetEGINs using generating functions based on the self-consistent probability approach. We first analyze the tipping points and determine the corresponding critical dependency strength $q_{2c}$ in ER–ER and SF–SF HomEGINs, where all supernodes in the same layer consist of the same number of ordinary nodes. We also find that networks with a larger number of supernodes (network $\tilde{B}$ in this paper) exhibit multiple phase transition properties, including continuous, discontinuous, and hybrid percolation behaviors, in different ranges of the dependency strength value. Then, by assigning different group size distributions to network $\tilde{B}$, i.e. the Gaussian distribution and the power-law distribution, we found that systems with higher heterogeneity are more robust for a given average group size. In general, based on the constructed EGIN model, we not only discuss the influence of any interdependence strength $q_2$ on the change of phase transition type of $Ã$ and $\tilde{B}$ but also introduce the heterogeneous distribution of group size into the model, which further proves the positive effect of heterogeneous network shape on the robustness of the network. These findings provide an essential reference for designing more robust networks in real-life edge-coupled network scenarios. In addition, we can continue to consider multiple dependencies, the correlation between edge coupling dependencies, extend the model toward higher-order networks, and so on. Future work still needs to be extended to understand edge-coupled interdependent networks.

Acknowledgments

This work is supported by Program of Song Shan Laboratory (Included in the management of Major Science and Technology Program of Henan Province) (221100210700-2).

ORCID

Junjie Zhang  https://orcid.org/0000-0002-0014-1650

Caixia Liu  https://orcid.org/0009-0006-5977-6171

Shuxin Liu  https://orcid.org/0000-0002-3366-4466

HaiTao Li  https://orcid.org/0009-0004-1004-0142

Lan Wu  https://orcid.org/0009-0003-6455-6673