A possible quadruple point in networks of directed networks under targeted attacks

Many real systems are known to interact with one another, forming networks of networks (NONs). Plenty of attention has been poured into the research on the robustness in NONs in the past decade. Previous studies focus on undirected networks, or directed networks under random attacks. While many real networks are directed and how networks of directed networks (NODNs) respond to targeted attacks remains unknown. We thus develop a general analytical tool for analyzing the robustness of NODNs under two kinds of targeted attacks: degree-based attacks and in-degree (out-degree)-based attacks. The analytical tool can perfectly predict the sizes of the final giant strongly connected components and the phase transitions on the NODNs in response to targeted attacks. By applying the tool to synthesis networks, we find that a quadruple point intersected by four different phase regions could appear in the random regular NODNs. To the best of our knowledge, it is the first time that a quadruple point is found in the studies of complex networks. In addition, we find triple points intersected by three phases in networks of directed scale-free networks, and critical points that connect two phases in networks of directed Erdös–Rényi networks. The discovery of these tipping points could help understand network robustness and enable better design of networked systems.


Introduction
As a powerful tool to study the characteristics and behaviors of complex systems, network science has attracted the attention of many researchers over the past two decades [1][2][3][4][5]. One of the key characteristics for complex networks is robustness, which is important for many fields [6][7][8][9][10]. Perturbations acting on networks may seriously endanger the systems' robustness to retain their basic functions [11][12][13]. Networks exhibit different states in response to different kinds of perturbations [7,[14][15][16]. Generally, networks are more vulnerable to targeted attacks compared to random failures [12].
It has long been recognized that the interaction along systems has become increasingly intensive [8,9,[17][18][19]. Due to the interaction between networks, systems can exhibit characteristics different from those in single-layer networks [8]. For example, Italy experienced a blackout caused by cascading failures between the power grid and communication network. Buldyrev et al [8] thus developed a framework of two fully interdependent networks to study the robustness of system subject to cascading failures. The fully interdependent networks exhibit a discontinuous first-order phase transition, which is quite different from a second-order phase transition in a single network. Considering that the networks are not necessarily fully interdependent, Parshani et al [20] proposed a more realistic model of two partially interdependent networks and found a change from second-order phase transition to first-order phase transition when the coupling strength between networks increases from zero to one. Later on, a series of more general frameworks are proposed to analyze the robustness of networks of networks (NONs) which are composed of more than two networks [21,22].
Studies show that the coupling between networks reduces the robustness of systems [23,24]. However, most of the above studies is based on random failures. Considering that initial failure is mostly not random, Huang et al [25] studied the robustness of two fully interdependent networks in response to targeted attacks and found that it is difficult to maintain the robustness by protecting nodes with high degrees. The change from first-order to second-order phase transition also appears in partially interdependent networks when the coupling strength decreases [26]. Then, the analysis framework of robustness is extended to the NONs [27].
All the studies discussed above focus on undirected networks, while many real networks are directed [28][29][30][31][32]. After considering the directionality, interdependent directed Erdös-Rényi (ER) networks show a hybrid phase transition which is not found in interdependent undirected ER networks [9]. Recently, a theoretical framework was developed to analyze the robustness of two interdependent directed networks in response to two types of targeted attacks, i.e. degree-based attacks and in-degree/out-degree-based attacks. The robustness of interdependent directed scale-free (SF) networks first increases then decreases as the degree heterogeneity increases, which is different from the case of random attacks [33]. However, the robustness of networks of directed networks (NODNs) under targeted attacks has not been addressed.
In this paper, we develop a general framework for analyzing the robustness of NODNs under targeted attacks. Here, we focus on a kind of NODN that has a random regular structure with loops. In the random regular NODNs, each network depends on l neighboring networks (figures 1(a) and (b)). For the first time, we find that a quadruple point appears in the phase diagram of the networks of directed ER networks with l < 5. The quadruple point is intersected by four phase transition regions, i.e. the second-order, the first-order, the hybrid, and the collapse phase transition. The observation of the quadruple point indicates that the collapse phase transition may appear even in two interdependent directed ER networks, which is not observed in previous studies. When l < 5, the second-order and the hybrid phase transition intersect at a duple point. Note that when l = 4, there are two duple points and the system does not exhibit the first-order phase transition. As for the networks of directed SF networks, the collapse phase transition only appears when l ⩾ 3. Besides, a triple point intersected by the hybrid, the first-order, and the collapse phase transitions is found when l ⩾ 3. Thus, networks with high heterogeneity first show a second-order, then a hybrid, and finally a collapse phase transition as the coupling strength increases.

Modeling of the random regular NODNs
We build a model of NODN, which is composed of n directed networks G i , with i = 1, 2, . . ., n. Each network G i consists of N i nodes and can be characterized by the joint degree distributions P i (k in , k out ), where k in and k out are the in-degree and out-degree of a given node respectively. As shown in figures 1(a) and (b), in a random regular NODN, each network is connected to the same number l of neighboring networks. For each connected network pair G i and G j , they are connected by dependency links, with q ij ⩾ 0 fraction of nodes in network G i depending on the nodes in network G j (nodes in network G j point to nodes in network G i ), and q ji ⩾ 0 fraction of nodes in network G j depending on the nodes in network G i ( figure 1(c)). Once the source node of a dependency link stops functioning, the target node of the dependency link fails immediately. In addition, the nodes from two networks are coupled following the no-feedback condition [34], that is, a node from one network can depend on no more than one node from the other network, and if a node is the source and target of dependency links, then the corresponding target and source of these dependency links are the same node.
Once 1 − p i fraction of nodes are removed from each network G i , the network fragments, and we assume that only the nodes located in the giant strongly connected component (GSCC) [35] are functional. A strongly connected component (SCC) is a subnetwork, where each node has at least one directed path to reach every other node and can be reached by every other node by a directed path. The GSCC is the largest SCC in the network and contains a significant fraction of the entire network. In addition, the couplings between networks enable failures to propagate back and forth between networks, forming cascading failures. When the cascading failures stop, only the nodes in the final GSCC remain functional. The final GSCC size is related to the initial attack method. Next, we will show how to theoretically analyze the effect of degree-based and in-degree-based attacks on network robustness.

Targeted attacks on isolated directed networks
We begin by intentionally removing 1 − p fraction of nodes from each network of the NODN according to nodes' degrees (in-degrees, or out-degrees). In the degree-based attack, each node in network G i is removed with the following probability [25,33] where k i j is the degree of node j of network G i , and α i is a parameter that controls the probability of each node to be attacked. We use a general technique that maps the targeted-attack problem to a random-attack problem in isolated directed networks.
In a directed network G with degree distribution P(k in , k out ), the generating function of the degree , where x and y are arbitrary complex variables. We remove 1 − p fraction of nodes with the probabilities following equation (1) and the links between the removed nodes, while keeping the links between the removed and the remaining nodes. In this way, the generating function of the degree distribution in the remaining network G 1 is where In the second step, we remove the links that connect the removed nodes and remaining nodes, then get remaining network G 2 . Since network G is randomly connected, the probability for a link that connects between two remaining nodesp is equal to the ratio of the number of links connect to the remaining nodes to the total number of links in the original network:p Removing the links that connect to the removed nodes is equivalent to randomly removing 1 −p fraction of links of the remaining nodes. The generating function of the degree distribution after randomly removing 1 −p fraction of links from G 1 is equal to [36] The next step is to find a network G eq with the generating function being Φ Geq0 (x, y), such that after randomly removing 1 − p fraction of nodes from network G eq , the generating function of the remaining nodes is the same as Φ G2 (x, y). That is, Φ Geq0 (x, y) = Φ G2 (1 − p + px, 1 − p + py). According to equation (4), we could get ) .
In a similar way, we can get the generating function of the equivalent networks after in-degree-based (or out-degree-based) targeted attacks is where Φ G1 (x, y) is the generating function of the degree distribution in the remaining network G 1 with initially removing 1 − p fraction of nodes and the links between the removed nodes but remaining the links between the removed and the remaining nodes(see section S1 of the supplemental material for further details).

