Interdependent networks can be defined as a sequence of graphs: G = {G _A , G _B …}. Usually, nodes in two or more monoplex networks are adjacent to each other via edges that are called dependency edges [1]. In our paper, interdependent directed networks contain two layers in terms of a layer A and a layer B, and each layer is a directed and unweighted scale-free network with the size N _A = N _B = N, as shown in Figure 1. Thus, the degree distributions of in-degree and out-degree are the power-law distribution in interdependent directed networks, where γ _in and γ _out are the power-law exponent of in-degree and out-degree, respectively. In Mathematics, it is sufficient to provide the adjacency matrix to formally characterize interdependent directed networks. For each layer (e.g., network A), and an asymmetric N × N matrix A whose generic entry a _ij = 1 if a link from node i to j exists, otherwise a _ij = 0.

One may observe cascades in interdependent directed networks, i.e., avalanches of failures triggered by the failure of one or more nodes, as the nodes are removed gradually with a specific order K. K is defined by K = max(k _A , k _B ), where the degree of nodes in the network A or B are set by k _A = k _A,in + k _A,out or k _B = k _B,in + k _B,out . In practice, we begin removing a fraction 1 − p in network A and a fraction 1 − p in network B, and removing all the links connected to these removed nodes. For interdependent nodes across layers, if node i fails to function due to being attacked or isolated, node i also fails in another layer. We continue this process until no further new failed nodes can occur.

To measure the robustness for interdependent directed networks under targeted attacks, we compute its mutually connected components (MCC) in each step of removing nodes with fraction 1 − p. Each layer network fragments into MCC, within which each pair of nodes can reach each other by a path [7, 20]. Some nodes in the MCC of the layer A network will play an important function in the layer A network, but they may not exist in the MCC of layer B. Thus, we define the MCC of interdependent systems to be the average value of all layers. A similar definition also applies to calculate the second maximum connected component (2nd-MCC). By doing this, when the reserved fraction p is tuned increasingly from zero to a unit, at a certain critical fraction p _c , the MCC of networks shifts from zero to non-zero. When p < p _c , the interdependent networks have no MCC, and otherwise p > p _c . The critical fraction p _c thus reveals the robustness of interdependent directed networks, i.e., the smaller p _c , the higher network robustness.

Intra-layer geometric correlations are composed of intra-layer radial correlations and intra-layer angular correlations in a certain layer, obtained by embedding interdependent directed networks into hyperbolic space. Therefore, we introduce the A-PSO (the asymmetric popularity and similarity optimization) model to map each layer of interdependent directed networks into hyperbolic space [24].

In this model, each node i is firstly split into two sets (a _i in the set a and b _i in the set b), and a directed link goes from a node i to a node j, which will be transformed a link between a _i and b _j . Then, each pair of nodes a _i and b _j correspond to polar coordinates ( θ a , i and r a , i ) and ( θ b , i and r b , i ) , respectively. The radial coordinates can be calculated by κ − r mapping: r = R − 2ln(κ/κ _min ), where hidden variable κ _*, _i (* ∈ {a, b}) is derived from ρ ( κ * ) = ( γ − 1 ) κ m i n ( γ − 1 ) κ * − γ , the minimum of hidden variable κ m i n = k ̄ ( γ * − 1 ) / ( γ * − 2 ) , and θ is drawn from uniform Probability Density Function (PDF). Finally, the directed link is created by any integrable function f ( χ ) = ( 1 + χ β ) − 1 in hyperbolic space, where hyperbolic distance χ = r _a,i + r _b,j + 2ln(d _ai,bj /2), β is a model parameter.

In practice, we do not know the nodes’ coordinates by given the adjacency matrix A of a layer. We are interested in the conditional probability P ( θ , κ | A ) that the possibility of assigning a coordinate to each node, giving our observed network data. Following Bayes’ rule, we have P ( { κ , θ } | a i j ) ∝ P ( a i j | { κ , θ } ) P ( { κ , θ } ) . (1) where the posterior distribution P ( { κ , θ } | a i j ) is proportional to two components: the likelihood P ( a i j | { κ , θ } ) of the network data a _ij , the prior probability P ( { κ , θ } ) κ and θ are obtained by following some constraints mentioned above, and we write P({κ, θ}) = 1 if the constraint is satisfied and otherwise P({κ, θ}) = 0. Thus, the likelihood can be calculated as followed: P ( a i j | { κ , θ } ) = ∏ 1 ≤ i ≠ j ≤ N f ( χ ) a i j [ 1 − f ( χ ) ] 1 − a i j , (2) where hidden variables are solved by κ _i = k _i − γ _/ β and angular coordinates are inferred by using the localized Metropolis-Hastings (LMH) algorithm [25, 26].

