The “weak” interdependence of infrastructure systems produces mixed percolation transitions in multilayer networks

Previous studies of multilayer network robustness model cascading failures via a node-to-node percolation process that assumes “strong” interdependence across layers–once a node in any layer fails, its neighbors in other layers fail immediately and completely with all links removed. This assumption is not true of real interdependent infrastructures that have emergency procedures to buffer against cascades. In this work, we consider a node-to-link failure propagation mechanism and establish “weak” interdependence across layers via a tolerance parameter α which quantifies the likelihood that a node survives when one of its interdependent neighbors fails. Analytical and numerical results show that weak interdependence produces a striking phenomenon: layers at different positions within the multilayer system experience distinct percolation transitions. Especially, layers with high super degree values percolate in an abrupt manner, while those with low super degree values exhibit both continuous and discontinuous transitions. This novel phenomenon we call mixed percolation transitions has significant implications for network robustness. Previous results that do not consider cascade tolerance and layer super degree may be under- or over-estimating the vulnerability of real systems. Moreover, our model reveals how nodal protection activities influence failure dynamics in interdependent, multilayer systems.

The robustness of a complex networked system to survive random component failures and/or intentional attacks is a significant problem with broad implications. Robustness, i.e., the likelihood that a network remains functional after losing nodes or links, has been a topic of interest since the beginning of modern network science [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] . In recent years, the robustness of interdependent, multilayered networked systems (sometimes referred to as networks-of-networks) has become a subject of focused study [20][21][22][23][24] . Examples of real-world interdependent, multilayer networks are ample, with a common and important case being urban infrastructure systems 25 such as transportation, communication, electric power, and water supply networks. The sharing of services within and between these infrastructure systems suggests that the loss of a single service like mobility can impact the provision of others like electric power and clean water-when a node or link in one infrastructure network fails, because of dependencies across infrastructures, its neighbors in other networks can also fail 26 . Since failures can propagate from one network to others, understanding the ways in which dependencies within and across these multilayer networks influences systemic robustness is a crucial first step to develop broadly influential resilience recommendations. The goal of this paper is to develop a model for the failure dynamics which is more realistic to interdependent infrastructure contexts than those in previous studies. To make the terminologies unambiguous, we use the multilayer network lexicon developed by Kivelä et al. 23 to describe our approach.
The phenomenon of failure propagation or spreading in a multilayer network depends on the dynamical or physical processes specific to the type of node or link failures. It is difficult to develop a general, dynamics based framework to address the robustness and resilience of systems like urban infrastructure where nodes, links, and layers may represent characteristically different objects and relationships. A viable approach is then to focus on