Cascading failure in NODNs
As we mentioned above, once a fraction of nodes are removed from each network, the failure could propagate between networks, forming cascading failures. Next, we develop an analytical tool to calculate the fraction of the final functional GSCC. Let N i denote the set of all the networks connected with network G i . Define r ij,t as the fraction of remaining nodes of network G i after the network G i suffers the damages from all networks connected to G i , except the network G j (G j ∈ N i ) at time step t. At time step t = 1, each network G i suffers damage due to the initial removal of 1 − p i fraction nodes but no damage from other networks. Thus, r ij,1 = p i for G j ∈ N i . At time step t ⩽ 1, G i suffer damage from all of its neighbor networks, and the damage from its neighbor network G j is q ji [1 − r ji,t−1 g j (ψ ′ j,t−1 )]. As a result, the fraction of the remaining nodes in .
After cascading failures, the system reaches a stationary state at time step τ and the fraction of the remaining nodes in network G i [10] is where and g j (x) is used to calculate the size of GSCC in isolated network G i (see section S2 of the supplemental material). Thus, the final size of GSCC in the stationary situation for each network G i G i [10] is In the random regular NODNs, each network depends on the same number l of neighboring networks. For simplicity and without loss of generality, we assume that q ij = q, p i = p, and g i (ψ ′ τ ) = g(ψ ′ τ ). Thus, the sizes of final GSCC for all the networks are the same, and follow. By simplify the equations (7) and (8), we can get And the final size of GSCC (see section S2 of the supplemental material) is where z in can be solved by As shown in figure 2, the final GSCC size is zero when the fraction of the remaining nodes p is small. When p is greater than the percolation thresholds p I c or p II c , the final GSCC size increases as p increases. The robustness of the system can either be defined by the values of the percolation threshold (the smaller threshold is, the higher robustness is), or defined as the integer size of the final GSCC during the entire attack process (the greater integer size is, the higher robustness is), shown as the area under the curve in figure 2.
For the network of directed ER networks(l = 3) with the same average degree ⟨k⟩, the final GSCC size under degree-based and in-degree-based targeted attacks (α = 1) are p As the coupling strength q increases from 0 to 1, the network of directed ER networks exhibits different phase transitions (see section S4 of the supplemental material). And we can find that the network of directed ER networks are less robust in response to targeted attacks compared to random attacks, but show similar robustness in response to degree-based and in-degree-based attacks. Similarly, the network of directed SF networks exhibits three types of phase transitions for degree-based attacks and in-degree-based attacks. Seen from figures 2(c) and (d), the network of SF networks under degree-based attacks is more vulnerable than that under in-degree-based attacks(see section S5 of the supplemental material).

