Spreading of localized attacks in spatial multiplex networks

Many real-world multilayer systems such as critical infrastructure are interdependent and embedded in space with links of a characteristic length. They are also vulnerable to localized attacks or failures, such as terrorist attacks or natural catastrophes, which affect all nodes within a given radius. Here we study the effects of localized attacks on spatial multiplex networks of two layers. We find a metastable region where a localized attack larger than a critical size induces a nucleation transition as a cascade of failures spreads throughout the system, leading to its collapse. We develop a theory to predict the critical attack size and find that it exhibits novel scaling behavior. We further find that localized attacks in these multiplex systems can induce a previously unobserved combination of random and spatial cascades. Our results demonstrate important vulnerabilities in real-world interdependent networks and show new theoretical features of spatial networks.

Here, we study localized attacks on a realistic spatial multiplex model that has been proposed recently [47,48]. The system is a model of multiplex with exponential link-length distribution of connectivity links in each of the two layers: Here ζ is a parameter determining the characteristic link length and thereby the strength of the embedding-a smaller ζ reflects a stronger embedding. Networks with links of characteristic length ζ appear in reality, for example, the European power grid and the inter station local railway lines in Japan [47,[49][50][51]. We further assume that the nodes require connectivity in each layer in order to function, a requirement which is equivalent to having dependency links of length zero with longer connectivity links. This is in contrast to the research based Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
on the model of Li et al [31] and Berezin et al [39] which considered the case where dependency links are longer than connectivity links. We suggest that the assumption of dependency links which are shorter than connectivity links is more natural, because, for example, it is more likely for a communication's station to receive power from its nearest power station than a distant one, though the communications and power networks are known to have potentially long links [47,52]. As we show here, the combination of spatially constrained connectivity links and multiplex dependencyboth ubiquitous features of real complex systems-makes these systems vulnerable to potentially catastrophic localized attacks. Such attacks are important and realistic because they can represent a local damage on two spatial networks that depend on one another to function in a very natural way: the nodes are either the same, or every node in one network layer depend on a close node in the other.
We find that for a broad range of parameters our system is metastable, meaning that a localized attack larger than a critical size-that is independent of the system size-induces a cascade of failures which propagates through the whole system leading to its collapse. We develop a theory which can predict this critical size of the initial local damage, and can explain the unique cascading process that makes the critical size independent of the system size. We find that when the localized attack is of the critical size-the cascade is at first random within a disc of radius of order ζ, and then it propagate spatially until it reaches the boundaries of the system. Using this theory we also find a new scaling exponent describing the critical nucleation (of damage) size.

Model
We model the multiplex composed of two layers in which the nodes are placed at lattice sites of a square lattice where the link lengths r are distributed with probability of equation (1) and average degree á ñ k (see figure 1). Here, we focus on the case in which both layers have the same characteristic length ζ and same á ñ k . In practice, we assign each node an (x, y) coordinate with integers , and construct the links in each layer as follows: (a) we randomly select a source node ( ) x y , s s and draw an angle α selected uniformly at random. (b) We draw a length r selected from the distribution P(r), equation (1). (c) We select the target node ( ) x y , t t , which is closest to satisfying, . This process is executed independently in each layer and is continued until we have a total of á ñ N k 2 links. The topological model is similar to the Waxman model [53] and recent work by Bianconi and Halu et al [24,54] with a key difference being that our model converges to a lattice as z  0.
For a node to remain functional it must be connected to the giant component in both layers. This reflects the assumption that in order for the node to continue to function it requires the two types of connectivity. Next, we perform a localized attack as follows: (a) we remove all nodes within a distance r h from a random location in the system.   (1)) with characteristic length z = 3 and are connected at random.
At the end of this cascade, the system is categorized as functional or non-functional depending on whether the MGC is of the order of the system size L 2 or not.

