Spreading of localized attacks on spatial multiplex networks with a community structure

We study the effect of localized attacks on a multiplex spatial network, where each layer is a network of communities. The system is considered functional when the nodes belong to the giant component in all the multiplex layers. The communities are of linear size $\zeta$, such that within them any pair of nodes are linked with same probability, and additionally nodes in nearby communities are linked with a different (typically smaller) probability. This model can represent an interdependent infrastructure system of cities where within the city there are many links while between cities there are fewer links. We develop an analytical method, similar to the finite element method applied to a network with communities, and verify our analytical results by simulations. We find, both by simulation and theory, that for different parameters of connectivity and spatiality --- there is a critical localized size of damage above which it will spread and the entire system will collapse.

We study the effect of localized attacks on a multiplex spatial network, where each layer is a network of communities.The system is considered functional when the nodes belong to the giant component in all the multiplex layers.The communities are of linear size ζ, such that within them any pair of nodes are linked with same probability, and additionally nodes in nearby communities are linked with a different (typically smaller) probability.This model can represent an interdependent infrastructure system of cities where within the city there are many links while between cities there are fewer links.We develop an analytical method, similar to the finite element method applied to a network with communities, and verify our analytical results by simulations.We find, both by simulation and theory, that for different parameters of connectivity and spatiality -there is a critical localized size of damage above which it will spread and the entire system will collapse.
In recent years, due to the advances in technology, many systems have become more and more integrated and interdependent.This dependence between these systems can cause a spread of damages, and lead to a cascade of failures and even entire system collapse.Therefore, many studies have been carried out to analyze cascading failures in interdependent networks [1][2][3][4][5][6][7].Further, in many real systems such as power grids and transportation systems, the links are of typical relatively short length due to the embedding in space.In such spatial systems, the initial failures or attacks can be localized to a specific region.Recent studies show that in different cases of spatial interdependent networks, localized attacks are significantly more damaging than random attacks [8][9][10][11].In addition, many real networks have a modular structure, such as biological networks [12] and many infrastructure systems [13,14].Therefore, recent studies have explored and compared the robustness of individual and interdependent modular non-spatial systems [15][16][17][18].Our study combines for the first time, three ubiquitous features of real complex systems -interdependence, spatiality and modularity.
Here, we analyze and predict the resilience of spatial multiplex networks with modular structure, see Fig. 1, under localized attacks, by developing tools based on percolation theory.An example of a realistic system that motivates our model is the infrastructure networks in a country, where each layer describes different infrastructure.The different infrastructures are dependent on each other, and in addition, each layer has high connections within the cities and a few long connections between nearby cities.We focus on localized failures because of two main reasons.First, a localized damage is a realistic scenario (due to flood or earthquake), and second, in such systems, a finite number of local failures concentrated in the same area might spread the damage throughout the system and cause significant damage and even to fully system collapse.with spatial and community properties, see Fig. 1.A multiplex network is a single network with at least two kinds of connectivity links.We assume here that the two types of links serve for two different functions, such as transportation and communication.In fact, a multiplex network with two kinds of connectivity links (for instance) can be regarded as a special case of interdependent networks with two layers with the same number of nodes, and every node in one layer has only one interdependent link with a single node in the other layer.For a node to remain functional in the multiplex, after the percolation process [19,20], it must be connected to the giant component in both layers.This reflects the assump-FIG.2: The critical attack size r c h as a function of α for different k total values.For every k total the lines represent the theory of Eq. ( 10) and the symbols represent simulation results.As seen, the simulations show excellent agreement with the theory.For a given k total , there is a metastable region that starts above a certain critical αc.In this region, r c h is smaller then the system size and increases only slightly with α.The figure suggests that for a certain k total one needs α above a critical αc to enable propagation of a localized damage.Below αc one needs to remove order of L 2 nodes for the system to collapse.Above αc a finite number of nodes (zero fraction) leads to system collapse.For the simulations we set L = 2100 and ζ = 100, with average over 5 runs for each data point.

