Cascading failures in anisotropic interdependent networks of spatial modular structures

The structure of real-world multilayer infrastructure systems usually exhibits anisotropy due to constraints of the embedding space. For example, geographical features like mountains, rivers and shores influence the architecture of critical infrastructure networks. Moreover, such spatial networks are often non-homogeneous but rather have a modular structure with dense connections within communities and sparse connections between neighboring communities. When the networks of the different layers are interdependent, local failures and attacks may propagate throughout the system. Here we study the robustness of spatial interdependent networks which are both anisotropic and heterogeneous. We also evaluate the effect of localized attacks having different geometrical shapes. We find that anisotropic networks are more robust against localized attacks and that anisotropic attacks, surprisingly, even on isotropic structures, are more effective than isotropic attacks.


I. INTRODUCTION
Many real-world systems are well correlated with a population density that is highly non-uniform and concentrates in large cities or is spread along seacoasts, rivers, or major transportation routes. Such systems are influenced by geographical and social features and usually combine both spaciality on large scales and randomness on small scales. More specifically, in many cases such as infrastructure networks, the connections within the cities are dense and almost uniform, while the connections between cities are mainly between nearby cities. These factors were the motivation behind the profound studies of spatial networks [1][2][3][4][5][6][7][8][9][10], either homogeneous or heterogeneous and composed of communities. However, the connections between the cities do not have to be isotropic. For example, if in one axis there are more topographic obstacles (e.g. mountains or rivers) then in that axis there will be fewer connections than in the other axes. Yet, most models of spatial networks are isotropic, i.e, have no preferred orientation in the network structure.
In addition, many realistic systems such as modern infrastructure networks can be considered as strongly interdependent multi-layered networks [11][12][13]. Therefore, great consideration was given to the study of interdependent networks and in particular their robustness, using percolation theory [14][15][16][17][18][19][20][21][22][23][24][25]. Such studies have been carried out to analyze cascading failures resulting from various attacks [26][27][28][29][30][31][32]. One of the key results in these types of studies is that often localized attacks are significantly more effective when compare to random attacks. To the best of our knowledge, in studies of localized attacks, the attack itself has always been uniform within a certain radius [33][34][35][36] although in the real world this "uniformity condition" usually does not occur. For example, natural The nodes are represented by gray circles. The connectivity links can be divided into three types, as follows: 'intra-links' which connect nodes within the communities (black links), 'horizontal-links' (orange) and 'vertical-links' (blue) which connect pairs of nodes from neighboring communities in different orientations. Nodes from the two layers are interdependent (dashed purple link). (b) The system is constructed as a m × m lattice of communities, each of which is an ER network of nodes located at the sites of the ζ × ζ square lattice. For clarity here we show 3 ER networks of size ζ = 3. In the simulations we set periodic boundary conditions and we solve the limit of infinitely large ER communities, i.e., the case of ζ → ∞. Anisotropy is modeled by different degree distributions of horizontal-links and vertical-links (three orange links and two blue links).
disasters (such as an earthquake) and malicious targeted attacks can be anisotropic and even contained within a single axis. Here, we take into consideration anisotropy in modeling both in the systems and in the localized attacks.

II. MODEL
Here we introduce an anisotropic model of interdependent spatially embedded networks with communities. The model assumes that the entire 2D territory is divided into squared communities of size ζ × ζ representing cities or densely populated areas, see Fig. 1. For simplicity, we assume a multiplex system with two topologically distinct network layers that share the same set of nodes. The interdependence between the layers is introduced as follows: if a node becomes dysfunctional in one layer it is also dysfunctional in the other layer. In each layer, the connectivity links within a community ('intra-links') are connected at random, i.e. forming ER networks, while the links connecting nodes in two different communities ('inter-links') can only connect neighboring squares in each network, horizontally or vertically, as illustrated in Fig. 1(b). Each node has a degree k intra of intra-links, a degree k H of horizontal-links, a degree k V of verticallinks, and the total degree is k total = k intra +k H +k V . We assume that k intra , k H and k V are independent random variables taken from three different degree distributions which are characterized by average degrees k intra , k H and k V , respectively. The anisotropy of the system is specified by the parameter which is the ratio between the degree of the horizontallinks and the inter-links (and thus γ = 1/2 is the isotropic case), and the heterogeneity of the system is specified by the ratio between the degree of the inter-link and the total degree α = ( k H + k V )/ k total .
Here we study the robustness of spatial anisotropic multiplex networks to different forms of localized attacks. Initial damage in a multiplex network spreads in a process in which failures of some nodes lead to failures of other nodes and so on. In the cascading failures process, a node fails if it is no longer connected to the giant component in one of the layers. Thus, the cascading stops when all the remaining nodes are part of the mutual giant component (MGC), i.e., all the remaining nodes are connected to the giant component in both layers. We consider a network as functioning if the MGC is of size O(N ), where N is the number of nodes in the network. Accordingly, in our spatial multilayer network model, the network is considered to be resistant to a localized attack if the damage spreads to a finite number of communities and does not reach the edges of the system.

