Two distinct transitions in spatially embedded multiplex networks
Michael M. Danziger, Louis M. Shekhtman, Yehiel Berezin, Shlomo Havlin
TTwo distinct transitions in spatially embedded multiplex networks
Michael M. Danziger, Louis M. Shekhtman, Yehiel Berezin, and Shlomo Havlin Department of Physics, Bar Ilan University, Ramat Gan, Israel (Dated: November 6, 2018)Multilayer infrastructure is often interdependent, with nodes in one layer depending on nearbynodes in another layer to function. The links in each layer are often of limited length, due to theconstruction cost of longer links. Here, we model such systems as a multiplex network composed oftwo or more layers, each with links of characteristic geographic length, embedded in 2-dimensionalspace. This is equivalent to a system of interdependent spatially embedded networks in two dimen-sions in which the connectivity links are constrained in length but varied while the length of thedependency links is always zero. We find two distinct percolation transition behaviors depending onthe characteristic length, ζ , of the links. When ζ is longer than a certain critical value, ζ c , abrupt,first-order transitions take place, while for ζ < ζ c the transition is continuous. We show that, thoughin single-layer networks increasing ζ decreases the percolation threshold p c , in multiplex networksit has the opposite effect: increasing p c to a maximum at ζ = ζ c . By providing a more realistictopological model for spatially embedded interdependent and multiplex networks and highlightingits similarities to lattice-based models, we provide a new direction for more detailed future studies. I. INTRODUCTION
Several models have been proposed for spatially em-bedded networks [1–10]. In lattice-based models, linksare only formed to nearest or next nearest neighbors. Inrandom geometric models, links are formed to all neigh-bors within some distance [11, 12]. In models of powergrid topology, links are formed with the m nearest neigh-bors, statically [13] or as a generative model [14]. Somemodels utilize a cost function [15–18] or a characteristicdistance distribution [19, 20] to determine link lengths.The model which we study here has spatiality expressedvia characteristic link lengths. We utilize exponentiallydistributed link lengths, similar to the Waxman model[21].Interdependent networks have been studied mainly onrandom topologies where analytic calculations are pos-sible [24–29]. However, since many real-world com-plex systems are embedded in space, it is importantto understand the properties of interdependent networkswith topology reflecting the dimensionality of the space[16, 30]. This is particularly important when dealing withcritical infrastructure which is heavily influenced by spa-tial constraints [9, 31–33]. An important first step inthat direction is the model presented by Li et al. [34]which models the networks as lattices and includes de-pendency links which are of limited geographic length(described by the parameter r ). This model was shownto include a number of surprising properties includingthree regimes of phase transitions depending on the valueof r . For r < r c , the phase transition is second-orderand appears to be in the same universality class as 2-dimensional lattices [35]. The percolation threshold in-creases as r increases and reaches a maximum value at r c ,where the transition switches to first-order [34, 36]. For r c < r < ∞ the transition is first-order and characterizedby a spreading process, with p c decreasing monotonically.For r = ∞ , the transition is an abrupt simultaneous firstand second order transition [36, 37]. The model was stud- l -4 -3 -2 -1 P ( l ) EU Power GridJapan Local Rail
FIG. 1:
Examples of real-world networks with links ofcharacteristic length.
We examine the distribution of thegeographic lengths of the edges in both the European powergrid [22] (1851 edges) and the inter-station local railway linesin Japan [23] (20745 edges). These networks have links ofcharacteristic length and longer links are exponentially un-likely, as indicated by the linear drop on the semi-logarithmicplot. To compare the two datasets, we rescale the lengths sothat they are measured in units of their own minimum length,which we determine as the peak of the distribution (modelength). The normalization value ( l = 1) corresponds to 3.7km (power) and 1.0 km (rail). The characteristic length asdetermined by the mean is 4.8 km (power) and 1.2 km (rail).As the slope of the exponential fit, it is 3.3 km (power) and2.0 km (rail). The Japan local railway data is formed fromthe complete railway network from [23] with bullet train linesand internal station tracks removed. ied under partial dependency with r = ∞ and it wasshown to have first-order transitions for any fraction ofdependency [38]. When there is partial dependency andfinite dependency link lengths, r c is shown to increaseas the fraction of dependent nodes decreases, divergingat q = 0 [36, 39]. This model was also studied for gen- a r X i v : . [ phy s i c s . s o c - ph ] M a y (a) (b) FIG. 2:
Spatially embedded multiplex networks. (a)
The nodes occupy regular locations in two-dimensional space whilethe links in each layer (blue and green) have lengths that are exponentially distributed with characteristic length ζ = 3 andare connected at random. (b) Cascading failures in multiplex networks. In the first stage the mutual giant connectedcomponent (MGCC) consists of the entire network since all of the nodes are in the giant component of both layers (blue linksand green links). An initial attack on Node 3 causes it and its links to fail. This detaches Node 5 from the giant component ofthe green links, and in the next step it and its links fail. After the failure of Node 5, Node 6 is no longer in the giant componentof the blue links. After Node 6’s failure we find that the remaining nodes are in the giant component of both layers. We notethat the MGCC is not simply the intersection of the giant components in the separate layers. For example, using such anapproach one would conclude that Node 6 is in the MGCC after the failure of Node 3, which is not the case. eral networks formed of interdependent lattices [40] forinterdependent resistor networks with process-based de-pendency [35] and in the presence of healing [41].However, in many real-world systems, the length of thedependency links may not be longer than the length ofthe connectivity links. For instance, in the example ofthe power grid and communications network, it is un-likely that a communications station will skip the nearerpower stations and depend on a power station that is far-ther away. Also, real-world networks including the powergrid do not have uniform link lengths like a lattice, butrather have links of a characteristic length, consistentwith an exponential distribution (see example in Fig. 1)[30]. To address these two issues, here we model interde-pendent networks where every node has a bi-directionaldependency link with the nearest geographic node in theother network. We treat each pair of nodes as singlenodes in a multiplex network, where the links in eachlayer are different but of the same characteristic length ζ . We thus interpret each node as a geographic entity(e.g. a city or neighborhood) which is linked via twotypes of links (e.g. electricity and communications) toother nodes. Each node requires both of its constituentsto function and each constituent requires connectivitywithin its layer.Our main focus in this paper is to examine the roleof the characteristic length ( ζ ) of the links on the ro- bustness of the multiplex. We find that, though increas-ing ζ decreases p c in single networks–making them morerobust–it has the opposite effect on multiplex networks.Increasing ζ increases p c for multiplex networks until acritical length ζ c where p c is maximal. At ζ c , the per-colation transition changes to first-order and the multi-plex network is susceptible to spreading cascading fail-ures similar to the ones observed in interdependent lat-tices [34, 36, 39, 42].By demonstrating that comparable critical behavioremerges in more realistic topologies, we show that thecritical behavior demonstrated in previous lattice-basedmodels is not limited to the specific implementation orthe lattice topology but is rather a generic property ofinterdependent spatially embedded networks. Further-more, by providing a topological model that more closelymatches real-world systems, we provide a more realistictopological framework for future studies of interdepen-dent critical infrastructure. II. MODEL
To model connectivity links of characteristic length ineach layer, we construct the network as follows. We be-gin by assigning each node an ( x, y ) coordinate with in-tegers x, y ∈ [0 , , . . . , L ). To construct the links in each p P ∞ ζ =0 . ζ =2 ζ =6 ζ =15 ζ =1000 (a) single-layer p P ∞ ζ =0 . ζ =2 ζ =6 ζ =15 ζ =1000 (b) multiplex FIG. 3:
Percolation of (a) single-layer and (b) multiplex networks with links of characteristic geographiclength.
The fraction of the network in the largest connected component ( P ∞ ) as a function of the fraction of nodes remainingafter random removal of a fraction 1 − p nodes. (a) In single networks, the transition is always second-order and the criticalbehavior above p c is the same for all finite values of ζ ( ζ (cid:28) L ). The value of p c decreases quickly and monotonically with ζ ( ζ > (b) In multiplex networks, the transition is comparable to single networks for ζ = 0 . p c ≈ . p c increases as ζ increases. This continues until the maximal value is reached (at ζ c ) and the transition becomes first-order. The case of ζ = 1000 is very close to purely random and has p c ≈ . / (cid:104) k (cid:105) = 0 . (cid:104) k (cid:105) = 4 and L = 4000. layer, we select a source node at random with coordinates( x , y ) and draw a length l according to P ( l ) ∼ e − l/ζ .We choose the permitted link length ( dx, dy ) which isclosest to fulfilling l = (cid:112) dx + dy , select one of theeight length-preserving permutations ( dx ↔ − dx, dy ↔− dy, dx ↔ dy ) uniformly at random and make a link tonode ( x , y ) with x = x + dx , y = y + dy . Thisprocess is executed independently in each layer and iscontinued until the desired number of links ( N (cid:104) k (cid:105) /
2) isobtained. For simplicity, we use the same characteristiclength ζ and average degree (cid:104) k (cid:105) in each layer. However,because they are constructed independently, the links ineach layer are different (as demonstrated in Fig. 2a andFig. 6) and this disorder enables the critical behaviorwhich we describe below.We then perform site percolation by removing a frac-tion 1 − p of the nodes from the system and finding themutual giant component [24]. When a node is removed, itcauses damage to nodes in both layers due to the connec-tivity links in each layer, which are also removed. How-ever, since the connectivity links are not the same in bothlayers, there will be nodes that are connected to the giantcomponent in one layer but are disconnected in the otherlayer. Since the node functionality requires connectivityin both layers, such nodes will fail, causing further dam-age in the system. This leads to the cascading failures asdemonstrated in Figs. 2b and Fig. 5, which are similarto those described in [24, 34, 36, 37, 42–44]. III. RESULTS
Since we are not aware of a discussion of the perco-lation properties of this topology for single networks,we briefly describe those properties here and in the Ap-pendix. In the limit of ζ →
0, the only permitted linkswill be to nearest neighbors (because links of length < (cid:104) k (cid:105) = 4, we recover the standard2-dimensional percolation behavior with p c ≈ . ζ increasesfurther, the robustness increases and in the limit ζ → ∞ ,all lengths are equally likely to be drawn and the topol-ogy reverts to purely random (Erd˝os-R´enyi topology)with p c = 1 / (cid:104) k (cid:105) = 0 .
25 (see Fig. 3a.) Thus, similar torewiring probability in the original small world model [4],we have a single parameter ζ which allows us to smoothlytransition from lattice to random topology. For all valuesof ζ , a single network undergoes a second-order transition(Fig. 4c and 4d).In multiplex networks, where connectivity to the giantcomponent in both layers is required, cascading failuresemerge [24, 26, 47, 48]. For ζ ≈
0, large cascading failuresdo not emerge. This is because the multi-layer structureis mostly redundant and the difference between connec-tivity in one layer or both layers is negligible. However,once ζ becomes long enough (above ζ c ), intensive cascad-ing failures emerge and the system undergoes an abrupt,first-order transition, similar to the transition in inter-dependent lattices [34], as shown in Fig. 4 and 5. As ζ
10 20 30 40 50 ζ p c (a) (cid:104) k (cid:105) = 3
10 20 30 40 50 ζ p c (b) (cid:104) k (cid:105) = 4
10 20 30 40 50 ζ P ∞ ( p c ) multisingle (c) (cid:104) k (cid:105) = 3
10 20 30 40 50 ζ P ∞ ( p c ) multisingle (d) (cid:104) k (cid:105) = 4 FIG. 4:
The effect of the characteristic link length ζ on the percolation threshold and size of the giantcomponent at criticality. (a-b) The percolation threshold in multiplex networks increases until it reaches a peak at ζ c andthen decreases slightly to its asymptotic value of 2 . / (cid:104) k (cid:105) . (c-d) The size of the order parameter P ∞ at p c . Single-layernetworks always have an order parameter close to zero at p c , which indicates a second-order transition. Multiplex networkshave a second-order transition for low values of ζ but at ζ c this jumps to a large fraction of the system size, indicating afirst-order transition. All panels have L = 4000 and are plotted based on at least 10 realizations of each system. The averagesand raw data are plotted for all systems. becomes even longer, p c decreases and slowly approachesits asymptotic value of 2 . / (cid:104) k (cid:105) as known from inter-dependent Erd˝os-R´enyi networks [24], (Fig. 4a and 4b).In interdependent lattices the mechanism of the first-order transition for dependency links of large finite lengthis a propagating spinodal interface [34, 36, 39, 42]. Ona microscopic level, the first-order transition takes placeby the emergence of a hole due to random fluctuationsof characteristic length ξ (the percolation correlationlength) which then propagates through the system. Thispropagation is enabled by the cascade dynamics (cf. Fig.2b) and dependency links which relay the damage causedby the hole into a concentrated area around the hole’sedge. It would seem that in order for this phenomenon totake place, long dependency links are needed and longerconnectivity link lengths would be insufficient to sustaindamage propagation. The reason that the dependencylinks have a stronger influence on robustness is that inorder for a node to fail due to connectivity, all of its linksto the giant component need to fail whereas with depen-dency, an otherwise well-connected node will fail if thesingle node that it depends on fails.We find that, in fact, this is not the case. Connectiv-ity link lengths which are above a critical length, ζ c , butmuch smaller than the total system size are sufficient tocause the cascading failures observed in lattice models,even when the dependency links have zero length. In-deed, the first order transition lacks scaling behavior justabove p c (Fig. 3b) and is characterized by a slow spread-ing process (Fig. 5), just like interdependent lattices de-scribed in Refs. [34, 36, 39]. However, the maximum p c in this system is lower than the comparable lattice sys-tem: ≈ .
64 ( (cid:104) k (cid:105) = 4 here) vs. ≈ .
74 (interdependentlattices). This is because, in order for the connectivitylinks to effectively relay the damage from the hole, thesystem must be closer to criticality. Indeed, only when the average degree of each node is 1 do the connectivityand dependency links have the same effect upon it.Surprisingly, the highly localized topology whichemerges from spatial embedding makes the system moresusceptible to cascading failures. This is due to the factthat, because the damage from an emergent (or induced[42]) hole is relayed by the cascade dynamics to the neigh-borhood of its interface, the nodes near the edge of thehole are far more likely to become disconnected. Thisis related to work on information diffusion in social net-works, where it was shown that high modularity makesviral cascades more likely to occur due to the increasedlikelihood of multiple exposure to the information [49–51].In single layer networks, p c decreases monotonicallyas ζ increases from ζ ≈ . ζ = ∞ . Incontrast, in multiplex networks, p c increases until ζ c andthen decreases monotonically thereafter. The peak in p c ( ζ ) is due to the fact that the size of the critical holethat is needed to trigger the transition scales linearly with ζ [42] and the size of emergent holes above the percola-tion threshold, ξ ( p ), decreases with p [45]. This wouldindicate that the smaller ζ is, the smaller the criticalhole needs to be and that p c would increase monotoni-cally as ζ decreases. However, when ζ < ζ c there is notenough space between the emergent hole and the extentof the damage propagation ( ζ ) for the network to dis-integrate and the small emergent holes remain in place[34, 36, 39, 42].For interdependent lattices with dependency links offinite length, the single network case is recovered whenthe dependency links have length zero ( r = 0) [34]. Inspatially embedded multiplex networks the same limitis recovered as ζ → (a) t = 30 (b) t = 100(c) t = 200 (d) t = 300 FIG. 5:
Dynamic evolution of cascading failures.
Herewe show the nodes of the multiplex network, colored blackif functional and white if not. After the emergence of alarge enough hole (due to random fluctuations), the dam-age is relayed predominantly to the area around the edge ofthe hole (due to the finite link lengths) which leads to thenodes around the interface becoming disconnected and prop-agation of the damage until the system disintegrates. Similardynamics have been observed in interdependent lattices withdependency lengths of finite length [34, 36]. In this figure L = 4000, ζ = 15, (cid:104) k (cid:105) = 4. plex network is called intersimilarity [43, 44] or overlap[52, 53]. The cascading failures and abrupt transitionswhich characterize interdependent networks decrease asoverlap increases. In the limit of total overlap, they dis-appear altogether because the system is composed of twoexact copies of the same network and the cascade shownin Fig. 2b does not take place at all.In multiplex networks with links of characteristiclength, the extent of the overlap can be estimating byconsidering the probability that, given the same sourcenode, two links lead to the same target node. In the con-tinuum limit this is proportionate to the probability thatthe links have the same length and the same direction.The system is isotropic by construction so the directionalcondition is simply 1 / πr , the size of a ring of radius r .We obtain P ( overlap ) ∼ (cid:90) ∞ P ( r )2 πr dr ∼ ζ (cid:90) ∞ e − r/ζ πr dr ∼ ζ (1)We find that the scaling in our system is scale-free, withan exponent of ≈ − . -1 ζ -7 -6 -5 -4 -3 -2 -1 o v e r l a p f r a c t i o n L =500 L =1000 L =5000 L =10000 FIG. 6:
The fraction of overlapping nodes.
The fractionof overlapping nodes is determined as the number of commonlinks across both layers divided by the total number of linksin each layer. When ζ ≈
0, the overlap is maximal and thenetworks are identical (for (cid:104) k (cid:105) = 4, as in this figure). As ζ increases, the fraction decreases with an exponent of ≈ − . the deviation from the continuum calculation is due tothe fact that with the discretization of space that weintroduce, links that would otherwise have been distinctare unified and the critical exponent is reduced from − ≈ − . IV. DISCUSSION
Previous models of spatially embedded interdependentnetworks [34, 36, 38–40, 42] have used two-dimensionallattices with dependency links connecting nodes fromone network to the other. The dependency links werealso affected by the spatial embedding via the restrictionthat they have length of up to r , a system parameter.This model led to many important results, but left sev-eral important issues unaddressed. First, the topology ofreal-world spatially embedded networks is not lattice-likeor even strictly planar and it was not clear that resultsderived on lattices would accurately describe real-worldtopologies. Second, the assumption that dependencylinks are longer than connectivity links does not corre-spond with what we would expect from critical infras-tructure: It is not reasonable to expect a communicationsstation to get power from a distant power plant and notthe one nearest to it. Here we address these problems bymodeling spatially embedded interdependent networks as
10 20 30 40 50 ζ p c (a) (cid:104) k (cid:105) = 3
10 20 30 40 50 ζ p c (b) (cid:104) k (cid:105) = 4 FIG. 7:
Dependence of p c on ζ for single networks Thepercolation threshold p c drops quickly as a function of ζ andby ζ ≈
10 it is already very close to p c = 1 / (cid:104) k (cid:105) , the valuefrom Erd˝os-R´enyi networks. multiplex networks where the dependency relationship isto the nearest node in the other layer and the connectiv-ity links are of finite characteristic length but not uniformor regular.We find that the most important features of the lattice-based models are reproduced by our new model: second-order and first-order transitions depending on the linklength, a first-order transition characterized by the emer-gence and spreading of a hole and substantially highervulnerability compared to single-layer networks. Fur-thermore, p c has a maximal value at the ζ value ( ζ c ) thatseparates the two types of transitions. This validates pre-vious work based on lattices and also shows a new wayforward for the modelling of critical infrastructure andother spatially embedded multilayer networks. Acknowledgments
We acknowledge the MULTIPLEX (No. 317532) EUproject, the Deutsche Forschungsgemeinschaft (DFG),the Israel Science Foundation, ONR and DTRA for fi-nancial support. We also thank Sergey V. Buldyrev forhelpful discussions and comments on the model.
Appendix: Percolation threshold in single networksand in multiplex and single networks for ζ < In single networks, for ζ > p c decreases monotoni-cally with increasing ζ . This is evident in Fig. 7, where p c is shown to rapidly approach 1 / (cid:104) k (cid:105) , the Erd˝os-R´enyivalue, as ζ is increased [54, 55]. Thus, we find that theeffect of increasing ζ on single networks is to make them more robust, the opposite of its effect in multiplex net-works. The behavior of single and multiplex spatiallyembedded networks is different for 0 < ζ <
3. Since, thisis not the main focus of this study and since it does notaffect the critical phenomena, we have avoided discussing ζ p c multisingle (a) (cid:104) k (cid:105) = 3 ζ p c multisingle (b) (cid:104) k (cid:105) = 4 FIG. 8:
Percolation behavior for very short character-istic lengths.
The percolation threshold increases slightlyfor all systems with very small ζ . For single networks itrapidly falls again but for multiplex networks, p c remainshigher for longer before falling to its ζ = 0 level. This is be-cause in single networks, the only effect is the relative weak-ness of lattices with complex neighborhoods. In multiplexnetworks, this effect is compounded with the weakness en-gendered by clustering [56, 57] and overlap [43, 44, 52, 53]. it in the main text. When ζ = 0, only the shortest pos-sible links will be drawn (those with l = 1). Thus, thetopology of the network layer will be a pure lattice for (cid:104) k (cid:105) = 4 and a diluted lattice for (cid:104) k (cid:105) <
4. However, for0 < ζ <
3, some of the nearest neighbor links are ex-changed for next nearest neighbor (diagonal) links whichare less robust than nearest neighbor only links. Thisis a well-known result from bond percolation on latticeswith complex neighborhoods, where it was found that (cid:104) k c (cid:105) for a next-nearest neighbor lattice ( z = 8) is higherthan (cid:104) k c (cid:105) for a nearest neighbor lattice ( z = 4) [58]. Forinstance, for ζ ≈ . ζ increases, the links are lessredundant and tend to strengthen the network, as seenin Fig. 7 and Fig. 8 for single networks.The robustness is more severely impacted in the mul-tiplex case. This is due to the effects of clustering anddisorder. Combining nearest and next-nearest neighborlinks leads to high clustering (number of triangles) which,though having only a minor impact on the robustnessof single networks, substantially weakens interdependentnetworks [56, 57]. Furthermore, the increased disorderdecreases the fraction of overlapping links substantially( ≈
50% as in Fig. 6). These effects combine in themultiplex case to cause a substantial increase in p c for0 . < ζ <
1. At ζ ≈
1, the links are no longer redun-dant, clustering decreases, and the system becomes morestable again. However, it is important to note that eventhough p c is increased in this interval, the transition isstill second-order. The first-order transition only takesplace when the links are substantially longer ( ζ > ζ c ),and the spreading process described above can take place. [1] M. Doar and I. Leslie, in INFOCOM’93. Proceedings.Twelfth Annual Joint Conference of the IEEE Computerand Communications Societies. Networking: Foundationfor the Future, IEEE (IEEE, 1993), pp. 82–89.[2] L. Wei and D. Estrin, Submitted to INFOCOM (1993).[3] E. W. Zegura, K. L. Calvert, and M. J. Donahoo,IEEE/ACM Transactions on Networking (TON) , 770(1997).[4] D. J. Watts and S. H. Strogatz, Nature , 440 (1998).[5] M. Penrose, Random geometric graphs , vol. 5 (OxfordUniversity Press Oxford, 2003).[6] J. Kleinberg, in
Proceedings of the thirty-second annualACM symposium on Theory of computing (ACM, 2000),pp. 163–170.[7] K. Kosmidis, S. Havlin, and A. Bunde, EPL (EurophysicsLetters) , 48005 (2008).[8] D. Li, G. Li, K. Kosmidis, H. E. Stanley, A. Bunde, andS. Havlin, EPL (Europhysics Letters) , 68004 (2011).[9] M. Barth´el´emy, Physics Reports , 1 (2011).[10] T. C. McAndrew, C. M. Danforth, and J. P. Bagrow,arXiv preprint arXiv:1501.05976 (2015).[11] A. F. Rozenfeld, R. Cohen, D. ben Avraham, andS. Havlin, Phys. Rev. Lett. , 218701 (2002).[12] M. Bradonji´c, A. Hagberg, and A. G. Percus, in Algo-rithms and Models for the Web-Graph (Springer, 2007),pp. 209–216.[13] P. Hines, S. Blumsack, E. Cotilla Sanchez, and C. Bar-rows, in
System Sciences (HICSS), 2010 43rd Hawaii In-ternational Conference on (IEEE, 2010), pp. 1–10.[14] D. Deka and S. Vishwanath, in (2013), pp. 591–598, ISBN 978-1-4799-3211-5/13.[15] S. S. Manna and A. Kabakioglu, Journal of Physics A:Mathematical and General , L279 (2003).[16] M. T. Gastner and M. E. Newman, The European Phys-ical Journal B-Condensed Matter and Complex Systems , 247 (2006).[17] T. Emmerich, A. Bunde, and S. Havlin, Physical ReviewE , 062806 (2014).[18] Y. Ren, M. Ercsey-Ravasz, P. Wang, M. C. Gonz´alez, andZ. Toroczkai, Nature Communications , 5347 (2014),ISSN 2041-1723.[19] Z. Wang, R. J. Thomas, and A. Scaglione, GeneratingRandom Topology Power Grids (Institute of Electricaland Electronics Engineers, 2008), pp. 183–183, ISBN 0-7695-3075-8.[20] P. Grassberger, Journal of Statistical Mechanics: Theoryand Experiment , P04004 (2013).[21] B. M. Waxman, Selected Areas in Communications,IEEE Journal on , 1617 (1988).[22] Q. Zhou and J. Bialek, Power Systems, IEEE Transac-tions on , 782 (2005), ISSN 0885-8950.[23] National Land Information Division, National SpatialPlanning and Regional Policy Bureau, MILT of Japan, National railway data , http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02.html (2012).[24] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, andS. Havlin, Nature , 1025 (2010).[25] J. Gao, S. V. Buldyrev, H. E. Stanley, and S. Havlin, Nature Physics , 40 (2012).[26] T. P. Peixoto and S. Bornholdt, Phys. Rev. Lett. ,118703 (2012).[27] F. Radicchi and A. Arenas, Nature Physics , 717 (2013).[28] M. Kivel¨a, A. Arenas, M. Barth´el´emy, J. P. Gleeson,Y. Moreno, and M. A. Porter, Journal of Complex Net-works , 203 (2014).[29] S. Boccaletti, G. Bianconi, R. Criado, C. Del Genio,J. G´omez-Garde˜nes, M. Romance, I. Sendina-Nadal,Z. Wang, and M. Zanin, Physics Reports (2014), ISSN0370-1573.[30] D. Li, K. Kosmidis, A. Bunde, and S. Havlin, NaturePhysics , 481 (2011).[31] S. Rinaldi, J. Peerenboom, and T. Kelly, Control Sys-tems, IEEE , 11 (2001), ISSN 1066-033X.[32] P. Hokstad, I. Utne, and J. Vatn, Risk and Interdepen-dencies in Critical Infrastructures: A Guideline for Anal-ysis , Springer Series in Reliability Engineering (Springer,2012), ISBN 9781447146612.[33] D. Helbing, Nature , 5159 (2013), ISSN 1476-4687.[34] W. Li, A. Bashan, S. V. Buldyrev, H. E. Stanley, andS. Havlin, Phys. Rev. Lett. , 228702 (2012).[35] M. M. Danziger, A. Bashan, and S. Havlin, New Journalof Physics , 043046 (2015).[36] M. M. Danziger, A. Bashan, Y. Berezin, and S. Havlin,Journal of Complex Networks , 460 (2014).[37] D. Zhou, A. Bashan, R. Cohen, Y. Berezin, N. Shnerb,and S. Havlin, Phys. Rev. E , 012803 (2014).[38] A. Bashan, Y. Berezin, S. V. Buldyrev, and S. Havlin,Nature Physics , 667 (2013).[39] M. M. Danziger, A. Bashan, Y. Berezin, and S. Havlin,in Signal-Image Technology Internet-Based Systems(SITIS), 2013 International Conference on (2013), pp.619–625.[40] L. M. Shekhtman, Y. Berezin, M. M. Danziger, andS. Havlin, Phys. Rev. E , 012809 (2014).[41] M. Stippinger and J. Kert´esz, Physica A: Statistical Me-chanics and its Applications , 481 (2014), ISSN 0378-4371.[42] Y. Berezin, A. Bashan, M. M. Danziger, D. Li, andS. Havlin, Scientific Reports (2015).[43] R. Parshani, C. Rozenblat, D. Ietri, C. Ducruet, andS. Havlin, EPL (Europhysics Letters) , 68002 (2010).[44] Y. Hu, D. Zhou, R. Zhang, Z. Han, C. Rozenblat, andS. Havlin, Phys. Rev. E , 052805 (2013).[45] A. Bunde and S. Havlin, Fractals and disordered systems (Springer-Verlag New York, Inc., 1991).[46] R. Ziff, Phys. Rev. Lett. , 2670 (1992).[47] G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, andJ. F. F. Mendes, Phys. Rev. Lett. , 248701 (2012).[48] M. M. Danziger, A. Bashan, Y. Berezin, L. M. Shekht-man, and S. Havlin, in Nonlinear Dynamics of ElectronicSystems , edited by V. Mladenov and P. Ivanov (SpringerInternational Publishing, 2014), vol. 438 of
Communica-tions in Computer and Information Science , pp. 189–202,ISBN 978-3-319-08671-2.[49] D. Centola, Science , 1194 (2010).[50] L. Weng, F. Menczer, and Y.-Y. Ahn, Scientific Reports , 2522 (2013), ISSN 2045-2322.[51] A. Nematzadeh, E. Ferrara, A. Flammini, and Y.-Y.Ahn, Phys. Rev. Lett. , 088701 (2014). [52] D. Cellai, E. L´opez, J. Zhou, J. P. Gleeson, and G. Bian-coni, Phys. Rev. E , 052811 (2013).[53] M. Li, R.-R. Liu, C.-X. Jia, and B.-H. Wang, New Jour-nal of Physics , 093013 (2013).[54] M. Newman, Networks: an introduction (OUP Oxford,2010).[55] R. Cohen and S. Havlin,
Complex Networks: Structure,Robustness and Function (Cambridge University Press,2010), ISBN 9781139489270. [56] Huang, Xuqing, Shao, Shuai, Wang, Huijuan, Buldyrev,Sergey V., Eugene Stanley, H., and Havlin, Shlomo, EPL , 18002 (2013).[57] S. Shao, X. Huang, H. E. Stanley, and S. Havlin, Phys.Rev. E , 032812 (2014).[58] X. Feng, Y. Deng, and H. W. J. Bl¨ote, Phys. Rev. E78