Model
We consider a percolation process on a system of M layers of networks A, B, C, …, each having the same number N of nodes. The networks form an interdependent, multilayer system with a possible structure shown in Fig. 1. A node in the system can be denoted by a pair of labels (x, X), with x(x = 1, 2, … N) and X(X = A, B, C…) being the nodal and layer indices, respectively. Nodes within the same layer X are connected with the degree distribution p k X . Nodes in different layers sharing a common nodal index are the replica nodes connected by dependency links if the underlying networks to which they belong are interdependent upon each other. That is, for every nodal label x, the set of M replica nodes corresponding to pairs of labels (x, X) can also be regarded as possessing a network structure through the dependency links. The number of interdependent networks, q X , for network X is the "super degree, " where every node in X has q X interdependent replica nodes. A randomly connected network layer X with degree distribution p k X is essentially a set of connected components 58 , where the nodes in the giant component are regarded as functional or viable and the remaining nodes are treated as failed or inviable. Our node-to-link failure propagation process can be described, as follows. Consider a pair of interdependent networks: A and B. If a node (x, A) in A fails, each link of its replica (x, B) in network B will be disabled with probability 1 − α. Similarly, if a node (y, B) in network B fails, the links of its replica (y, A) in network A will be cut off with the same probability. A initial failure due to isolation from the giant component in a certain network can then spread across the whole system through an iterative process. In each iteration, disconnecting certain nodes from the giant component of network A will cause some nodes to be isolated from the giant component of network B, which in turn will induce more nodal failures in A. This recursive or cascading process occurs in all pairs of interdependent networks synchronously at each iteration. When no failures are possible, the whole system reaches a stable steady state.
From the perspective of a single node, the parameter α quantifies the tolerance to its failed interdependent replicas, which controls the impacts that a node will endure if one of its replicas fails and represents the probability that interdependent systems will be unsuccessful at preventing inter-layer cascades. From the standpoint of the whole networked system, α determines effectively the strength of interdependence among the networks. For α → 1, the interdependence between nodes is the weakest so that, practically, failures cannot spread from one network to another. The opposite limit α → 0 signifies the case where interdependence is the strongest. In this case, our model reduces to previous percolation models with the node-to-node failure propagation mechanisms 35 . In our study, we use the size of the giant components S X for the final network layers X (X = A, B, …) to measure the robustness of the system, as in previous works 21 .
We consider our percolation model to be a simplified representation of urban infrastructure systems with buffers against interdependent infrastructure failures. A representative system following the general form of our model is multi-modal transportation linking multiple cities. To convert Fig. 1 into a transportation model, we would represent each colored node as a separate city and each layer as a different mode of transportation. Within cities, people travel via urban transportation (e.g., cars and buses) to different kinds of hub locations (e.g., coach stations, airports, and railway stations) to facilitate mode switching and interactions. One transportation layer is then linked via one kind of hubs in cities, as connecting highways, planes or railways link urban regions separated by large geographic distances. Although the current model is insufficient to consider specific travel dynamics within each city, our model provides a general form to approach the question of how dynamic buffers may lead to cascades within and across modes of transportation. In particular, the tuning of α from weak to strong interdependence captures situations where the hubs in different modes of travel have either an excess or a lack of capacity to handle the congestion caused by additional passengers, respectively. This model is general, such that cascade tolerance can be due to physical constraints (e.g., number of highway lanes), temporal constraints (e.g., early morning vs. rush hour traffic), and/or dynamic actions (e.g., lane shifting and emergency procedures). Thus, although simplified, robustness analysis of the multilayer system is informative to the ways in which transportation across urban regions may cascade and disconnect interdependent traffics. Moreover, the general form of our model can also capture buffers across other systems by treating layers as different infrastructures that link to transportation hubs within and across cities (e.g., electric power transmission substations).

Theory
General formalism. We solve our model analytically in terms of the final state after a cascading process using the standard method of generating functions [59][60][61] . Let R X be the probability that a randomly chosen link in network layer X belongs to the giant component, for is the generating function that gives rise to the degree distribution for random nodes in layer X, and G x is the generating function for the underlying branching processes of layer X, which generates the distribution for the number of outgoing links of randomly chosen links in layer X, where 〈k〉 X is the average degree of network X.
For simplicity, we consider the case where the M network layers within the multilayer system have an identical degree distribution: , and 〈k〉 X ≡ z to simplify the notations. Assuming that the network of M network layers has no loop, we aim to obtain the equation governing the probability S X that a random node in a given layer X belongs to the giant component -the viable probability for this node. Suppose the node has t failed replicas. Each link of this node is preserved with the probability α t R X . The viable probability for this node is then 1−∑ k p k (1−α t R X ) k . Let f X (t) be the probability distribution for the number of failed t replicas of a random node in network layer X. We get the probability S X in terms of the generating function G 0 as From the probability distribution f X (t), we can also derive the equation for R X based on the branching process in network X. Following a randomly chosen link in layer X, we arrive at a node (x, X) of degree k, where k follows the distribution kp k /z. With t failed replicas, each link of the node (x, X) is preserved with the probability α t R X . The probability that node (x, X) has at least one outgoing link to the giant component is 1−(1−α t R X ) k−1 , which is also the probability that a randomly chosen link can lead to the giant component. Using the probability distribution f X (t) for a random node in layer X, we have To solve Eqs (1) and (2), we need the probability distribution function f X (t). We obtain a solution based on the level-by-level calculating process with respect to a hierarchical structure 62-64 , as illustrated in Fig. 2. In particular, network layer X is assigned to the top level, whose nearest neighbors constitute the 2 nd level, and the neighbors of the nearest neighbors belong to the 3 rd level, and so on. This way, for a given layer Z with a super degree q Z , it has q Z neighbors in the lower level if it is at the top level. Otherwise it has q Z − 1 neighbors in the lower level.
We first calculate the viable probabilities for the nodes in the bottom-level layers. For a random node in such a network layer Z, it has no lower-level replicas. The probability distribution function f Z (t) for the number of failed replicas of a random node is given by The viable probability S Z for a random node in network Z is where q Z − 1 is the number of neighboring networks in the lower level of network Z. After calculating the viable probabilities for the nodes in all bottom-level network layers, we can get the probability distribution f Y (t) for a random node in a higher level network layer Y. Substituting f Y (t) into Eq. (4), we can get the the viable probability S Y . Repeating this process level by level, we can get the probability distribution function f X (t) for the top level.