III. ANALYTICAL APPROACH
In our recent study [36], we developed a mathematical framework that consists of general equations for the MGC size of a multiplex model of m × m ER communities connected as a lattice. However, in that study, we applied the equations only for the case of an isotropic network and we analyzed only the effect of isotropic localized attacks on the functionality of the network. Here, we study a more general and realistic case, of an anisotropic model with different average degrees for the vertical links k V and the horizontal links k H , as well as anisotropic localized attacks. For a network with a given levels of heterogeneity and anisotropy, we determine whether it remains functional after the cascading failure induced by isotropic or anisotropic initial damages. The initial damage in community i (for i from 1 to m 2 ) is expressed by the parameter p i , which is defined as the fraction of nodes that survived as a result of the damage. In the analytical model we consider infinitely large ER communities (ζ → ∞) that are characterized by Poisson degree distribution. Using the analytical formalism from our recent study [36], we obtain the following equation for the MGC size of each community i, P ∞,i , where Next we use this analytical formalism (Eq. (2)) to evaluate the robustness of interdependent spatial systems characterized by heterogeneity and anisotropy by calculating the steady state of each community. In the Appendix, we verify these analytical solutions through simulations. Fig. 2 shows that for a given k total , the robustness of the network to a removal of one community is highly dependent on both parameters α and γ representing the heterogeneity and anisotropy, respectively. In particular, in the phase diagrams for k total = 2.48, 2.49, 2.5 (Figs. 2(a)-(c)), the steady state of the system can be in one of two extreme states: stable (in yellow) -in which the system remains functional, and unstable (in blue)in which the entire system collapses. It is important to note that the regions of the stable and unstable states do not depend on the system size m (see also Fig. 7 in the Appendix). In addition, the transition from stable to unstable is non-monotonic with the anisotropy of the network. For instance, for k total = 2.5 and α = 0.5 (see Fig.2(c)), the unstable region is only in a specific narrow region of γ that is located between the extreme cases of γ = 0.5 (isotropic network) and γ = 0 or γ = 1 (highly anistortopic network).

IV. RESULTS
The blue and orange curves added to the phase diagrams (Figs. 2(a)-(c)) show that the unstable region is contained within a specific area which is between the two (a)-(c) Phase diagrams of the size of the functioning component, P∞, at steady state after initial removal of one community from a system of 441 communities (i.e. m = 21). The color of each pixel represents the relative size of the system's functioning component. In all phase diagrams we show the analytic solutions for k total − kV /2 = kc (blue curve) and k total − kH /2 = kc (orange curve), where kc = 2.4554 [14]. (d) The network at steady state for the case of α = 0.8, γ = 0.15 and k total = 2.5. Each pixel represents a community, and the color of the pixel represents the relative size of the functioning component of the community. Since γ < 1/2, the network is anisotropic with more vertical-links than horizontal-links (see Eq. (1)). The initial damage of removing one community causes cascading failures which spread more in the vertical axis than in the horizontal axis, i.e., along the direction of the higher degree.
curves and to the right of their intersection. This area is defined by the following constraints: k total − k H /2 < k c and k total − k V /2 < k c , where k c ≈ 2.4554 is the critical average degree below which a single ER multiplex collapses without any initial damage [14]. These constraints describe a case when removing a community causes its neighboring communities to reduce their average total degree below k c . In this mentioned area, there is also a stable region because the communities are not isolated and therefore a community with a total degree smaller than k c can be sustained by the communities next to it. In Fig. 2(d) we show the behavior of cascading failures in this stable region. We show that the initial damage of removing one community spreads further along the direction of the higher inter-degree.
Next, we determine the robustness of the system to various forms of isotropic and anisotropic localized attacks. We find that our system is metastable for a broad range of parameters, meaning that for a localized attack of a given shape there is a critical size above which the induced cascade of failures propagates through the whole system leading to its collapse. We observe that the size of the critical attack, which does not depend on the num- We find a metastable region, where a finite-size localized attack larger than r c h , which does not depend on the system size m, causes cascading failures and system collapse. The metastable region reduces as k total increases, and for a fixed k total , there exists a critical value of α above which the damage of finite radius causes the collapse of the network. The critical α increases both by k total and anisotropy.
ber of communities m, is significantly different for various shapes of attacks.
In Fig. 3 we analyze the case of isotropic localized attacks in a form of a circle with a radius r h . We show that for a given k total and γ, there is a critical α c which has the following properties. For heterogeneity of α < α c the critical size r c h is of size ≈ m/2 (system size), and for α > α c , r c h does not depend on the system size and may contain more than one community. In Fig. 4(a) we present the critical attack size for attacks in a form of a strip, i.e, removing communities connected in a row. We find that the attack is more efficient when the removed strip is orthogonal to the direction of the higher degree so it destroys more interlinks per one removed community. In addition, we find that increasing the anisotropy of the network causes the critical length of the stripe attack to The critical size (in units of a community size) for an attack in a form of a strip with a thickness of 1 community. For γ < 1/2 the strip is orthogonal to the direction of the higher degree interlinks and for γ > 1/2 the strip is parallel to that direction. This critical size increases with γ and thus it can be concluded that for an efficient attack on an anisotropic network it is more advisable to remove a strip perpendicular to the direction of the higher degree interlinks. For these simulations m = 21 and k total = 2.58. (b) Demonstration of the evolution of cascading failures after an attack of length 10, which is larger than the critical size. Each pixel represents a community and the color represents the time t when the community completely collapsed. For this simulation m = 100, k total = 2.58, α = 0.8 and γ = 0. 3. decrease. Note that we have chosen to focus on the range of anisotropy 0.3 < γ < 0.7 because outside this range the system is always stable, i.e. the critical attack size is in the size of the system (see Fig. 3). Note also that the small fluctuations are caused by the fact that the network is not continuous but consists of discrete communities, i.e., there are jumps when the damage crosses neighboring communities.
In addition, in Fig. 4(b), we analyze and demonstrate the evolution of cascading failures after a strip attack whose length is above the critical size. For these simulations, we solved the general iterative equations from [36], which represent the stages of the cascading failures, for our specific anisotropic model (see Eq. (5) in the Appendix). We find that the damage propagates first in the axis with the higher inter degree and only after the damage reaches a certain distance in this axis, it begins to propagate also in the axis with the lower inter degree. From this we conclude that the cascading along the axis is affected more by the size of the front of the attack and less by the depth of the attack.
Finally, we analyze the robustness of our anisotropic model with respect to localized attacks in a form of an ellipse with different flattening f e , defined as f e = (a − b)/a where a and b are the semi-major and semi-minor axes, respectively. The case f e = 0 corresponds to our earlier circular attack while f e → 1 represents a strip-like attack. Thus, the general ellipse case covers a wide range of anisotropic attacks between a strip attack and a circle attack. In Fig. 5 we show that the ellipse attack is more efficient for ellipses with higher flattening f e , and that a strip attack is more efficient than any ellipse attack. We therefore conclude that for a given area of the attack, it is more efficient that the shape of the attack will be as anisotropic as possible.