To this end, we have the angular coordinate (θ _a ; θ _b ) and the radial coordinates (r _a ; r _b ) in according with giving a directed network layer. Then, we use mutual information to describe intra-layer geometric correlations. Formally, the mutual information about two random various X, Y is obtained by [27]. I ( X ; Y ) = ∫ Y ∫ X p ( x , y ) l n p ( x , y ) p ( x ) p ( y ) d x d y , (3) where p(x, y) is the joint probability density function of X, Y, and p(x), p(y) are marginal PDF of X and Y. In this paper, the intra-layer angular correlation of each layer is quantified by the normalized mutual information NMI _θ = I(θ _a ; θ _b )/max{I(θ _a ; θ _a ), I(θ _b ; θ _b )}. Similarly, the intra-layer radial correlation is defined as NMI _r = I(r _a ; r _b )/max{I(r _a ; r _a ), I(r _b ; r _b )}. The higher the NMI (NMI ∈ [0, 1]), the stronger are the intra-layer geometric correlations.

We simulate targeted attacks on artificial networks to investigate the relationship between the robustness and intra-layer angular correlations. The geometric multiplex model (GMM, Ref. [19]) is applied to generate the artificial undirected interdependent networks with across-layer geometric correlations. Inspired by it, we use this framework to develop a single-layer directed network with an intra-layer geometric correlation. The difference between the GMM and our work is that we aim to obtain the out-direction and the in-direction coordinates in a specific correlation. In Figure 2, each directed layer is generated according to the following steps:

Step 1. Determine the initial parameters: the network size N, the exponent β of the connection probability, the power-law exponent of out-degree γ _a , the power-law exponent of in-degree γ _b , the average degree k ̄ , the intra-layer radial correlation v ∈ [0, 1], and the intra-layer angular correlation g ∈ [0, 1].

Step 2. Determine the hyperbolic coordinates with a certain correlation in each layer of interdependent directed networks.

First of all, each node is assigned in-direction hidden variables κ _b , θ _b in the set b, as sampled from ρ ( κ b ) = ( γ b − 1 ) κ b , m i n ( γ b − 1 ) κ b − γ b and uniform PDF, respectively.

Secondly, the out-direction angular coordinates are chosen from θ _a = mod[θ _b + 2πl _i /N, 2π], where l _i is an arc length of radius R in a hyperbolic disc, which is satisfied by zero-mean truncated Gaussian PDF, defined as f σ = 1 / σ ϕ ( l / σ ) Φ ( N / 2 σ ) − Φ ( − N / 2 σ ) , σ = σ ₀(1/g − 1), σ ₀ = min[100, N/4π]. ϕ(x) is normal distribution. Φ(x) is the PDF of ϕ(x).

Thirdly, each node of out-direction radial coordinates r _a is assigned. Notice that r _a is taken the place of the hidden variables κ _a to implement the algorithm easily. Specifically, the κ _a is derived from the copulas function C η ( F ( κ a ) , F ( κ b ) ) = e − [ − l n ( F ( κ b ) ) η + ( − l n ( F κ a ) ) η ] 1 / η , where F ( κ a ) = 1 − κ a ( 1 − γ a ) κ a , m i n ( γ a − 1 ) , F ( κ b ) = 1 − κ b ( 1 − γ b ) κ b , m i n ( γ b − 1 ) , η = 1/(1 − v), the minimum of hidden variable κ a , m i n = k ̄ ( γ a − 1 ) / ( γ a − 2 ) and κ b , m i n = k ̄ ( γ b − 1 ) / ( γ b − 2 ) . We then transform hidden variables to radial coordinates r _a = R − 2ln(κ _a /κ _a,min ) and r _b = R − 2ln(κ _b /κ _b,min ).

In particular, when g = 1 and v = 1, the coordinates of each node is identical in the two directions (that is, θ _a = θ _b and κ _a = κ _b , respectively), and a generated network degenerates into a undirected network. To overcome this problem, we regard v = 0.99 and g = 0.99 as the full correlations (i.e., v = 1 and g = 1) in this paper and make sure to generate a directed network.

Step 3. Determine the artificial networks. Links are created by the connection probability, i.e., each node pair i, j is connected by probabilities f ( χ i j ) = ( 1 + χ β ) − 1 , χ i j = Δ θ N k ̄ β s i n ( π / β ) κ i a κ j b and Δθ = |π − |π − |θ _ia − θ _jb ‖|. To do so, a layer of the artificial network has been constructed. The steps 1–3 are repeated to generate another layer.

In this section, all artificial networks are double-layer directed networks, where a pair of nodes across layers are interdependent. The artificial geometric model is used to generate each layer which is a heterogeneous directed network with the power-law exponent of out-degree γ _a = 2.6, the power-law exponent of in-degree γ _b = 2.6, average node degree k ̄ = 6 , and the parameter β = 3.5. The intra-layer angular correlation and the intra-layer radial correlation are denoted by the symbols ( g _A , v _A ) in layer A and ( g _B , v _B ) in layer B.

