Abrupt transition due to non-local cascade propagation in multiplex systems
AAbrupt transition due to non-local cascade propagation in multiplex systems
Oriol Artime and Manlio De Domenico
CoMuNe Lab, Fondazione Bruno Kessler, Via Sommarive 18, 38123 Povo (TN), Italy (Dated: September 21, 2020)Multilayer systems are coupled networks characterized by different contexts (layers) of interactionand have gained much attention recently due to their suitability to describe a broad spectrum of em-pirical complex systems. They are very fragile to percolation and first-neighbor failure propagation,but little is known about how they respond to non-local disruptions, as it occurs in failures inducedby flow redistribution, for example. Acknowledging that many socio-technical and biological systemssustain a flow of some physical quantity, such as energy or information, across the their components,it becomes crucial to understand when the flow redistribution can cause global cascades of failuresin order to design robust systems,to increase their resilience or to learn how to efficiently dismantlethem. In this paper we study the impact that different multiplex topological features have on therobustness of the system when subjected to non-local cascade propagation. We first numericallydemonstrate that this dynamics has a critical value at which a small initial perturbation effectivelydismantles the entire network, and that the transition appears abruptly. Then we identify that theexcess of flow caused by a failure is, in general, more homogeneously distributed the networks inwhich the average distance between nodes is small.Using this information we find that aggregatedversions of multiplex networks tend to overestimate robustness, even though to make the systemmore robust can be achieved by increasing the number of layers. Our predictions are confirmed bysimulated cascading failures in areal multilayer system.
I. INTRODUCTION
Power grids, trophic relations, online micro-blogging services or scientific collaborations: complex networks areubiquitous in nature [1, 2]. Since the nascence of modern network science more than two decades ago, efforts havebeen made towards finding common characteristics among complex systems of different sort. One of the most salientexamples is the robust-yet-fragile property [3]: it suffices to have a large variance in the distribution of the number ofconnections to observe that networks are very good at handling random failures but are extremely weak when facinglocalized attacks [4, 5]. The study of network robustness sheds light on the characteristics that play a central roleon maintaining the proper functioning of these networks and how their different properties are altered when somecomponents are corrupted. This has been a topic of paramount importance in the research agenda of network sciencedue to its conceptual, modelling challenge [6–14] and its wide range of applications [15–21].From the modeling side, the simplest conceptual way to consider the problem of network robustness is to disregardany intrinsic dynamics in the nodes or in the links, and look for the geometrical properties of the remaining networkafter the failure or the attack from the viewpoint of classical percolation [22], such as the size of the giant connectedcomponent, the cluster size distribution and the nature of the phase transition. The next level of complexity includesthe spreading of failures across the network, for example, by letting a node collapse if its fraction of failed neighbors ishigher than a threshold [23], by mapping the propagation to a branching process [24, 25], or even by combining differentrules [26]. All these models share the limiting feature that the failure propagates via first-neighbors connections. Thisallows mathematical tractability but does not capture the non-local propagation observed in some real-world cascadefailures, e.g. in the 1996 disturbance of the Western Systems Coordinating Council (WSCC) system [27], in the2003 blackout in northeast USA [28] or in the air-traffic disruption due to the eruption of the Icelandic volcanoEyjafjallaj¨okull [29, 30]. The family of models that offers a plausible explanation for this non-local chain of failuresis the one that considers the flow of some quantity across the network [9, 10, 31–34], such as the electrical energy in a r X i v : . [ phy s i c s . s o c - ph ] S e p the power grid or the number of passengers between public transport stations. When the structure of the network isperturbed, for instance by removing a small subset of nodes, an overall flow redistribution is observed and subsequentfailures might occur, not necessarily close to the previous failure.The flow redistribution can occur in various ways depending on the characteristics of each system. Embracingthis heterogeneity, several types of models have been proposed in the past. As it is usual in mechanistic models ofstatistical physics, we have simple, toy models that seek plausible physical mechanisms that are able to explain themost remarkable features of the empirical observations, providing qualitative explanations of emergent behaviors viawell-isolated cause-effects relations [35–37]. One of the most paradigmatic examples of this family of models is the oneproposed by Motter and Lai [9], where nodes have a binary state — operative or failed — and the flow is completelydetermined by the topology. Since our goal is not the accurate reproduction of behavior of a particular system but thequalitative understanding of the role played by several multiplex properties in the propagation of non-local cascades,throughout the article we adopt the aforementioned framework. Further details of our model are given in Section II.It is worth mentioning, however, that there is a plethora of models proposed to describe the phenomenon of non-localcascades, with different level of abstraction and whose suitability depends on the particular scientific question onewants to address, see for example [38–42]. The general tendency is that the more details introduced in a model, thebetter it can reproduce empirical observations but at the expense of less transparent interpretations of the role playedby each of its parameters.The topological substrates on which we consider the non-local cascade propagation driven by load redistribution aremultiplex networks [43, 44]. These are graphs composed by a same set of nodes interacting through different contexts,the layers. The multilayer framework is particularly useful to model systems with mutual interdependencies, a real-world feature, and they display a very rich phenomenology that cannot be otherwise observed in aggregated, single-layer networks [43, 45–48]. Regarding the issue of the robustness in multilayer networks, the same division mentionedabove is valid: percolation analyses [49–51] and first-neighbor cascade propagation models are abundant [16, 52–55],while non-local cascade propagation is a promising, but still largely unexplored topic, except for a few exceptions [56].Similarly, the characterization of non-local failures in real multilayer systems is certainly interesting to understandhow failures can propagate within and across layers. This is indeed possible: for instance, in the examples of realsystems displaying non-local cascades given above (power grid and European Air Transportation network) accepta description as multilayer structures [57, 58]. To the best of our knowledge this has not been done yet, althoughit would be a very nice contribution because it would provide a benchmark to compare with the outcome of thedynamical models.We work under the assumption that connectivity is a first proxy for functionality, therefore the robustness ofthe network is related to the size of the largest connected component once the cascade stops. We numericallyshow that overload cascades can display an abrupt phase transition [59] that separates the regimes in which theinitial perturbation completely dismantles the network or leaves it almost intact. The transition is characterized bytechniques borrowed from the theory of critical phenomena. Beyond the order of the phase transition, which certainlyhas implications for network fragility assessments, we also address the question of how the robustness is modified ifthe multiplexity of a network is disregarded and, instead, the aggregated system is considered. In the same vein, weexplain how the number of layers impact the overall robustness. To answer these questions, we present compellingevidence that when the perturbed flow can be more homogeneously redistributed across the network, the tendency ofrobustness is to increase, and quantify this effect by means of the average path length [2]. To check the validity of ourresults, we apply our model to the European air transportation multiplex network, finding a good qualitative match. II. MODEL DEFINITION
The model that we use is inspired by the load-capacity model of Ref. [9]. In our case we deal with a multiplextopology formed by M independent networks ( M ≥ N nodes and created from the same degree distribution but with different degree sequences. Besides the intralayerlinks, every node is one-to-one connected with its replica node in the M − i in layer m is assigned a load L mi and a capacity C mi . Loads are computed as the number of shortest paths crossing each node(the interlayer connections also count in the computation) and vary as the cascade propagates. The capacity is afixed quantity, proportional to the initial load L m,inii , that is, C mi = (1 + α ) L m,inii . (1)The parameter α is called the tolerance, and gives an idea of the ability of the nodes to readapt to flow changes.When the load of a node exceeds its capacity, L mi > C mi , that node automatically fails by breaking all its connections.Moreover, the replicas of the overloaded node in the other layers also fail by disconnecting themselves from theirneighbors. We consider cascades that propagate in a simultaneous manner, that is, at each time step all overloadednodes fail at once, and then the loads are recalculated. To simplify the analysis, the initial perturbation is performedby deleting the node with highest load, i.e., we are dealing with intentional attacks that assume global knowledge ofthe network. Despite being single-node attacks, we will see that the consequences can be catastrophic for the integrityof the network.The load-capacity model by Motter and Lai was originally introduced in single-layer networks and it providedinsights on how the initial load distribution and thus, the topology, can influence the robustness of the system.However, the nature of the transition — continuous or discontinuous — is a problem that was left unexplored in theoriginal article and has been ignored in other works that followed up. To the best of our knowledge, the only attemptto characterize it is Ref. [60], that finds a first-order transition but for massive initial random attacks, instead of single-node attacks. Here our focus is on the response of multiplex systems to the overload cascades, but, additionally, weshed light on the overlooked problem of the order of the phase transition in monoplexes subjected to single-nodetargeted attacks. In both cases we find an abrupt transition, evincing that it is the dynamics and not the topologywho determines the nature of the phase transition.So far we have been using the concepts of robustness and fragility in an intuitive way. To be more precise, we employthem following the next criterion. On the one hand, a network experiencing a first-order transition is considered bydefault more fragile than one experiencing a continuous transition, because the transition is more abrupt, hard toanticipate and of more destructive power. On the other hand, when dealing with the same type of transition, wesay that the lower the critical point (critical tolerance α c ), the more robust the network is. When critical points aresimilar, then the faster the order parameter saturate to its maximum value, the more robust the network is considered.That is because the maximum value corresponds to the state of no damage. III. RESULTS
We start our analysis by studying duplexes of random regular (RR) graphs and Erd˝os-R´enyi (ER) networks, whichare the simplest non-trivial network topologies. The degree distribution of the former is a Dirac delta P ( k ) = δ ( k − k c ),i.e., all nodes share the same number of connections k c , while in the latter case P ( k ) is well-approximated by a Poissondistribution of parameter (cid:104) k (cid:105) . In Fig. 1 we show the effects of the cascade propagation in these types of graphs. Wetake as order parameter the average of the size of the largest connected component (LCC) at the end of the cascadepropagation, (cid:104) S (cid:105) . We see that, for both topologies, the destructive effects of the initial perturbation are barely . . . . Tolerance α . . . h S i ( a ) RR . . . L ini ( × )036 P ( L i n i ) × − . . . . Tolerance α . . . h S i ( b ) E R . . . L ini ( × )10 − − − P ( L i n i ) FIG. 1. Size of the largest connected component once the cascade propagation stops, for duplexes of random regular graphs ( a )and Erd˝os-R´enyi networks ( b ). In ( a ) all nodes have the same degree k c = 4, while in ( b ) the mean degree is (cid:104) k (cid:105) = 4 and aminimum degree k = 2 has been imposed to guarantee global connectivity. In both cases the number of nodes per layer is N = 5000. Error bars correspond to 1 standard deviation, although in most of the points the error it is smaller than the markersize. Averages are taken over 1000 realizations of the dynamics. In the insets, the probability density of the initial load L ini computed from 100 independent networks. appreciated for almost any tolerance α . The only possible disruption occurs when the network is initially operatingvery close to its maximum capacity, i.e., when α →
0. The reason behind this high robustness is that nodes arestatistically equivalent among them, so they all have an initial load of the same order of magnitude — see the insetsof Fig. 1 which show that the initial load distributions are bounded and have well-defined mean and variance. Whena node is attacked, be it the one with the largest load or one chosen at random, the excess of load that needs to behandled by other nodes is a small fraction of their own load. Recalling that the failure condition depends on individualinitial loads, it suffices a small tolerance to avoid the cascade propagation.We next move our focus to cascading dynamics in a duplex of scale-free networks. This system displays a veryheterogeneous load distribution [61], hence we expect a different qualitative behavior in its robustness properties.Indeed, it is immediate to see [Fig. 2( a )] that the transition from a completely destroyed network to an almostintact one occurs at a finite tolerance. This critical value α c separates two actual regimes, or phases in the parlance ofcritical phenomena: one in which the initial perturbation propagates at a network-wide scale, dismantling the networkfunctionality, and another in which the perturbation barely affects a few nodes. Therefore scale-free multiplex networksare weaker than RR and ER multiplexes, because, when an attack occurs, the nodes need more resources (tolerance)to absorb the redistributed flow and avoid the global failure.Knowing whether the nature phase transition is continuous or discontinuous is useful not only to shed light on thedynamics of cascade propagation but also for robustness assessments. Indeed, a discontinuous transition is much moresudden, difficult to anticipate and devastating than a continuous one. The order parameter in Fig. 1 jumps abruptlyfrom 0 to 1, thus we cannot rule out the possibility that the phase transition in multiplex scale-free networks is alsodiscontinuous, even though the growth of (cid:104) S (cid:105) is much smoother in this case. The first thing to note is that the errorbars of Fig. 2 in the critical region are surprisingly large. Their origin can be understood by plotting the probabilitydensity function of the size of the LCC, P ( S ) [see the inset of Fig. 2( a )]. For low values of the tolerance, P ( S ) issingle-peaked very close to zero. By increasing α we enter in the critical region and another peak develops, away fromthe first one. The position of these peaks is barely dependent on α , but their relative height is. As the tolerance keepsincreasing, the second peak becomes higher, until eventually the first peak disappears. This bimodality is a typical . . . . Tolerance α . . . h S i ( a ) S P ( S ) .
075 0 .
150 0 . Tolerance α − − B i nd e r U N ( b ) /M N m i n ( α ) FIG. 2. Characterization of the abrupt transition for a duplex composed by scale-free networks (degree distribution p ( k ) ∝ k − )with mean degree (cid:104) k (cid:105) = 3, minimum degree k = 2 and maximum degree k max ∼ O (cid:16) N / (cid:17) . In ( a ), average of the size of thelargest connected component as a function of the tolerance parameter. The error bars correspond to one standard deviation.In the inset, probability of S for α = 0 . , . , .
25 (see the video in the Supplemental Material for a better appreciation ofthe transition). The network size per layer is N = 1000. In ( b ), the Binder cumulant, Eq. (2), for different system sizes. Thedashed vertical line is located at the critical point. In the inset, position of its minima, with a linear fit (dashed line). In ( c ),probability of S , from N = 500 to 10000. In the inset, distance between peaks as a function of the system size. All points andhistograms are computed by at least 20000 independent realizations of the dynamics. feature of abrupt phase transitions. Note, moreover, that our model is genuinely out of equilibrium, the ergodicity isbroken due to an absorbing state, and the bimodality is not a signature of hysteresis or phase coexistence but remindsthe absence of self-averaging [62].To further characterize the transition we use the Binder cumulant [63] U N ≡ − (cid:104) S (cid:105) (cid:104) S (cid:105) , (2)where (cid:104) S n (cid:105) is the n th moment of the size of the largest connected component once the cascade has stopped. Thebehavior of this quantity for different system sizes provides the value of the critical point and can be used to infer thenature of the phase transition. In Fig. 2( b ) we show U N for different system sizes, and we can appreciate that theydisplay a non-monotonic behavior with a minimum that becomes lower and more narrow as N increases, which isagain a signature of an abrupt transition. The position of the minima can be fitted against the inverse of the systemsize [64], finding α c as the intercept of the fit [see the inset of Fig. 2( b )]. For the parameters used in the simulationsof Fig. 2, we obtain α c = 0 . ± . N → ∞ and is a bona fide discontinuous transition, we verify that, at the value of the tolerance for which the susceptibilityof S is maximum, the value of P ( S ) between the peaks vanishes as the size is increased [Fig. 2( c )], and that distancebetween the maxima ∆ S does not decrease [see inset of Fig. 2( c )] [65, 66]. Notice that ∆ S ( N ) ≈ .
85, hence thetransition is devastating and very hard to anticipate. In the light of the previous findings, we devote the rest of thearticle to study the phase transition precisely in scale-free multilayer networks, having in mind that they represent aworst-case scenario when compared to random regular graphs and Erd˝os-R´enyi networks.The dynamics of different sort of models, not only those regarding the robustness, suffer non-trivial changes depend-ing on whether the topology is layered or aggregated (single-layer). For example, when it comes to phase transitions,it is often observed in percolation and in first-neighbor cascade propagation that continuous transitions in monolay-ers become first-order when the dynamics is considered in the multilayer. This prompts us to compare the level ofrobustness and the type of transition for non-local cascade propagation in multiplex networks and their aggregatedrepresentation. For the aggregation process, we collapse all intralayer links into an aggregated network of N nodes.For the sake of simplicity, in the aggregated structure we do not consider multiedges, although this is not really aproblem as far as the number of layers of the original network is kept low, the networks are sparse and the systemsize is large enough, since the probability of observing an existing link in another layer scales as (cid:104) k (cid:105) /N . The firstobservation is that the aggregated network is more robust than the layered one, since (cid:104) S (cid:105) saturates to 1 for lowervalues of the tolerance [see Fig. 3( a )]. Put another way, the robustness of an aggregated network obtained from anoriginally layered structure will tend to be overestimated. Although the aggregated robustness is higher, the metricsto characterize the transition indicate that the system behaves as in the multiplex. Indeed the Binder cumulantbecomes negative in the critical region, with a narrower and deeper minimum as the system size increases [Fig. 3( b )],and the distribution of cascade sizes maintains the desired bimodality [Fig. 3( c )].The robustness overestimation effect can be understood from topological arguments. The aggregated network, byconstruction, is less heterogeneously connected than the multiplex one. It is well-known that degree heterogeneityamplifies the cascade propagation (see Ref. [9], but also our above results) hence we can expect a higher robustness inthe former. Moreover, aggregated networks are less sparse than layered ones: in the former the mean degree is (cid:104) k (cid:105) M ,while in the latter is (cid:104) k (cid:105) + M −
1. This sparsity is reflected in the average shortest path length (cid:96) , which is smaller foraggregated networks: for Fig. 3( a ) we obtain (cid:96) agg = 3 . ± .
05 and (cid:96) ML = 5 . ± .
1. We will see that (cid:96) correlateswell with the level of robustness of the network.Most real multilayer networks have been formed or constructed in a non-optimal way, in the sense that there existredundant layers to the overall functioning of the network [67–69]. This is natural, as multilayer systems usually arenot designed at once but their structure is a result of several growing and shrinking processes occurred at differenttime scales. This prompts us to ask how the robustness of a layered system depends on the number of layers. Wesee in Fig. 3( d ) that adding more layers causes an increase of robustness. This effect is not rooted in a change ofthe load heterogeneity due to the addition of layers: see the collapse of the load distributions in Fig. 3( e ). Instead,the global connectivity is enhanced as M increases, so the load excess caused by a perturbation can be absorbedmore easily by the rest of the network. This is reflected in the decreasing tendency of the average path length, seeFig. 3( f ). An apparent paradox is posed: why both collapsing the layers into an aggregated network and unfolding amultiplex network in an even more layered structure turn the system more robust? The answer is simply that thesetwo procedures generate less sparse networks, i.e., it is observed a decrease in the average path length (compare thetwo curves in Fig. 3( f )). Note, however, that here we are dealing with synthetic multiplex networks, with statisticallysimilar topologies across layers: same number of nodes N and same degree distribution P ( k ). Moreover, they havebeen constructed with no structural correlations, neither at the intralayer level nor between layers. Therefore, thefurther we are from this setting, the more uncontrollable effects might occur in the processes of aggregation or additionof layers and the less reliable (cid:96) might become to assess the robustness.To provide further evidence of the relation between the average path length and the robustness of the system, wenext study the case of partial multiplexity, which is another feature present in real networks and is known, as well, to . . . . Tolerance α . . . h S i ( a ) Agg.ML . . Tolerance α − B i nd e r U N ( b ) (Agg.) . . . . Tolerance α . . . h S i ( d ) L/M − − − P ( L ) M ( e ) M ‘ ( f ) AggregatedMultilayer
FIG. 3. In ( a ), comparison of the robustness between a duplex ( M = 2) of scale-free networks and its aggregated counterpart.In ( b ), Binder cumulant for aggregated networks of size N = 250 ,
500 and 1000. In ( c ), probability density function of thesize of the largest connected component at the end of the cascade, for aggregated networks. Values of the tolerance are α = 0 . , . , .
11 and 0 .
12 . In ( d ), robustness of a multiplex network as a function of its number of layers. In ( e ), loaddistribution for the networks in ( d ). In ( f ), average path length for multiplex networks formed with different number of layers M and for their aggregated counterparts. The network of each layer is constructed with the parameters indicated in the captionof Fig. 2. Each point in ( b ) and the histograms in ( c ) are computed from 20000 independent realizations of the dynamics. Allthe other results are averaged over 1000 realizations. strongly modify the dynamics of models running on multilayer networks [46, 48]. Let us focus in the scenario of twolayers, for simplicity. In this case, the multiplexity parameter q is the fraction of nodes that simultaneously participatein both layers, that is, they share the same state —functioning or failed— and are connected via an interlayer link,hence contributing to the load. This way, when q is low, the common nodes in the multilayer are highly loaded sincethey act as a bridge between the two layers. At the global level, higher multiplexity implies lower average path length,as it can be seen in Fig. 4( a ). Thus, for a fixed and finite N , we expect systems with large q to be more robust thanthose with a lower value. Indeed, as displayed in Fig. 4( b ), this is the behavior that we find: the value of the tolerancesuch that a cascade has barely affected the network ( S ∼
1) for q = 0 . q = 0 . . . . Multiplexity q ‘ ( a ) Tolerance α . . . h S i ( b ) . . . FIG. 4. In ( a ), average path length as a function of the multiplexity parameter. In ( b ), average of the size of the largestconnected component for different values q , indicated in the inset. Each point is averaged over 1000 independent realizationsof the dynamics, in networks with the same properties as those of Fig. 2. links can be interpreted as passenger layovers.Simulating the propagation of cascades in the air transportation multiplex and its aggregated version, we find that,as anticipated, the process of aggregation increases the robustness [see Fig. 5( a )]. Moreover, the transition in thelayered version of the network occurs at a finite tolerance and the jump is quite significant, agreeing with the expectedresult of a discontinuous transition.We also look at the role played by the multiplex parameter q . Recall that for a duplex, q represents the numberof nodes participating simultaneously in the dynamics and sharing the same state. When M >
2, however, newpossibilities unfold, since a node might be present in some layers but absent in other ones. Indeed, this occurs in ourexample of the European air transportation network: for a given airport, there are some companies that operate thereand others that do not, although technically they would be allowed to do so. We can generalize q in order to includethis scenario as follows. Let n i be the number of times node i is present (with intralayer neighbors) in a multiplexnetwork. In our case, this is the number of airlines operating in a given airport. The actual number of interlayerlinks between that node and its replicas is then n i ( n i − /
2. Let n T be the total number of possible interlayer linksof a node, even if it does not have intralayer connections, i.e., n T = M ( M − /
2. Then the multiplexity parameterfor such a network is defined as the ratio between the actual number of interlayer links and all the possible interlayerlinks, that is, q = 1 N M ( M − N (cid:88) i =1 n i ( n i − . (3)Note that the duplex case is included in this new definition.We can increase the q of a multiplex network by taking an inactive node of a layer (e.g, an airport in which acompany does not operate) and creating some intralayer connection(s) to it. This way that node becomes active inthe layer, permitting interlayer flow (e.g., passengers can take a layover from another layer and then fly using thatcompany). Although many strategies can be devised to increase q taking into account topological considerations, wedo so in the simplest possible manner: we randomly pick a node with no neighbors in a random layer and connectit to a node with neighbors on that layer, randomly selected as well. The more nodes we add, the higher will bethe multiplexity q . In Fig. 5( b ) we confirm that the air transportation network becomes more robust as q is raised,coinciding with the results obtained in synthetic networks. Certainly the robustness can be even further increased byapplying strategies other than random, but this is out of the scope of the present article. . . . . . α . . . S ( a ) AggregatedLayered . . . α . . . h S i ( b ) q FIG. 5. In ( a ), comparison of the robustness of the layered and aggregated versions of the European air transportation network.In ( b ), size of the LCC as the multiplexity is increased with the method explained in the main text. The curve for q = 0 . IV. CONCLUSIONS
In this article we have tackled the problem of cascades of failures that propagate in multiplex network in a non-localway. That is, the failure of a node can be induced by the failure of other nodes not necessarily linked to it. To thisgoal we define the rules of a load-capacity model in multilayer topologies. In such a model, nodes are assigned avariable load and a time-independent capacity, and they fail when, for some reason, their load exceeds the capacity,i.e., we encounter cascades of overload entities. Due to the non-local dynamics and the non-trivial topologies on whichthe dynamics runs, the computational approach to the problem is more accessible than an analytical treatment, thuswe offer a numerical characterization of the phenomena. The robustness of a network is related to purely structuralproperties, i.e., to what is the largest portion of the network that remains connected once the cascade stops.We have given compelling evidence of an abrupt phase transition, rooted in the intrinsic dynamics of non-localpropagation of the cascades. Hence, we have confirmed that non-local flow redistribution in multiplex networksshares the same fragility as the one observed in local spreading phenomena in these topologies, and, additionally,we have shown that this qualitative similarity between processes is broken for single-layer networks, being the non-local failure transition discontinuous too. To the best of our knowledge, this is the first time that such transitionis characterized for a load-capacity model in layered systems. We observe that this transition has a dependence onthe network topology. To shed some light on this relation, we work with the hypothesis that the amount of loadthat needs to be redistributed once a node fails is more evenly distributed across the network if the average shortestdistance between nodes is small. We check this hypothesis via three independent processes that induce a decrease inthe average path length: the collapse of a layered topology into a single layer, the addition of new layers, and theincrease of the common nodes participating simultaneously in the layers. In all cases, we find that the robustness isincreased. Finally, we have tested our predictions on the European air transportation multilayer network. There wehave shown that the transition is discontinuous, and that the aggregate version of the network overestimates the realrobustness of the original system. We have also defined a way to generalize the multiplexity parameter q to morethan 2 layers, and have shown that when it increases, the robustness does it as well.Our work is a first approximation to the study of non-local failure propagation processes in multiplex structures.As such, we have identified the interesting phenomenology that one can find, i.e., an abrupt transition, and we have0studied how to increase or decrease the robustness of the network by tuning a topological variable such as the averagepath length. However, several interesting questions and extensions have been left out of this first characterization ofoverload cascades on multiplexes, although we believe that our findings can serve to pave the way for further researchon this topic. Some of them are a better description of the phase transition in terms of the load distribution, whatis the impact of link weights on the critical properties, and a detailed analysis of the robustness for statisticallydissimilar networks in the different layers, as well as including structural correlations both in the intralayer and inthe interlayer sense. Note also that an underlying assumption of the model is that nodes communicate among themvia the shortest path, which is an approximation that might fail in some cases. Hence studying the same dynamicsof load redistribution where the load is computed with other centrality measures is certainly an interesting extensionfor the future. The overarching goal would be to develop an analytical theory of non-local propagation of overloadcascades that is able to support the results found computationally and flexible enough to encompass the open questionsmentioned above. [1] A.-L. Barab´asi, Linked: How everything is connected to everything else and what it means (Plume, 2003).[2] M. Newman,
Networks (Oxford University Press, 2018).[3] R. Albert, H. Jeong, and A.-L. Barab´asi, Nature , 378 (2000).[4] R. Cohen, K. Erez, D. Ben-Avraham, and S. Havlin, Physical Review Letters , 4626 (2000).[5] R. Cohen, K. Erez, D. Ben-Avraham, and S. Havlin, Physical Review Letters , 3682 (2001).[6] D. S. Callaway, M. E. Newman, S. H. Strogatz, and D. J. Watts, Physical Review Letters , 5468 (2000).[7] M. E. Newman, Physical Review E , 016128 (2002).[8] P. Holme, B. J. Kim, C. N. Yoon, and S. K. Han, Physical Review E , 056109 (2002).[9] A. E. Motter and Y.-C. Lai, Physical Review E , 065102 (2002).[10] A. E. Motter, Physical Review Letters , 098701 (2004).[11] S. Shao, X. Huang, H. E. Stanley, and S. Havlin, New Journal of Physics , 023049 (2015).[12] J. Gao, B. Barzel, and A.-L. Barab´asi, Nature , 307 (2016).[13] L. Zhao, K. Park, and Y.-C. Lai, Physical Review E , 035101 (2004).[14] X. Guan and C. Chen, Proceedings of the National Academy of Sciences , E8125 (2018).[15] R. Kinney, P. Crucitti, R. Albert, and V. Latora, The European Physical Journal B-Condensed Matter and ComplexSystems , 101 (2005).[16] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin, Nature , 1025 (2010).[17] Y. Yang, T. Nishikawa, and A. E. Motter, Science , eaan3184 (2017).[18] R. M. May, S. A. Levin, and G. Sugihara, Nature , 893 (2008).[19] R. H. Heiberger, Physica A: Statistical Mechanics and its Applications , 376 (2014).[20] M. Dubois, V. Rossi, E. Ser-Giacomi, S. Arnaud-Haond, C. L´opez, and E. Hern´andez-Garc´ıa, Global Ecology and Bio-geography , 503 (2016).[21] O. Artime, V. d’Andrea, R. Gallotti, P. L. Sacco, and M. De Domenico, arXiv preprint arXiv:2004.14879 (2020).[22] D. Stauffer and A. Aharony, Introduction to percolation theory , 3rd ed. (Taylor & Francis, 1991).[23] D. J. Watts, Proceedings of the National Academy of Sciences , 5766 (2002).[24] C. D. Brummitt, R. M. DSouza, and E. A. Leicht, Proceedings of the National Academy of Sciences , E680 (2012).[25] D.-S. Lee, K.-I. Goh, B. Kahng, and D. Kim, Physica A: Statistical Mechanics and its Applications , 84 (2004).[26] Y. Yu, G. Xiao, J. Zhou, Y. Wang, Z. Wang, J. Kurths, and H. J. Schellnhuber, Proceedings of the National Academy ofSciences , 11726 (2016).[27] North American Electric Reliability Council (NERC), 1996 System Disturbances. Review of Selected
Electric SystemDisturbances in North America , Tech. Rep. (2002).[28] NERC Steering Group,
Technical Analysis of the August 14, 2003, Blackout: What Happened, Why, and What Did WeLearn? , Tech. Rep. (2004) Report to the NERC Board of Trustees. [29] EUROCONTROL, “Ash-cloud of April and May 2010: Impact on Air Traffic,” (2010), Download from ; Accessed: 2020-03-17.[30] See the section “Airspace restrictions in Europe” in https://en.wikipedia.org/wiki/Air_travel_disruption_after_the_2010_Eyjafjallaj\unhbox\voidb@x\bgroup\let\unhbox\voidb@x\setbox\@tempboxa\hbox{o\global\mathchardef\accent@spacefactor\spacefactor}\accent127o\egroup\spacefactor\accent@spacefactorkull_eruption ; Accessed: 2020-03-17.[31] A. E. Motter and Y. Yang, Physics Today (2017).[32] P. Crucitti, V. Latora, and M. Marchiori, Physical Review E , 045104 (2004).[33] A. Moussawi, N. Derzsy, X. Lin, B. K. Szymanski, and G. Korniss, Scientific Reports , 11729 (2017).[34] M. Waniek, G. Raman, B. AlShebli, J. Chih-Hsien Peng, and T. Rahwana, arXiv preprint cs/2003.03723 (2020).[35] R. J. Baxter, Exactly solved models in Statistical Mechanics (Elsevier, 2016).[36] C. Castellano, S. Fortunato, and V. Loreto, Reviews of Modern Physics , 591 (2009).[37] O. Artime, A. Carro, A. F. Peralta, J. J. Ramasco, M. San Miguel, and R. Toral, Comptes Rendus Physique , 262(2019).[38] I. Dobson, B. A. Carreras, V. E. Lynch, and D. E. Newman, Chaos: An Interdisciplinary Journal of Nonlinear Science , 026103 (2007).[39] S. Pahwa, C. Scoglio, and A. Scala, Scientific Reports , 3694 (2014).[40] S. Ruan, in Mathematics for Life Science and Medicine (Springer, 2007) pp. 97–122.[41] S. Zapperi, P. Ray, H. E. Stanley, and A. Vespignani, Physical Review Letters , 1408 (1997).[42] I. Simonsen, L. Buzna, K. Peters, S. Bornholdt, and D. Helbing, Physical Review Letters , 218701 (2008).[43] M. De Domenico, A. Sol´e-Ribalta, E. Cozzo, M. Kivel¨a, Y. Moreno, M. A. Porter, S. G´omez, and A. Arenas, PhysicalReview X , 041022 (2013).[44] M. Kivel¨a, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, Journal of Complex Networks , 203(2014).[45] J. G´omez-Gardenes, I. Reinares, A. Arenas, and L. M. Flor´ıa, Scientific Reports , 620 (2012).[46] M. Diakonova, V. Nicosia, V. Latora, and M. San Miguel, New Journal of Physics , 023010 (2016).[47] M. De Domenico, A. Sol´e-Ribalta, S. G´omez, and A. Arenas, Proceedings of the National Academy of Sciences , 8351(2014).[48] O. Artime, J. Fern´andez-Gracia, J. J. Ramasco, and M. San Miguel, Scientific Reports , 7166 (2017).[49] B. Min, S. Do Yi, K.-M. Lee, and K.-I. Goh, Physical Review E , 042811 (2014).[50] G. Bianconi, Journal of Statistical Mechanics: Theory and Experiment , 034001 (2017).[51] F. Radicchi and G. Bianconi, Physical Review X , 011013 (2017).[52] C. D. Brummitt, K.-M. Lee, and K.-I. Goh, Physical Review E , 045102 (2012).[53] K.-M. Lee, C. D. Brummitt, and K.-I. Goh, Physical Review E , 062816 (2014).[54] S. D. Reis, Y. Hu, A. Babino, J. S. Andrade Jr, S. Canals, M. Sigman, and H. A. Makse, Nature Physics , 762 (2014).[55] M. Turalska, K. Burghardt, M. Rohden, A. Swami, and R. M. D’Souza, Physical Review E , 032308 (2019).[56] D. Zhou and A. Elmokashfi, PloS one , e0189624 (2017).[57] P. H. Nardelli, N. Rubido, C. Wang, M. S. Baptista, C. Pomalaza-Raez, P. Cardieri, and M. Latva-aho, The EuropeanPhysical Journal Special Topics , 2423 (2014).[58] A. Cardillo, J. G´omez-Gardenes, M. Zanin, M. Romance, D. Papo, F. Del Pozo, and S. Boccaletti, Scientific Reports ,1344 (2013).[59] R. M. D’Souza, J. G´omez-Garde˜nes, J. Nagler, and A. Arenas, Advances in Physics , 123 (2019).[60] Y. Kornbluth, G. Barach, Y. Tuchman, B. Kadish, G. Cwilich, and S. V. Buldyrev, Physical Review E , 052309 (2018).[61] M. Barthelemy, The European Physical Journal B , 163 (2004).[62] S. Wiseman and E. Domany, Physical Review Letters , 22 (1998).[63] K. Binder, Zeitschrift f¨ur Physik B Condensed Matter , 119 (1981).[64] M. M. de Oliveira, M. da Luz, and C. E. Fiore, Physical Review E , 060101 (2018).[65] P. Grassberger, C. Christensen, G. Bizhani, S.-W. Son, and M. Paczuski, Physical Review Letters , 225701 (2011).[66] L. Tian and D.-N. Shi, Physics Letters A , 286 (2012). [67] M. De Domenico, V. Nicosia, A. Arenas, and V. Latora, Nature Communications , 6864 (2015).[68] L. Lacasa, I. P. Mari˜no, J. Miguez, V. Nicosia, ´E. Rold´an, A. Lisica, S. W. Grill, and J. G´omez-Garde˜nes, Physical ReviewX , 031038 (2018).[69] A. Ghavasieh and M. De Domenico, Physical Review Research2