Percolation thresholds and phase transitions in NODNs
As the fraction of remaining nodes after the initial attack p decreases, the NODNs show various percolation behaviors under different coupling strengths q: (1) the final size of the GSCC continuously decreases to zero at a percolation threshold p II c , showing a second-order phase transition, when q < q c2 ; (2) the final size of the GSCC discontinuously jumps to zero at a percolation threshold p I c , showing a first-order phase transition, when q c1 < q < q max ; (3) the final size of the GSCC first discontinuously jumps to a small value at p I c , and then continuously decreases to zero at p II c , showing a hybrid phase transition; (4) collapse when q > q max , where the system collapses even if only one node is removed. The parameters q c2 , q c1 , and q max are called critical coupling strengths which separate the second-order, hybrid, first-order, and collapse phases. The percolation thresholds p II c and p I c and the critical coupling strengths q c1 , q c2 and q max can be analytically computed using the function R (s) (z, q) and its derivations (see section S3 of the supplemental material).
We draw the curves of percolation thresholds p I c and p II c for network of the directed ER networks (l = 3) with ⟨k⟩ = 6, ⟨k⟩ = 8 and ⟨k⟩ = 12, respectively. In figures 3(a) and (b), we can see that p II c increases and disappears at the critical coupling strength q c1 , and p I c appears at the critical coupling strength q c2 and increases as the coupling strength q increases from 0 to 1. In addition, only the percolation threshold p II c exists when ⟨k⟩ = 6, which represents there are no hybrid and first-order phase transitions(see section S4 of the supplemental material).
The robustness of the networks of SF networks is related to their degree heterogeneity, which can be characterized by the exponent of the degree distribution λ. The lower λ is, the higher heterogeneity of the degree distribution is. The percolation thresholds for networks of directed SF networks under degree-based attacks are shown in figure 3(c). When the coupling strength is very low, the robustness of NODNs is close to that of isolated networks. In isolated SF networks with low λ, there are a few nodes that have very high degrees, leading to low robustness under degree-based attacks [11]. Thus, the robustness of the networks of SF networks with higher degree heterogeneity tends to be lower when the coupling strength is low, as shown in figure 3(c). When the coupling strength is high, the heterogeneity in degree distribution could help improve the system's robustness, which is in accord with the finding in interdependent networks [33]. In the case of in-degree-based attacks, the networks of SF networks with higher degree heterogeneity have lower percolation thresholds, which means higher robustness, as shown in figure 3(d).
We compute the critical coupling strengths q c1 , q c2 , and q max and draw the phase diagram for the network of directed ER networks and SF networks. Figures 4(a) and (b) show the phase diagram of the network of ER networks where l = 3. The phase diagram is divided into the region of the second-order phase transition(labeled 'Phase II'), the region of the hybrid phase transition(labeled 'Hybrid'), the region of the first-order phase transition(labeled 'Phase I') and the region of collapse(labeled 'Collapse') by critical coupling strength q c2 , q c1 and q max . In addition, a quadruple point and two duple points appear in the phase diagrams, which will be discussed in the next sections.
In addition, we draw the phase diagrams(figures 4(c) and (d)) of the network of directed SF networks (l = 3) attacked according to nodes' degrees and in-degrees respectively. We find that the critical coupling strengths q c1 and q c2 decrease as the exponent of degree distribution λ increases. While the value of the critical coupling strength q max is represented by two lines, which intersect at a triple point where λ t = 2.2343. When the exponent of the degree distribution λ < λ t , the critical coupling strength q max decreases as λ increases. When λ is greater than λ t , q max is a constant with q max = 1 l−1 = 0.5 (see section S5 of the supplementary material for details).