Star-like multilayer networks.
To be concrete, we consider a star-like interdependent system of four network layers, as shown in Fig. 2(a), where layer B is the hub and the other three layers A, C, D are peripheral networks (on the same footing). We thus have We first calculate the viable probability S B for a randomly chosen node in the hub layer B. The viable probability for a node in one of the peripheral (bottom-level) network layers is 1 − G 0 (1 − R), and the probability distribution f B (t) for a random node in B to have t failed replicas is binomial: (1), we obtain the probability for a random node in layer B to belong to its giant component: Next, we calculate the viable probability for a random node in network A, which depends the state of its replica in layer B. If the replica is viable, we have t = 0; Otherwise t = 1. According to the updating process in Fig. 2, the viable probability of a node in layer B is determined by the number t of its failed replicas in the other two peripheral network layers (excluding layer A), which is The inviable probability for a random node in B is thus given by (1), we obtain the viable probability for a random node in layer A as The same equation holds for S C and S D .
To solve Eqs (5) and (6) requires the equations for the quantities R and R′, which we get through the branching process in network layers A and B, respectively. In particular, based on the probability distribution function f B (t), we arrive at a self consistency equation in terms of the generating functions: Similarly, the self consistency equation for R is It is not feasible to get closed-form expressions for R′ and R from Eqs (7) and (8). We thus resort to numerical solutions. For a given degree distribution p k and a fixed value of the tolerance parameter α, we plot the curves for Eqs (7) and (8)  Substituting the values of R′ and R into Eqs (5) and (6), we can get the sizes of the giant components S A and S B . Figure 3 shows, for random networks with a Poisson distribution 65 p k = e −z z k /k!, graphical solutions of R and R′ for different values of α and z. For α = 0.9, there is a trivial solution at the point (R = 0, R′ = 0) for z = 1.2, indicating that the hub and the peripheral network layers are completely fragmented. For z = 1.3, a nontrivial value of R arises but R′ remains to be zero, suggesting a continuous phase transition for R. Increasing z to 1.5, we observe that R′ also becomes nontrivial which changes growth rate for R. These results imply that, for α = 0.9, the peripheral network layers percolate firstly, and then the hub percolates as z is increased, leading to double phase transitions for the peripheral network layers. For α = 0.7, the phenomenon of mixed percolation transitions occur, where the peripheral network layers percolate in a continuous manner, as shown in Fig. 3(e), but the hub percolates in a discontinuous fashion, as demonstrated by the tangent point shown in Fig. 3(f). Specifically, for z < 2.2348, the solutions for (R, R′) are given by the crossing point of the curves with the R-axis. However, for z ≈ 2.2348, their solutions are given by the tangent point (0.5367, 0.3403), giving rise to a discontinuous change in both R and R′. For α = 0.6, we find that the crossing point for solutions of R and R′ change abruptly from (0, 0) to the tangent point (0.5874, 0.4253), indicating also a discontinuous percolation transition.
The sizes S A , S B , S C , and S D of the giant components for the four network layers in the star-like system versus the average degree z are shown in Fig. 4. (See Sec. 4 for an explanation of the numerical methods.) For α = 0.9, as z is increased from the value of one, all four layers exhibit a continuous percolation transition, where the peripheral layers (layers A, C, and D) percolate first at the critical point z c1 ≈ 1.2345 [ Fig. 4(a)], followed by the hub network (layer B) at the critical point z c2 ≈ 1.4872, as shown in Fig. 4(b). We see that, while the curves of S A , S C , and S D versus z are continuous, they exhibit a discontinuity in their derivative at both z c1 and z c2 . In this sense we say that the peripheral network layers exhibit double percolation transitions. Note that, in the whole range of z values considered, the hub layer exhibits only a single percolation transition -continuous (or second order) for relatively large values of α but discontinuous (of first order) for smaller α values. The feature of double percolation transitions for the peripheral network layers persists for α = 0.8 and α = 0.7. Remarkably, for these two values of α, the two percolation transitions for the peripheral networks are characteristically different: the first one is continuous while the second is discontinuous. In contrast, the transition for the hub network layer is discontinuous. Thus, from the standpoint of the whole interdependent, multilayer system, for these values of α, both continuous and discontinuous percolation transitions exist, leading to the phenomenon of mixed percolation transitions. For α = 0.6, the hub and peripheral layers percolate discontinuously at the same point. These results indicate a richer variety of percolation transition scenarios than revealed in previous studies for multilayer systems with strong interdependence (i.e., α = 0) 35 .
A practical implication is that the tolerance parameter can be exploited for modulating or controlling the characteristics of the percolation transition. In particular, for relatively small values of α (e.g., α = 0.6), the percolation transitions are abrupt and discontinuous -the defining characteristic of first-order phase transitions. In this case, the interdependent system will collapse suddenly as a system parameter is changed in response to random nodal failures or intentional attacks. The system is thus not resilient. To improve the resilience of the system, a larger value of α can be chosen (e.g., α = 0.9). In this case, the percolation transitions in the interdependent networks are continuous -signature of second-order phase transitions. While the whole system still collapses, the manner by which the collapse occurs is benign as compared with that of first-order phase transitions. Tree-like multilayer networks. We consider a tree-like interdependent system of five random network layers, labeled as A, B, C, D and E respectively, as shown in Fig. 2(b), and obtain theoretical solutions for the sizes of the giant components and for the nodal viable probabilities for all the layers. Figure 5 shows the theoretical and numerical solutions of the sizes of giant components S A , S B , S C , S D and S E versus the average degree z (see Sec. 4 for an explanation of the numerical methods). For α = 0.9, we find that the peripheral layers A, C, E percolate first, followed by the sub-center layer D. The size of the giant component of D's nearest neighboring layer E is then increased, leading to a continuous phase transition for S E . The central network layer B percolates last, giving rise to a distinct continuous transition for its nearest neighboring networks A, C, D. As α is decreased, the hub network percolates discontinuously, at which the giant components of layers A, C, E and D increase abruptly in size, giving rise to the phenomenon of mixed percolation transitions, as can be seen from the curves for α = 0.8. For α = 0.7, the percolation transition points for some network layers with large super degrees begin to merge but the phenomenon of mixed percolation transitions persists. As α is decreased further, a first-order phase transition occurs at which all layers percolate at the same point in a discontinuous manner, as indicated by the curves for α = 0.6.
These transition phenomena suggest that the super degree of a layer in the multilayer system plays a critical role in its percolation transition. In particular, the peripheral network layers with the lowest super degree percolate first, followed by the layers with a moderate value of the super degree, and finally by the layers with the highest super degree. When a layer with a larger super degree percolates continuously, it will increase the giant component sizes of its neighboring network layers that have already percolated, leading to multiple percolation transitions. In contrast, if a layer with a larger super degree percolates discontinuously, it will lead to a sudden and discontinuous increase in the sizes of the giant components of all network layers that have already percolated, leading to mixed percolation transitions. These results indicates that the percolation type of the hub layers plays a critical role in the robustness of the whole system. In particular, if the hub layer percolates continuously, at the transition point the sizes of the giant components of the nearest network layers are continuous but their derivatives are discontinuous. However, if the hub layer percolates discontinuously, an abrupt and discontinuous increase in the giant component sizes of the neighboring network layers can occur.