By doing this, three kinds of artificial networks have been generated to simulate targeted attacks, as shown in Figures 3A,B. Results reveal two geometric contributions to the robustness. One, the value of p _c in the orange line is larger than others, which shows that intra-layer angular correlations increase the vulnerability of interdependent directed networks. Notice that our results are in contrast to the situation on across-layer correlations between interdependent networks [20], which reveals that intra-layer angular correlations are hidden factors to understand complex systems. Two, such vulnerability will be exacerbated by the increase in the number of layers. Additionally, we also provide the behavior of cascading failures for the 2nd-MCC in different size systems. The largest 2nd-MCC achieves its extremum near the critical point, which is a way to estimate and compare p _c . Figure 3C illustrates the extreme value point p _c for interdependent directed networks with full angular correlations is always significantly higher than the case of the non-angular correlations. Multi-subsystem interaction and its hidden geometric structure thus should be considered designing network systems more robust. Additional, we analyze the fraction of coupling strength q ∈ [0, 1], where q = 0 represents that network systems become two single and independent networks, and q = 1 represents the mapping relationship of nodes between two layers is one to one. Figures 3D–F illustrates that the results of Figure 3A can expand to general cases (the inter-layer coupling of arbitrary proportions). Figure 3F shows their percolation behaviors with the fraction of remaining nodes p changing from 0 to 1 under intra-layer angular correlations and different coupling strengths q. The results show that decreasing coupling strength can mitigate the vulnerability of interdependent directed networks with the intra-layer angular correlation against targeted attacks.

To study this issue further, we examine the impact of different angular correlations on robustness by several variations of artificial networks, as shown in Figures 3G,H. As intra-layer angular correlations decrease, the vulnerability of directed systems is mitigated, irrespective of the effect of intra-layer radial correlations. This means that, although the contribution of in-out degree correlations is positive to robustness for interdependent directed networks, intra-layer angular correlations play an essential factor in undermining the robustness. Thus, intra-layer angular correlations have an early-warming function when interdependent directed networks face a sudden extreme attack. In addition, we also found that such an increasing trend is not apparent in low-correlation situations. To analyze the cause, we checked the relationship between parameter g and the Pearson correlation between intra-layer angular coordinates. Figure 3I suggests that the nonlinear relationship induces the phenomenon mentioned above.

The influence of intra-layer angular correlations on the robustness in real-world interdependent networks is simulated in this subsection. Empirical networks are all derived from open databases and describe in detail, as followed. 1) C. elegans neural dataset describes the neural interconnection via chemical synapses and gap junctions, which can be obtained from the Wormatlas database [28]. The nodes are neurons, and each layer corresponds to a different type of synaptic connection. 2) International trade dataset considers different types of trade relationships among countries, obtained from Ref. [29]. The worldwide food import/export network is an economic network in which layers represent products, nodes are countries, and edges at each layer represent import/export relationships of a specific food product among countries. Each layer is directed and weighted networks with 214 nodes. 3) Arabidopsis interdependent Genetic networks are obtained from the Biological General Repository for Interaction Datasets (BioGRID, thebiogrid.org), a public database that archives and disseminates genetic and protein interaction data from humans and model organisms [30, 31]. Each layer is directed and unweighted networks with 1,449 nodes after removing the isolated nodes. 4) Social networks consist of 3 kinds of (Co-work, Friendship, and Advice) between partners and associates of a corporate law partnership [32, 33]. Each layer is directed and unweighted networks.

Secondly, we construct reshuffled counterparts (so-called reshuffled networks) from real-world networks (so-called original networks). The reshuffled counterpart is a variant of the original network to alter intra-layer geometric correlations. Specifically, each layer (a directed network) is transformed into a bipartite structure, as shown in Figure 4, and randomly reshuffled nodes’ ID of the set b in a way. Notably, interdependent nodes are also reshuffled in the same way in other layers if nodes’ ID is reshuffled at a layer. Node b ₁ and node b ₄ are also reshuffled in layer B, when node b ₁ and node b ₄ are reshuffled in layer A, To this end, the reshuffled counterparts are destroyed the intra-layer geometric correlation and preserved across-layers geometric correlations.

Then, each layer for these networks can be embedded into a hyperbolic space, where each layer is represented by a group of angular coordinates (θ _a , θ _b ) and a group of radial coordinates (r _a , r _b ). To validate the influence of intra-layer geometric correlations on this real-world multiplex network, we implement targeted attacks on the original networks and reshuffled networks for empirical networks, respectively. Figure 5 displays that the p _c of original networks is smaller than their reshuffled counterparts under targeted attacks. We also observe that the NMI _θ,ori is always larger than the NMI _θ,re for different real-world networks, as shown in Table 1. Linking those results, we find that the larger the value p _c , the stronger intra-layer angular correlations. It is the reason why interdependent directed networks are more robust after the reshuffle. Results also suggest our arguments and highlight the importance of intra-layer angular correlations.