V. DISCUSSION
Many spatial networks are influenced by geographical features that can lead to an anisotropic structure in which the amount of links differs in different directions. In this work, we have introduced a new realistic model, which considers, for the first time, the aspect of network anisotropy. Using tools from percolation theory, we systematically study anisotropic systems and analyze their robustness to various isotropic and anisotropic localized attacks. We determine how the robustness of our network model depends on its heterogeneity and anisotropy. We also show that anisotropic attacks reveal significantly increased vulnerability compared to the considered earlier, simple circle-shaped attacks. Specifically, we find that a localized attack causes larger damage if the area of the attack is more elongated in the direction of the smaller inter degree, i.e. when the attack cuts more links that are in the direction of the higher inter degree. We find that even in isotropic networks the anisotropic localized attacks are more efficient than isotropic attacks, see the case γ = 0.5 in Fig. 5. In our future work, we will determine how the velocity of the cascading failures in the different axes depends on the parameters of the anisotropic network, and will try to answer the question of whether this velocity exhibits a scaling behavior.  In Fig. 6, we compare between the theory Eq. (2) and numerical simulations for P ∞ after initial removal of one community and observe excellent agreement between them.

B. The cascading failures equations
Here we present, for the convenience, the general equations from our previous article [36], and apply them to the anisotropic model of this manuscript. For a multiplex model with two layers A and B, the vector equations of the cascading failures starting from t = 0 are: where the parameters are defined as follows, (1) f i,j is the probability that a randomly selected link, that passes from a node in community i to a node in community j, does not lead to the giant component (GC), (2) g i is the fraction of nodes in community i which belong to the GC, and (3) p(t) is the fraction of survived nodes at stage t of the cascade of failures. For our anisotropic multiplex model, in which the degree distributions are the same in both layers, and in the limit of infinitely large ER communities, the functions Φ and Ψ are obtained by f j = Φ j ( f , p) ≡ (1 − p j ) + p j · e i ki,j ·(fi−1) g j = Ψ j ( f , p) ≡ p j (1 − e i ki,j ·(fi−1) ), where f j ≡ f i,j and k i,j is defined in Eq.(3) in the main text. Figures 7 and 8 show that there is no effect of the system size m on our analytical results. The critical size attack for attacks in an ellipse form and for removing a strip, on a log-linear graph, for a wide range of γ for k total = 2.58 and α = 0.8 (as in Fig.  5(b)). The continuous lines represent m = 21 and the dashed lines represent m = 41. For γ < 0.25 the system size m affects the size of the critical attack, which is the whole system. In contrast, for γ ≥ 0.25 the critical size does not depend on the system size m.