Model. Our model is generated as a multiplex system
tion that in order for a node in the system to functionit requires both resources provided by the two layers.
Our multiplex model is composed here, for simplicity and without loss of generalization, of two layers in which the nodes are placed at sites of a square lattice of size L × L. The multiplex is constructed as m × m Erdős-Rényi (ER) multiplex networks (communities), each of which of size ζ × ζ, that are tiled and connected to each other as a square lattice (see Fig. 1).We assume, that the links within a community (intra-links) are connected at random, while the links connecting nodes in two distinct communities (inter-links) can only connect neighboring tiles.Each node has a degree k inter of inter-links and a degree k intra of intra-links.The total degree is: k total = k inter + k intra .In addition, the heterogeneity of the system is specified by the interconnectivity parameter α = k inter / k total .It should be noted that the homogeneous case (without communities) has been previously studied both for single-layer [21,22] and multi-layer [11,23] networks.In that model, all links have a characteristic length ζ with no distinction between inter and intra links, and therefore representing an homogeneous systems, without communities.In contrast, the present new model can describe systems with a spatial structure of communities, where the heterogeneity of the system is controlled by the α parameter.Thus, this model enables us to expand the previous model to a more general and realistic one for systems such as interconnected cities.
Theory and Simulation Results.Here we develop a framework for understanding the cascading process for a general case of a multiplex network.The constrains are that each layer has a spatial structure of communities that are connected in a form of a lattice.The cascading failures in multiplex networks is a process in which failures of some nodes lead to failures of other nodes and so on.When the cascading failures stops we are left with the mutual giant component (M GC) of the network.In the cascading process, at first we remove all nodes that are not in the giant component (GC) of the first layer.Then, from the set of the remaining nodes, we remove all the nodes that are not in the GC of the second layer.We repeat this two steps until there are no nodes to remove, and we are left with the M GC.The existence of a M GC of volume O(N ) where N is the number of nodes in the network, expresses the functionality of the system.
In order to obtain analytically the size of the M GC, we use a method similar to a finite element approach [24] in which we introduce non-linear equations for each community and for each intercommunity link treating the entire system as a network of communities.We begin with deriving equations for the size of the GC after a percolation process in a single layer and precisely defining the parameters of these equations.We assume that the number of links k ν,i,j linking any node ν in community i to nodes in community j is statistically independent from number k ν,i, , linking node ν to any other community .These numbers are randomly taken from a given degree distributions P i,j (k).We define the generating functions of these distributions: and the generating functions of the excess degree distribution [25]: where k i,j are the average degrees of distributions P i,j (k).We define f i,j as the chance that a link passing from a node in community i to a node in community j does not lead to the GC.The link is an "intra-link" if i is equal to j otherwise it is an "inter-link".If we assume that p j is the fraction of nodes survived in community j after an initial attack or as a result of cascading process, f i,j must satisfy recurrent equations: where the index goes over the set of neighboring communities of community j including community j itself.When = i, one of the links leading from a node in community j back to community i is used by an incoming link, hence we must use the generating function of the excess degree distribution.Finally, the fraction of nodes in community i which belong to the GC is: where the index j goes over the set of neighboring communities of community i including community i itself.If we introduce vectors f with components f i,j , p with components p i and g with components g i , then Eq. ( 3) can be written in a symbolic vector form: This equation can be solved using the iteration method starting with f = 0, and it will uniquely define the vector f ( p) as function of vector p.Analogously, Eq. ( 4) can be presented in a vector form: For generality, we will assume that in Eq. ( 6) the vectors f and p, are two arbitrary vectors, independent of one another.Now we will obtain equations for the M GC of the multiplex.Suppose that the survival probability vector after initial attack is p(0) with components p i (0), and the vector of survival probabilities after stage t of the cascade is p(t).In principle, for the layers of the multiplex A and B the degree distributions can be different, and hence the functions Φ and Ψ and the vectors f and g should be different.Therefore, we will distinguish them by adding indexes A and B. Using the same logic as in Ref. [1], the equations of the cascade of failures starting from t = 0 will be: where g A (t) and g B (t) are the fraction of nodes of each community in the giant component at stage t and p(t) is the effective fraction of survived nodes representing stage t of the cascade of failures as a percolation process after a random attack.As t → ∞ the vectors g A (t) and g B (t) will converge to the mutual giant component P ∞ .
If all distributions P i,j (k) are Poisson distributions as in ER graphs, then G i,j (x) = H i,j (x) = exp[ k i,j (x − 1)] and all probabilities f i,j for the same community j but different i satisfy the same equation and hence they must be equal and we define f j ≡ f i,j .Thus, Eq. (3) and Eq. ( 4) are significantly simplified and become respectively: and From this follows that g When the two layers of multiplex A and B are with average degrees k i,j A , and k i,j B respectively, Eqs.(7) give: (10) which is a generalization of the mutual giant component formula for a multiplex of homogeneous graphs given in Eq. (40) of Ref. [3].We find the numerical solution of the system (10) by iterations starting from the initial values P ∞j = 1.
We next analyze the robustness of our community multiplex model with respect to localized attacks or failures.To this end, we consider the case where all nodes within a radius r h (radius-hole), from the center of the multiplex, are removed from the network.When m is an even number, then the center of the multiplex is at the corner of 4 neighboring ER communities, and else it is in the center of one ER.Note that r h translates into the value of p i (0) by counting the fraction of lattice sites outside the hole of radius r h in the damaged communities.For example, if the center of the hole is in the center of the community i and r h < 1/2, p i (0) ≈ 1−πr 2 h , alternatively, when the center is at the corner of the four communities j, p j (0) ≈ 1 − πr 2 h /4 in each of the four damaged communities.
We find for networks with different system parameters, by simulations and theory, what is the critical radius r c h needed to cause a system collapse.We find the accurate value of r c h through a binary search, where increasing or decreasing of the radius attack is determined by the M GC size at the end of the cascading.At a given radius attack r h , if the M GC size (after the cascading) remains in order of the system size then we increase r h , and otherwise -we decrease it.We define a threshold condi- tion for the M GC size, below which we assume that the M GC is 0. For the numerical calculations of the theory we set the threshold to be 10 −12 , and for the simulations (after some tests) a fraction of 0.1 of the system size seems to provide a good threshold condition.Furthermore, for the numerical calculations we divide each community into ζ 2 of points.Thus, we can calculate numerically with good approximation the fraction of nodes that fail in each community for an attack r h .In addition, since we study the case of a symmetrical two-dimensional square lattice, for the theoretical calculations (using Eq. ( 10)) we set k i,j A and k i,j B to be k inter /4 for i = j.We find that for a network with structure parameters within a certain parameter range of L, ζ and k total , there are two regimes that are divided by a critical α c , see Fig. 2. Hence, a different ratio α between k inter and k total -for a fixed k total -can completely change the system's resilience to localized attacks.For α > α c , we have a metastable regime, where a finite size localized attack larger than r c h causes cascading failures, leading to v system collapse.The critical radius r c h in this regime, for a given k total , depends weakly on the interconnectivity parameter α.Note that this metastable regime located in the narrow interval of k total above k c ≈ 2.4554, where k c is the critical average degree (independent of α) below which the homogeneous ER multiplex collapses without any initial damage [1].In order to interpret our model as a model of a real infrastructure, we must assume that the infrastructure is designed to be as economical as possible: i.e. k total is minimized in a such a way that the multiplex retains the mutual giant component.Moreover, we assume that the initial state of the system before the localized attack is the MGC of the multiplex with given initial k total .Since in the mutual percolation, the final state of the system does not depend on the order in which the damage was made, the final result is the same as if the localized attack was produced simultaneously with the creation of the multiplex as given by Eqs.(7).Remarkably, in this metastable regime networks with the same k total but larger interconnectivity ratio α are more vulnerable to localized attacks than networks with small α where the communities are not well connected, but more self-sufficient.Moreover, in this metastable regime, r c h is independent of the number of communities (Fig. 3).In marked contrast, for α < α c , the critical attack r c h is ∼ 0.5m in ζ units.Namely, the system remains functioning for any attack that is less than removing the entire system.In addition, we obtain numerically based on Eq. ( 10) a phase diagram of r c h , for a large system with m = 100 in Fig. 4 (a).The phase diagram is the same for different m values, see for instance Fig. 4 (b), except for r c h that are in the order of the system size.As we noted above, the initial state of the system must be understood as the MGC obtained for given k total before the localized attack.Therefore, to understand how the damage produced by the localized attack spreads with time, we first produce the cascade of failures with p i (0) = 1.After this cascade stops we produce the localized attack of a given radius r c h .For example, for the simulations in Fig. 5 we perform the attack on step time = 19.After the attack, there is a long latent period during which only a few nodes fail at every time step, and they are located mostly in the vicinity of the attack area.Then, the damage quickly spreads until it reaches the edges of the system.The spreading process explains why the attack size does not depend on the system size.
Discussion.Typically, real-world engineered systems are designed to be as efficient and cost effective as possible.Therefore, such systems are slightly above the critical state, and hence, these systems are very vulnerable to various local failures.Here we have investigated the stability of realistic interdependent networks, consisting of interconnected communities embedded in space, against local failures.We develop a theory for calculating the magnitude of the critical damage needed to destroy the entire system for different parameters of connectivity and spatiality.Our approach is similar to the finite element method which is applied here to the network of communities, where each community is treated as an element, participating in a system of equations.We find that for the same k total (and, hence, the same cost) the networks with low interconnectivity α are more robust against localized attacks than the system in which the communities are well connected.If α is large, the damage produced by the localized attack spreads over the entire system.For small α, the damage does not spread.Thus, the interlinks connecting neighboring communities could serve as vehicles of damage propagation rather than for stabilizing the system.This finding explains why islanding, the strategy employed the electrical engineers by dividing the system into almost isolated self-sustained islands is an efficient strategy against cascading failures in the power grid.In addition, we study the dynamical process of cascading and find a long latent period during which the number of failed nodes is very small and they are localized close to the initial attack.During this period, a relatively small intervention by reinforcing a few nodes can stop the propagation of the cascade of failures.After the latent period is over, the damage quickly spreads over the entire system and there is no economic way to stop it.