Results
We analyzed the damage spreading of the localized attack on the multiplex with different á ñ k , ζ and r h . Our simulations suggest the existence of r h c , a minimum radius of damage needed to cause the system to collapse. Below r h c the damage remains localized while for a radius above r h c the damage propagates indefinitely and destroys the whole multiplex. When we calculate the critical attack size r h c for different á ñ k and ζ, we discover three regions, as shown in the phase diagram in figure 2. The regions are: (a) stable (in red)-in this region the system remains functional after a localized attack of any finite size. (b) Unstable (in blue)-in this region the system is non-functional even if no nodes are removed. (c) Metastable (between the above-mentioned regions) -in this region only attacks with radius larger than r h c propagate, through cascading failures, the entire system and makes it non-functional.
To understand these phenomena we consider the network as being composed of regions of size of order ζ that are tiled on a 2D lattice, each of which can be approximated as a random network. The localized attack of size r h can then be approximated as a random attack of size pr h 2 in an interdependent random network with z 2 nodes. In this case, the percolation threshold for random removals is known to be = á ñ p c k k c , where » k 2.4554 c [18]. It has also been shown that localized attacks (formed by shells surrounding a root node) in Erdős-Rényi multiplex networks have the same percolation threshold as random attacks [40,41]. Based on this, we can predict the critical attack size r h c close to k c as follows: where a is the constant of proportionality for the effective random network (radius) size, which we determine numerically. This r h c is the minimal expected size of the hole that destroys the entire random network regime ( z a ). However, since there are links between the tiled Erdős-Rényi sub-networks, the collapse propagates toward the surrounding sub-networks and we see a typical spreading cascade in an embedded network.
For the limit of ζ of the order of L, the multiplex can be well approximated as two interdependent Erdős-Rényi networks, and therefore we can calculate r h c as follows: from which, We show that equations (2) and (4) predict the simulation results in figure 3(a) with » a 9. Because of the long links and since the sub-networks are not isolated-a is relatively big. In supplementary section II is available online at stacks.iop.org/NJP/19/073037/mmedia we can see similar phenomenon on a system with average degree á ñ k with links that connected slightly different. In this alternative model we choose a node randomly and link it to another node with link-length distribution of step function up to ζ, for each of the two layers. In this case, there are no long links (but the small Erdős-Rényi networks are still not isolated) so a is found to be smaller than in our model (approximately 3.2).
For multiplex networks, near criticality, ¥ P (the size of the MGC) fulfill the scaling~á ñ -b ¥ ( ) P k k c (in lattice for example b = 5 36 see e.g. [57]). In our model, in analogy to ¥ P , we find theoretically (equations (3) and (5)) that r h c scales as á ñ -( ) k k c 1 2 , suggesting that 1 2 is a critical exponent for r h c . Indeed the simulations shown in figures 3(b) and (c) support this exponent. Generally, it is difficult to find evidence for universality in the absence of a second-order transition. This new scaling, related to nucleation type processes, may provide an alternative approach which can be useful to understand universality properties in critical phenomena associated with a first-order transition where nucleation processes are involved [58,59]. We also find a new dynamical process of cascading when the localized attack is near the critical size, that is consistent with our theory. To understand this process for a given á ñ k and ζ, we follow the standard cascade process until the MGC reaches a steady state. At this time (which we call t a ), we remove a hole with radius r h c which initiates a new cascade. Figure 4(a) shows the whole spatial-temporal process of cascading and figures 4(b) and (c) demonstrate the different types of number of iterations (NOI) in the two regimes as described below. The graph in figure 5(a), of á ñ r , the average distance from the center of the nodes that failed in every iteration, reveals explicitly the three main stages of the whole process shown in figure 4(a): (i) before the localized attack (until the dashed line at t a ), there are a few steps where the cascade describes the removal of nodes that are not in the MGC,  (3) and (5)). For this figure L=2000 with averages over at least five runs for each data point. so á ñ r is close to the average distance from the center to all nodes (∼1500). (ii) Random branching process [60,61] in limited annulus around the hole (demonstrated in figure 4(b)) so á ñ r is fixed for many iterations at distance z » · a 2 . (iii) Spatial spreading process that propagates the whole system (demonstrated in figure 4(c)), so á ñ r increases linearly as a function of NOI. Indeed we can see the effect of the three above processes in figure 5(b)-the size of the MGC, ¥ P , at first decreases sharply, then, after the attack in t a , it decreases very slowly in a plateau, and then parabolically as a function of NOI. The transition from phase (ii) to (iii) can be discerned by identifying a transition in á ñ r from constant to linear increase ( figure 5(a)), or from a transition in ¥ P from approximately constant to parabolically decreasing (figures 5(b) and (c)). Additionally, the processes are also described in supplementary section I in the discussion about the branching factor.
In figure 5, we see that the cascade begins random-like, with no spatial influence within the neighborhood of the failure and a random branching process with expected branching factor of »1, as established for interdependent random networks [28,60,62]. However, this random-like behavior is constrained to the neighborhood of radius z a . Once the damage spreads beyond this neighborhood, it expands linearly in space, with a constant rate and a parabolic decrease in ¥ P , as documented for spatially embedded interdependent networks [31,32,35]. A similar coexistence of random and spatial properties, differentiated by scale, has been observed in the single-layer case [48,63].
In the spreading process we can see the cascading dynamics both in the simulations for á ñ( ) r t and in the simulations for ¶ ¶ ¥ ( ) P t t in figure 5(c). The connection between them is expressed in the equations below so that t expresses the NOI and v, that sets the speed of the cascading, is z0.6 , Understanding the dynamical process of cascading can explain why in the metastable region, when the size of the network crosses the size of our approximated random network (around point z » L a 2 in figure 6), there is no correlation between the critical attack size r h c and the system size. This is because once the network is large enough for a damage spreading process to take place, the hole will spread until the damage reaches the edges of the system, regardless of its size.

Discussion
We have presented a study of interdependent spatial networks with a novel and realistic combination of spatially localized damage and connectivity links which are longer than the dependency links. This combination is ubiquitous in nature, and yet has not been studied methodically, to our knowledge. We find that a nucleation phenomenon can be triggered by local damage, with failures spreading through the entire system. The cascade itself has random behavior on a small scale but spatial behavior on a large scale, similar to what has been observed in the single-layer case [48,63]. We further find that the critical nucleation size has novel scaling features. Future research will determine whether this indicates a general, universal feature of nucleation transitions. Figure 6. Dependence of the critical attack size r h c on the system size L. We see that above a certain value of L the critical attack size r h c is constant. For this figure á ñ = k 2.5, with 5 runs for each data point.