Numerical verification and effects of interdependence heterogeneity on percolation transition
Mixed and multiple percolation transitions. To verify the phenomena of mixed and multiple percolation transitions directly, we perform bond percolation 20 times using the Newman-Ziff algorithm 66 and measure the average relative size of the viable component of each network layer. As shown in Figs 4 and 5, the percolation transitions and the behaviors of the giant components as predicted theoretically agree well with the numerical results. Figure 6(a) and (b) show, for the star-like and tree-like interdependent systems, respectively, the percolation transition points z c versus α for each layer in the multilayer system, which are obtained by the graphical solution as illustrated in Fig. 3. An alternative way to identify the transition points is to examine the fluctuations in the size of the giant component which, for a finite system, tend to be relatively large at the transition 42 . From Fig. 6(a), we

Heterogeneity in interdependence.
In realistic situations, heterogeneity in interdependence across the full network can be expected. In the pioneering works on multiplex networks 35,36 , the authors considered unbalanced interdependence, where both strong and weak interdependence can arise in opposite directions of network interaction. In our model, the degree of interdependence is characterized by the parameter α for 0 ≤ α ≤ 1, where α = 0 and α = 1 correspond to the two extreme cases of the strongest possible and the weakest possible interdependence, respectively. Mathematically, heterogeneity in interdependence can then be modeled by a probability distribution of α. For example, for the star-like interdependent system analyzed in detail in our work, we can assume that the hub network has a dependence on each peripheral network determined by the tolerance parameter α 1 , while the dependence of each peripheral network on the hub network is characterized by another tolerance parameter α 2 , with α 1 ≠ α 2 . The difference |α 1 − α 2 | thus quantifies the degree of heterogeneity of interdependence. In this setting, the equations determining the sizes of the giant components of the hub network B and any peripheral network, e.g., A, can be modified to and Similarly, the quantities R′ and R can be obtained based on the branching processes by using Eqs (9) and (10), respectively, as Figure 7 presents results from direct simulation and theory. We find that the phenomenon of mixed percolation persists in the presence of interdependence heterogeneity.
Another form of heterogeneity we consider to relate our model to a broader set of multilayer infrastructure systems is in the number of interdependent links between network layers. We study this form of heterogeneity through the construction of a general model with more than one interdependent link for each node based on the "weak" interdependence mechanism. In particular, for the star-like interdependent system, this situation can be modeled as follows. We first choose a random node in network A and another random node in network B, and place an interdependent link between them. We then repeat this step until Nl AB interdependent links have been placed, where l AB denotes the average interdependent link degree for nodes in network A or B. Finally, for the other two pairs of interdependent networks, namely, C and B as well as D and B, we add the same number of interdependent links. In this revised model, the interdependent link degree for the nodes in each layer follows the Poisson distribution with the average degree l AB = l CB = l DB = m. Figure 8 shows the simulation results on percolation transisions, which exhibit features characteristic of mixed transitions, giving further support for the robustness of the phenomenon.