FIG. 1 :
FIG. 1:A schematic representation of the model.The nodes are at the lattice sites of a two-dimensional square lattice of size L × L with L = 15.The system is constructed as m × m Erdős-Rényi (ER) networks.Here m = 3, where each ER is of size ζ × ζ with ζ = 5.In each ER network the links are between pairs of sites chosen at random, and there are also links between nearby ER networks.The green and blue lines represent the links in the first and second layer of the multiplex.In our simulations we set periodic boundary conditions, not shown for clarity.

FIG. 3 :
FIG. 3: Analytical results -the dependence of the critical attack size r c h on the number of communities.In (a) and (b) we show three behaviors of r c h for α = 0.4 and α = 0.8.The number of communities is m × m, therefore the dependence on m expresses the dependence of r c h on the number of communities.For relatively small k total , r c h does not depend on the number of communities.For intermediate values of k total , r c h initially grows with m and then from a certain m, reaches a stable value.This stable value oscillates between two values of r c h which correspond to even and odd values of m.For large k total values, r c h grows linearly with m and is approximately 0.5m.In (c) we show (for m = 100) that r c h increases with k total , when for α = 0.4 the transition is sharper than for α = 0.8.This result brightens why the jump in (a) larger than in (b).

FIG. 4 :
FIG. 4: Analytical results -Phase diagram of the critical attack size r c h .Dependence of r c h on the average degree k total and the interconnectivity parameter α.For the phase diagrams in this figure, we sample with equal intervals 16 values for k total and 20 values for α.(a) Contour for the phase diagram with m = 100.The color bar in the right represents the size of r c h in ζ units (in log scale).(b) Here we show the same contour lines as in (a) for r c h values for m = 100 and m = 150.We see that both m values give identical results except for near to the border where r c h ∼ 0.5L (in the last contour line).For lower m values, the misfit of the contour lines starts in lower values of r c h .

FIG. 5 :
FIG. 5: The cascading failures near the critical point.Propagation of a local damage with radius slightly above the critical size r c h .(a) The average distance from the center, r , of the nodes that fail at every iteration.The inset figure is an illustration of the network.(b) The continuous lines represent the size of P∞ j , for 4 communities having different distances d from the center (where the critical hole was removed), as a function of time.The colors of the lines correspond to the colors of the painted communities in the inset figure.The dotted line shows the M GC size over the whole multiplex ( j P∞ j ).For the simulations we set L = 4500, ζ = 300, k total = 2.5 and α = 0.4.The critical size r c h for this simulation, which obtained through a binary search, is r c h = 0.57 in ζ units.