Emergence of a possible quadruple point intersected by four phases
The quadruple point represents a set of conditions under which four phases of a system can exist in equilibrium. In the phase diagram of the network of ER networks, we find a quadruple point that represents the intersection of four phase transition states, and this point should have all the characteristics of these four phase transition states (see section S4 of the supplemental material). Therefore, at the quadruple point, the system satisfies The values of ⟨k * 1 ⟩ and q * 1 for different l are shown in table 1. According to table 1, there are four-phase points in the network when l < 3. As shown in figure 5, four phase transition regions intersect at the quadruple point (9.3341, 0.4016). Choosing four points around the  quadruple point in different regions, we can see the simulation results agree well with the theoretical results, which confirms the existence of the quadruple point (shown in figure 5). The duple point (7.1421, 0.3092) are intersected by the hybrid region and second-order region. The quadruple point and duple point divide the system into three categories: A. ⟨k⟩ < ⟨k * 2 ⟩ = 7.1421; B. ⟨k * 2 ⟩ < ⟨k⟩ < ⟨k * 1 ⟩; C. ⟨k⟩ > ⟨k * 1 ⟩ = 9.3341. When the system belongs to class A, that is, the average degrees of the networks are small, the system only exhibits second-order. Thus, the percolation threshold p I c is absent from figure 3(a) when ⟨k⟩ = 6. When the system belongs to class B, the system first shows a second-order phase transition, then a hybrid phase transition, then a second-order phase transition again, and finally collapses as the coupling strength increase from 0 to 1. When the system belongs to class C, i.e. the average degree of the networks are large, the system can exhibit all four phase transitions(see section S4 of the supplemental material).

Triple points appearing in networks of directed SF networks
Referring to the triple point, we can think of the numerical point in thermodynamics that makes the gas phase, liquid phase, and solid phase coexist. Around this point, small fluctuations can lead to an abrupt change of network GSCC. It was found that three-phase points are found in the α − p phase diagram, where α represents the interdependence fraction of high-degree nodes, and p represents the size of functional nodes. And it indicated that the triple point is characteristic of nontrivial patterns of interdependency [37].
And we also find a triple point (λ * , q * ), which is intersected by the hybrid, the first-order, and the collapse regions, appear in the phase diagram of the network of directed SF networks. Power-law exponent λ * and coupling strength q * can be solved by The values of λ * and q * are shown in table 2. When l < 3, there is no collapse region and triple point, which is in accord with the result of two interdependent directed SF networks. When l ⩾ 3, the system robustness rapidly decreases. The collapse region and the triple point appear in the phase diagram. The larger l is, the larger λ * is while the smaller q *  Figures 4(c) and (d) show the phase diagrams of the network of directed SF networks with l = 3 under degree-based attacks and in-degree-based attacks, respectively. The phase diagram is also divided into four regions: collapse region (labeled by Collapse), first-order region (labeled by Phase I), hybrid region (labeled by Hybrid), and second-order region (labeled by Phase II). For high heterogeneous networks, there is no first-order phase transition. In addition, we find that the critical coupling strengths are the same for the two targeted attacks, which is caused by the fact that the phase is determined by the network topological properties and coupling strengths, regardless of the attacking method.

Conclusion
We develop a theoretical framework for analyzing the robustness of the NODNs under degree-based and in-degree (out-degree) based attacks. Our framework could perfectly predict the sizes of GSCC in steady state, percolation thresholds, and critical coupling strengths of the NODNs. The systems' robustness decreases dramatically as the number of coupled networks increases. A quadruple point, which is not discovered before, is found in the phase diagram of the networks of directed ER networks with l < 4. The quadruple point is intersected by first-order, second-order, hybrid, and collapse phase transition regions, indicating the system may collapse even without removing nodes (p = 1) when the networks are sparse. Besides, a duple point intersected by the regions of second-order and hybrid phase transitions appears in the phase diagram when l < 5. Especially, there are two duple points when l = 4, and the system cannot exhibit the first-order phase transition. When l ⩾ 5, both the quadruple point and the duple point disappear, and the system only exhibits second-order and collapse phase transitions. For the networks of directed SF networks, the system exhibits collapse phase transition only when l ⩾ 3. In addition, a triple point is discovered in the phase diagram. Thus, the phase transition changes from a second-order, through a hybrid to a collapse phase transition. Our findings could help us understand the complex system's robustness and indicate ways of pushing the system into a more robust region.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).