Discussion
The results of this paper demonstrate that the pervasive assumption of "strong" interdependence in percolation models is limiting our understanding of the robustness of multilayer systems. When tuning the interdependence between networks via a cascade tolerance mechanism, the surprising phenomenon of mixed percolation transitions occurs: network layers with large super degrees percolate discontinuously while those with small super degrees percolate continuously. That is, multilayer networks with weak interdependence across layers, represented by moderate values of the tolerance parameter, experience both second-order and first-order percolation transitions. This result is unlike those in any study that assumes either strong or no interdependence between network layers, where all network layers experience either simultaneous and abrupt percolation or continuous  percolation, respectively. These previous models imply that catastrophic cascades were either inevitable or impossible, yet our results demonstrate that both the strength of interdependence and layer position play a critical role in the functioning of interdependent systems.
In cases where the connectivity of the central layer is critical to the functioning of the interdependent system, sufficient tolerance must be included to prevent sudden, system-wide failures. This result implies that studies that do not consider interdependence may underestimate the amount tolerance needed to ensure cascades do not occur. On the other hand, where layers with high super degree only serve as bridges to critical periphery layers, far less tolerance is needed to maintain interdependent function. Likewise, studies that assume strong interdependence may be underestimating the likelihood that the system will survive. In both cases, the under-and over-estimation of multilayer network robustness may produce unwanted cascading dynamics if used to study and design real networks.
Since cascade tolerance provides a generic metric for the physical, temporal, and dynamic buffers found in real infrastructure to prevent losses across interdependent systems, the results of this paper have broad implications for infrastructure design. Infrastructure design often uses risk analysis to give the individual assets that comprise large-scale infrastructure systems tolerances to known failure modes. For example, individual assets have an oversized capacity to handle sudden shifts in service flows, have uninterruptible power supplies to ensure asset functioning during blackouts, and can be built a certain number of feet above the ground to prevent flooding. Despite significant work protecting individual assets in infrastructure systems, there is far less understanding how risk adverse practices in a single facility may impact the robustness of an entire infrastructure system or interdependent, multilayered system. The tolerance parameter introduced in our model is generic to capture the probability that these risk adverse practices influence percolation dynamics in a wide variety of infrastructure systems. Although previous models would indicate that interdependent systems could experience more sudden, global cascades, mixed percolation transitions may be an indicator of an even more precarious situation. With buffers that appear sufficient within a single network layer, a central layer with high super degree can still experience sudden, catastrophic failures that inhibits functionality of the multilayer network. This situation is common in disaster situations, where individual assets and entire infrastructure systems may survive the initial failure event, but remain unusable because interdependent systems that they require to function did not. The values of α that generate mixed percolation transitions in our model provide a heuristic measure of when these interdependent failures may occur that is characteristically different from other percolation models. Thus, our results suggest that developing a metric for buffering capacity and identifying the super degree relationships among infrastructure systems may reveal a disconnect between local, risk adverse practices, and global, interdependent vulnerability.
Although mixed percolation transitions have important implications for future percolation models and infrastructure design, we do concede that our current model is too simplified to represent many real-world systems. Percolation dynamics provide an important, generic method for considering the interactions within multilayer networks. However, the dynamics that occur within real-world systems are not captured in this model, suggesting that our assumptions for tolerance and cascading dynamics remain unrealistic. For example, modeling the failure dynamics in a multi-modal transportation network like the one discussed above (see Section 2) requires detailed information regarding the infrastructure located within each city and models for how people choose between different transportation modes. The current percolation model assumes that the loss of an airport or railway station in a periphery layer affects losses in interdependent layers equally. However, real travel information may show greater heterogeneity in congestion and failures. Heterogeneity can be introduced by assuming that the hub network layer depends on each peripheral with a tolerance parameter α 1 , while each peripheral layer depends on the hub layer with another tolerance parameter α 2 in, e.g., a star-like system. Our theoretical and simulation results show that the phenomenon of mixed percolation persists in the presence of heterogeneous interdependence. A deficiency our model is that it cannot depict the system in which one node has several interdependent partners in another layer. We also extend our model by allowing the nodes to have multiple interdependent links, with simulation results supporting the emergence of mixed percolation. For either heterogeneous or multiple interdependence or both, exactly capturing the detailed interdependence across different infrastructures needs empirical data. Future work should focus on linking discipline-specific infrastructure models that incorporate these dynamics with measures of cascade tolerance across network layers to produce more realistic percolation transitions.