Redundant interdependencies boost the robustness of multilayer networks

In the analysis of the robustness of multiplex networks, it is commonly assumed that a node is functioning only if its interdependent nodes are simultaneously functioning. According to this model, a multiplex network becomes more and more fragile as the number of layers increases. In this respect, the addition of a new layer of interdependent nodes to a preexisting multiplex network will never improve its robustness. Whereas such a model seems appropriate to understand the effect of interdependencies in the simplest scenario of a network composed of only two layers, it may seem not suitable to characterize the robustness of real systems formed by multiple network layers. It seems in fact unrealistic that a real system, evolved, through the development of multiple layers of interactions, towards a fragile structure. In this paper, we introduce a model of percolation where the condition that makes a node functional is that the node is functioning in at least two of the layers of the network. The model reduces to the commonly adopted percolation model for multiplex networks when the number of layers equals two. For larger number of layers, however, the model describes a scenario where the addition of new layers boosts the robustness of the system by creating redundant interdependencies among layers. We prove this fact thanks to the development of a message-passing theory able to characterize the model in both synthetic and real-world multiplex graphs.


I. INTRODUCTION
Multilayer networks [1][2][3][4][5] are emerging as a powerful paradigm for describing complex systems characterized by the coexistence of different types of interactions or coupling among different types of networks. Multilayer networks represent an appropriate descriptive model for real networked systems in disparate contexts, such as social [6,7], technological [8][9][10] and biological systems [11][12][13] , just to mention a few of them. For example, global infrastructures are formed by several interdependent networks, such as power grids, water supply networks, and communication systems [4]. Cell function and/or malfunction (yielding diseases) cannot be understood if the information on the different nature of the interactions forming the interactome (protein-protein interactions, signaling, regulation) are not integrated in a general multilayer scenario [11]. Similarly, the complexity of the brain is encoded in the different nature of the interactions existing at the functional and the structural levels [12,13].
General multilayer networks as well as multiplex networks [1][2][3][14][15][16] are composed by nodes belonging to different layers and links connecting nodes within and across layers. In general multilayer networks there are no restriction in the way the links across different layers can be placed. In multiplex networks, however, the nodes of each pair of layers are mapped one to one and the links across different layers can only be present among corresponding nodes. Therefore multiplex networks are a specific case of multilayer networks with a well defined struc-ture. In multilayer networks, and multiplex networks as well, nodes belonging to different layers are often interdependent on each other, in the sense that a failure of one node might cause the failure of a node in another layer. For example in global a infrastructure network, a power plant might be interdependent on a node of the communication system that is controlling its function [4].
Percolation models are generally used as proxies to quantify the robustness of networked systems under local failures, by monitoring how the connectedness at the macroscopic level changes as a function of the amount of microscopic damages of the individual elements of the system. Although different percolation models can be suitably defined and studied on multilayer networks (e.g. k-core percolation [17], weak percolation [18], and bond percolation [19]), here we focus our attention on the case of the so-called site-percolation model, where the units that can be potentially damaged are the nodes of the network and the order parameter of interest is the macroscopic connectedness of the system. When interdependencies are taken into account, the resulting percolation theory [4,[20][21][22][23][24][25][26][27][28][29] provides a general framework to characterize the dramatic avalanches of failure that can affect multilayer networks. A large volume of publications has pointed out that, according to this model, a multilayer network is much more fragile than each of the various network layers taken in isolation [4,[20][21][22][23][24][25]. In particular, the overall fragility of the system increases as the number of layers increases [23,[26][27][28][29]. Such a feature has an intuitive explanation. In this percolation model, a node is damaged if at least one of its interdependent nodes is damaged. As the number of layers increases, the probability of individual failures grows thus making the system more fragile. This scenario leads, however, to the conundrum: if the fragility of a system is increased by the number of layers of interactions, why are there so many real systems that display multiple layers of interactions? Further, the addition of new layers of interactions in a pre-existing multilayer network has generally a cost, so it doesn't seem reasonable to spend resources just to make the system less robust. For instance, according to a recent study focusing on diffusion dynamics on multimodal transportation networks [9], the presence of multiple interconnected modes of transportation makes the system more navigable and more robust than its individual layers considered in isolation. However, the aforementioned conundrum has not been fully addressed in terms of purely topological properties. The purpose of the current paper is to provide a potential explanation by introducing a suitable model for percolation in multiplex networks composed of multiple interacting layers. In the model, we will assume that a node is damaged only if all its interdependent nodes are simultaneously damaged. The model is perfectly equivalent to the one currently in use when the number of layer equals two. Additional layers, however, provide the system with redundant interdependencies, generating backup mechanisms against the failure of the system, and thus making it more robust.
Here, these effects are observed despite of the presence of interdependencies. We provide a comprehensive study of percolation in presence of redundant interdependencies thanks to the development of an exhaustive message-passing theory [20,[33][34][35][36] (also known as the cavity method). We build on recent advances obtained in the "standard" percolation theory for multilayer networks [26,[37][38][39][40] to propose a theory that is valid for arbitrary systems (thus also including overlap among layers [7,16] ), as long as the network structure is locally tree-like. This limitation is common to all message-passing approaches for studying critical phenomena on networks. Corrections have been recently proposed [41] on single networks to improve the performance of message-passing theory and similar approximations valid for loopy multilayer networks might be envisaged in the future. We remark that our model represents a starting point to address an obvious-yet-neglected feature that makes percolation more realistic as a model to study the robustness of real multilayer networks. Eventual modifications and/or the addition of further ingredients to the model presented here may be still necessary to deal with specific scenarios to make the model even more realistic.

II. REDUNDANT PERCOLATION MODEL ON MULTIPLEX NETWORKS
A multilayer network structure is not equivalent to a large single network. As a network is ultimately a way to encode information about a complex system, distinguishing between different types of links and nodes may significantly alter the characteristics of the structure and the dynamical behavior of the system. This fact is particularly evident when we associate a different role to links within each layer (intralinks) and links across different layers (interlinks). For example when describing the diffusion in a multilayer networks [42], we might reasonably assume that the diffusion constant along intralinks is different from the diffusion constant among interlinks and this changes significantly the typical timescale of the dynamics. Similar considerations are also valid for spreading processes in multilayer networks [43]. In the context of percolation theory, attributing a different role to interlinks and intralinks can yield a scenario significantly different from percolation in single layers. Specifically, if interlinks describe interdependencies between the nodes, the percolation transition becomes discontinuous and hybrid [4,21,25], and close to the percolation transition the multilayer network is affected by large avalanches of cascading failures. In this case, we will refer to the multilayer network as a set of interdependent networks. Therefore "interdependent networks" is a term that refers specifically to the response of the system to the damage of the nodes more than to the actual structure of the multilayer network.
Depending on type and number of interlayer connections present in the system, different classes of multilayer networks [1][2][3]14] can be considered. Here, we deal with one of the simplest classes characterized by the fact that every node is connected to one and only one node in every of the other layers. These systems are generally named as "multiplex networks" [1][2][3][14][15][16]. Very often the linked nodes across different layers (also called replica nodes) actually describe different realization of the same node. For example in the London transportation network, replica nodes could represent Oxford Circus tube station and Oxford Circus bus station.
However, multiplex networks can be also used to indicate the scenario where nodes belonging to the various layers represent physically distinct units as long as the nodes of the different layers are mapped one-to-one and the links cross layers are placed among all the corresponding (replica) nodes and nowhere else. For instance multiplex networks can be used to model interdependent infrastructures when the interdependencies occur exclusively between replica nodes. Systems of this type are interdependent multiplex networks [4,21,22,44] also named one-to-one interdependent networks.
We consider a multiplex network composed of M layers G α with α = 1, 2, . . . , M . We indicate the set of all layers as G = (G 1 , G 2 , . . . , G M ). Every layer contains N nodes. Exactly one node with the same label appears in every individual layer. Nodes in the various layers sharing a common label are called replica nodes, and they are considered as interdependent on each other. Nodes in the network are identified by a pair of labels (i, α), with i = 1, 2, . . . , N and α = 1, 2, . . . , M , the first one indicating the index of the node, and the second one standing for the index of the layer. For every node label i, the set of replica nodes is given by the M nodes corresponding to pairs of labels (i, α) with with α = 1, 2, . . . , M (see Figure 1). When at least two replica nodes (i, α) and (i, α ) are connected to two corresponding replica nodes (j, α) and (j, α ) we say that the multiplex network displays link overlap. Given a multiplex network as described above, we consider a percolation model where some of the nodes are initially damaged. We assume that interlinks represent interdependencies among replica nodes, but we consider the case in which such interdependencies are redundant, i.e., every node can be active only if at least one its interdependent nodes is also active. We refer to this model as "Redundant percolation model." As an order parameter for the model, we define the so-called Redundant Mutually Connected Giant Component (RMCGC). The nodes that belong to the RMCGC can be found by following the algorithm: (i) The giant component of each layer α is determined, evaluating the effect of the damaged nodes in each single layer; (ii) Every replica node that has no other replica node in the giant component of its proper layer is removed from the network and considered as damaged; (iii) If no new damaged nodes are found at step (ii), then the algorithm stops, otherwise it proceeds, starting again from step (i).
The set of replica nodes that are not damaged when the algorithm stops belongs to the RMCGC. The main difference with the percolation model introduced in Ref. [4], and the consequent definition of Mutually Connected Giant Component (MCGC), is that step (ii) must be substituted with "Every replica node that has at least a single replica node not in the giant component of its proper layer is removed from the network and considered as damaged, i.e., if a replica node is damaged all its interdependent replica nodes are damaged" [4,[20][21][22][23][26][27][28]. In particular, the RMCGC and the MCGC are the same for M = 2 layers, but they differ as long as the number of layers M > 2. In the latter case, the RMCGC naturally introduces the notion of redundancy or complementarity among interdependent nodes.
Note that the choice of considering just two operating layers, i.e. assumung that a replica node can remain functional as long as there is at least another interdependent replica node that is also functional, is a simplification. In a real scenario redundant interdependencies can include more than two operating layers or even a different number of operating layers for each node.
As a proof of concept to demonstrate the difference between the notions of MCGC and RMCGC, in Figure 2 we present results of numerical simulations for the two percolation models applied to the air transportation network within the US [39]. We remark that this analysis represents only a starting point to illustrate the key impact of the setting of redundant interdependency. Airports are the nodes in network. Two airports are connected if at least a flight is connecting them. Layers correspond to flights operated by the three major carriers in the US: American Airlines, Delta, and United. In this system, having a connected component that is shared by at least two carriers could be important for several reason. For example, competition on similar itineraries may favor a market for plane tickets that is more fair than the one that would be present in case of monopoly by a single carrier. If only two layers are considered, MCGC and RMCGC coincide and both metrics indicate the relative size of the system where no monopoly is present. The addition of a third layer in the system should have beneficial effects for the system by increasing the size of the system where monopoly is absent. This scenario is described by the redundant percolation model. Beneficial effects of the complementarity among routes offered by the various carriers are not only visible when the system is fully functional (i.e., parameter p = 1), but also when a relatively large fraction of airports are considered as removed from the system (approximately for p ≥ 0.5). The percolation model introduced in Ref. [4] instead describes a much more restrictive situation, where the largest cluster  [39]. We consider only US domestic flights operated in January 2014 by the three major carriers in the US (American Airlines, Delta, and United), and construct a multiplex network where airports are nodes, and connections on the layers are determined by the existence of at least one flight operated by a given carrier between the two locations. The number of nodes in the network is N = 183. Please note that some of the nodes appear as connected only in one layer, therefore the relative size of the largest cluster is always smaller than one. (a) We consider a single realization of the random damage by assigning to every replica node (i, α) a random number extracted uniformly in the interval [0, 1]. Replica nodes are considered as damaged if their associated random number is smaller than p. Note that the same exact configuration of random damage is considered in all cases. In the percolation diagram, we consider the size of the MCGC and RMCGC as a function of p. Large symbols are results of numerical simulations, whereas the tick line is obtained from the solution of our mathematical framework [system of Eqs. is formed by airports that are connected simultaneously by all three carriers. The considerations reported above are valid for results valid both for a single instance of the percolation model (Fig. 2a), and for average values over multiple instances of the model (Fig. 2b).
As it appears clear from Figure 2, our theory is able to reproduce with high accuracy the results of numerical simulations. The next sections are devoted to the description of a complete mathematical framework that allows us for the description of the redundant percolation model. We stress that the framework is devised for arbitrary topologies, and can be therefore applied safely also to real multilayer networks as long as their structures are sufficiently compatible with the locally treelike approximation.

III. MESSAGE-PASSING ALGORITHM
We assume that interactions within each layer α are described by elements a [α] ij of the adjacency matrix of the layer, indicating whether the replica nodes (i, α) and (j, α) are connected (a Additionally, we consider a specific realization of the initial damage to the replica nodes indicated by the set {s iα }. The generic element s iα = 0 indicates that the replica node (i, α) has been initially damaged, whereas s iα = 1 indicates that the replica node (i, α) has not been initially damaged. Under these conditions, as long as the multiplex network is locally treelike, the following message-passing algorithm identifies the replica nodes that are in the RMCGC.
Each node i sends to a neighbor j a set of messages n [α] i→j in every layer α where node i is connected to node j, i.e., with a (a) node i is connected to node j in layer α, and both nodes (i, α) and node (j, α) are not damaged, i.e., s iα = s jα = a [α] ij = 1; (b) node i is connected to the RMCGC through at least one node = j in layer α; (c) node i belongs to the RMCGC assuming that also node j belongs to the RMCGC. This conditions is satisfied if and only if, assuming that node j belongs to the RMCGC, node i is connected in at least two layers to the RMCGC.
If the previous conditions are not simultaneously met, then n [α] i→j = 0. Put together, the former conditions lead to Here N α (i) indicates the set of nodes that are neighbors of node i in layer α. The term 1 − ∈Nα(i)\j 1 − n [α] →i therefore will equal one if at least one message is arriving to node i from a neighboring node = j, while it will be equal to zero, otherwise. θ(v i→j , 2) = 1 for v i→j ≥ 2 and θ(v i→j , 2) = 0, otherwise. v i→j indicates in how many layers node i is connected to the RMCGC assuming that node j also belongs to the RMCGC, i.e., Therefore v i→j indicates the number of initially undamaged replica nodes (i, α) that either receive at least one positive messages from nodes ∈ N α (i) \ j or are connected to the undamaged replica nodes (j, α). Finally, the replica node (i, α) belongs to the RMCGC if (i) it is not damaged, (ii) it is connected to the RMCGC in layer α, and (iii) it receives at least another positive message in a layer α = α. These conditions are summarized by The average number S of replica nodes belonging to the RMCGC is computed as The system of Eqs. (1), (2), (3), and (4) represents a complete mathematical framework to estimate the average size of the RMCGC for a given network and a given initial configuration of damage. The solution can be obtained by first iterating Eqs. (1) and (2) to obtain the values of the messages n i→j . Those values are then plugged into Eqs. (3) to compute the values of the variables s iα , and finally these variables are used into Eq. (4) to estimate the average size of the RMCGC. We stress that, being valid for a given network and for a given configuration of damage, the values of the variables n i→j and s iα are either 0 or 1. The variables v i→j can assume instead integer values in the range [0, M ]. The mathematical framework works properly also in presence of edge overlap among layers. This is an important feature that can change dramatically change the robustness properties of multilayer networks [26,[37][38][39][40].

IV. MULTIPLEX NETWORKS WITHOUT LINK OVERLAP
A. General results

Simplification of the message-passing equations on a single realization of the initial damage
In absence of link overlap, a given pair of nodes i and j may be linked exclusively along a single layer α. Nontrivial messages potentially different from zero will therefore exist only on a specific layer for every pair of connected nodes i and j. It can be easily seen that the messagepassing Eqs. (1) and (2) reduce to We further notice that in this situation the result of the message-passing algorithm does not change if we consider messages that depend exclusively on the state s iα of the node i that sends the message. Even if we drop the factor s jα in Eq. (5), the message will be allowed anyways to propagate further at the next iteration step, if the replica node (j, α) is not initially damaged. Therefore, we can further simplify Eq. (5) and consider Eqs. (6) replace Eqs. (1) and (2) in the case of a multiplex network without link overlap. The rest of the framework is identical, so that Eqs. (3) and (4) remain unchanged.

Message-passing equations for random realizations of the initial damage
Eqs. (6), (3), and (4) determine the average size of the RMCGC in a multiplex network without link overlap for a given realization of the initial damage {s iα }. These equations can be, however, extended to make predictions in the case of a random realization of the initial damage when the replica nodes are damaged independently with probability 1 − p, i.e,. such that the the initial damage {s iα } is a random configuration obeying the probability distribution To this end, we denote the probability that node i sends a positive message to node j in layer α byn i→j , and the probability that the replica node (i, α) belongs to the RMCGC byσ iα = σ iα . The message-passing algorithm determining the values ofn [α] i→j andσ iα is given byn This algorithm can be applied to a given network, and provides the average number of replica nodes S belonging to the RMCGC for a random realization of the initial damage obeying Eq. (7). Specifically the value ofŜ is related toσ iα byŜ 3. Message-passing equations for random multiplex networks A multiplex network where every layer is a sparse network generated according to the configuration model is a major example of a multiplex network without link overlap in the limit of large network sizes. It is therefore natural and important to characterize the RMCGC in this case. We assume that every network layer G α is a random graph taken from the probability distribution where k indicates the pre-imposed degree of node i in layer α, δ(x, y) = 1 if x = y and δ(x, y) = 0, otherwise, and Z is the normalization factor indicating the total number of networks in the ensemble. Averaging over the network ensemble allows us to translate the messagepassing equations into simpler expressions for the characterization of the percolation transition.
Let us consider a random multiplex network obeying the probability of Eq. (10), and a random realization of the initial damage described by the probability of Eq. (7). The average message in layer α, namely ij = 1 , and the average number of replica nodes of layer α that are in the RMCGC, denoted by S α = σ i,α , obey the equations where P (k) indicates the probability that a generic node i has degrees k i = k, i.e. (k [1] i , k [2] i , . . . , k If there are no correlations between the degrees of a node in different layers, the degree distribution P (k) can be factorized as where P [α] (k) is the degree distribution in layer α. In this case, Eqs. (11) can be expressed in terms of the generating function of the degree distribution in each layer. Specifically, we have where the generating functions H Finally the average number S of replica nodes in the RMCGC is given by If we consider the case of equally distributed Poisson layers with average degree z, we have that Eq. (12) is for every layer α = 1, 2, . . . , M . Then, using Eqs. (13), one can show that S α = S α = S, ∀α, and S is determined by the equation This equation has always the trivial solution S = 0. In addition, a nontrivial solution S > 0 indicating the presence of the RMCGC, emerges at a hybrid discontinuous transition characterized by a square root singularity, on a line of points p = p c (z), determined by the equations where For p > p c there is a RMCGC, for p ≤ p c there is no RMCGC. The entity of the discontinuous jump at p = p c in the fraction S of replica nodes in the RMCGC is given by S = S c . The percolation threshold p c as a function of the average degree z of the network is plotted in Figure 3 for M = 2, 3, 4, 5. It is shown that as the number of layers M increases the percolation threshold decreases for every value of the average degree z. Also, the discontinuous jump S c decreases as the number M of layer increases for very given average degree z (see Figure 4). Therefore, as the number of layers increases the multilayer networks becomes more robust.

B. Comparison between the RMCGC and the MCGC
In this section, we compare the robustness of multiplex networks in presence of ordinary interdependencies and in presence of redundant interdependencies. To take a concrete example, we consider the case of a multiplex network with M Poisson layers, each layer having the same average degree z. In this case the fraction S of replica nodes in the RMCGC is given by the solution of Eqs. (17) while the fraction of replica nodes in the MCGC is given by [45] S =p 1 − e −zS M .
In Eq. (20), it is assumed that every replica node (i, α) of a given node i is damaged simultaneously (with probabilityf = 1 −p). On the contrary, in presence of redundant interdependencies it is natural to assume that the initial damage is inflicted to each replica node independently (with probability f = 1 − p). Therefore, in order to compare the robustness of the multiplex networks in presence and in absence of redundant interdependencies, we set p =p = 1, i.e., replica nodes are not initially damaged, and compare the critical value of the average degree z = z at which the percolation transition occurs respectively for the RMCGC and for the MCGC. Additionally we will characterize also the size S = S of the jump in the size of the RMCGC and the MCGC at the percolation transition. In Fig. 5, we display the values of z and S as a function of the number of layers M for the RM-CGC and the MCGC. For M = 2, the two models give the same results as they are identical. For M > 2, differences arise. In presence of redundant interdependencies, multiplex networks become increasingly more robust as the number M of layers increases. This phenomenon is apparent from the fact that the RMCGC emerges for multiplex networks with an average degree of their layers z which decreases as the number of layers M increases. On the contrary, in ordinary percolation the value of z for the emergence of the MCGC is an increasing function of M . Additionally, the size of the discontinuous jumps S at the transition point decreases with M for the RM-CGC, while increases with M for the MCGC showing that the avalanches of failures have a reduced size for the RMCGC. We expect that the beneficial effect of the addition of new layers will extend also to the scenario in which the number of operating layers necessary for a node to be functional is assumed to be lager than two. Also, we believe that the same conclusion will apply to more general multilayer network structures where nodes have different number of interlinks (redundant interdependencies) such as a the topologies considered in Ref. [28] as long as the number of operating layers necessary for a node to be functional is the same for every node.

C. Comparison with numerical simulations
In this section, we compare the results obtained with Eqs. (1), (2), (3), and (4) on a single instance of damage with the predictions the message-passing algorithm described in Eq. (13) characterizing the size S of the RMCGC in an ensemble of networks. Specifically, we consider the case of a multilayer network with M = 3 Poisson layers with the same average degree z. In order to draw the percolation diagram for single instances of initial damage as a function of the probability of damage 1 − p, we associate each replica node (i, α) with a ran-dom variable r iα drawn from a uniform distribution and we set In isolated networks, two nodes can be either connected or not connected. In multiplex networks instead, the complexity of the structure greatly increases as the ways in which a generic pair of nodes can be connected is given by 2 M possibilities. A very convenient way of accounting for all the possibilities with a compact notation is to use the notion of multilink among pairs of nodes [7,16]. Multilinks m = m [1] , m [2] , . . . , m [M ] with m [α] = 0, 1, describe any of the possible patterns of connections between pairs of nodes in a multilayer network with M layers. Specifically, m [α] = 1 indicates that a connection exists in layer α, whereas m [α] = 0 indicates that the connection in layer α does not exists. In particular, we can say that, in a multiplex network with M layers, two nodes i and j are connected by the multilink m ij = (a [1] ij , a [2] ij , . . . , a In order to distinguish the case in which two nodes are not connected in any layer with the case in which in at least one layer the nodes are connected, we distinguish between the trivial multilink m = 0 and the nontrivial multilinks m = 0. The trivial multilink m = 0 indicates the absence of any sort of link between the two nodes. Using the concept of multilinks, one can define multiadjacency matrices A m whose element A m ij indicates whether (A m ij = 1) or not (A m ij = 0) a node i is connected to node j by a multilink m. The matrix elements A m ij of the multiadjacency matrix A m are given by We note that multiadjacency matrices are essentially equivalent to rank-3 tensors as those considered in Ref. [15] for multiplex networks, and generalized to the case of arbitrary multilayer networks in Ref. [14]. Using multiadjacency matrices, it is straightforward to define multidegrees [7,16]. The multidegree of node i indicated as k m i is the sum of rows (or columns) of the multiadjacency matrix A m , i.e., and indicates how many multilinks m are incident to node i. Using a multidegree sequence {k m i }, it is possible to build multiplex network ensembles that generalize the configuration model. This way, overlap of links is fully preserved by the randomization of the multiplex network. These ensembles are specified by the probabilityP( G) attributed to every multiplex network G of the ensembles, whereP( G) is given bỹ withZ normalization constant equal to the number of multiplex networks with given multidegree sequence.
B. General discussion of the message-passing equations for the RMCGC Our goal here is to generalize the message-passing algorithm already given by Eqs. (1), (2), (3), and (4) for a generic single instance of a multiplex network and single realization of initial damage to the cases of (i) random multiplex networks with given multidegree sequence and/or (ii) random realizations of the initial damage. The extensions for both cases has been already considered for the case of multiplex networks without link overlap. In presence of link overlap, however, a more complex formalism is needed. For two nodes i and j in fact, the messages n i→j given by Eq. (1) and sent from node i to node j over the different layers α = 1, 2 . . . M are correlated because they all depend on the value of the variable v i→j given by Eq. (2). Such correlations require particular care when averaging the messages to treat the percolation transition for random initial damages. Similar technical challenges are also present in the treatment of the MCGC model where interdependencies are not redundant [26,37]. In presence of redundant interdependencies there is an additional precaution that needs to be taken. In fact, the messages n [α] i→j are explicitly dependent on the state of all replicas (j, α ) of node j. This state is indicated by the variables s j = (s j1 , s j2 , . . . , s jα , . . . s jM ) where s jα specifies whether the replica node (j, α ) is initially damaged or not. As a consequence of this property, when averaging over random realizations of initial damage, message-passing equations are written in terms of the messagesσ mij , n i→j ( s j ) explicitly accounting for the probability that node i is sending to node j the set of messages n = (n [1] i→j , n [2] i→j . . . n i→j ), given that node j is in state s j and node i and node j are connected by a multilink m = m ij . We have derived these equations for a general multiplex network with M layers. However, the message-passing algorithm has a very long expression. To make the paper more readable, we decided to place the exact treatment of the general case in the SM, and consider here only the special case of ensembles of random multilayer networks with overlap. For these ensembles in fact, the message-passing equations are written in terms of average messages sent between nodes with given multilinks m, i.e., S m, n ( s j ) = σ mij , n i→j ( s j )| m ij = m , and the equations greatly simplify. Two specific cases of multilayer network ensembles are discussed below, for the cases of M = 2 and M = 3 layers.

C. Ensembles of multilayer networks link overlap and M = 2 layers
In this case, every replica node is in the RMCGC if and only if also its interdependent node in the other layer is in the RMCGC. Therefore, the only messages that are different from zero are the messages S m, n ( s j = (1, 1)) sent to nodes j in state s j = (1, 1). Specifically, we consider the case of a random multiplex network with Poisson multidegree distributions characterized by the averages The messages S m, n ( s j = (1, 1)) only depend on the multiplicity of overlap of the multilinks m and n given respectively by The fraction S of replica nodes in the RMCGC is determined by the variables x µ,ν = S m, n ( s j = (1, 1)), where for example the value of x 2,2 indicates the probability that node i to sends a message n = (1, 1) to its neighbor j with s j = (1, 1) connected by a multilink m = (1, 1). The values of the variables x µ,ν and S are determined by the following set of equations (See Supplementary Material for details) These equations are the same equations as those that determine the value of the MCGC as long we make the substitution p 2 → p [26,37] taking into account that the damage in each replica node is independent in the present model. Notably in this case the percolation phase transition is discontinuous and hybrid been characterized by a square-root singularity above for p approaching the percolation threshold p c from above.

D. Ensembles of multilayer networks link overlap and M = 3 layers
We consider now the case of a random multiplex network with M = 3 layers. The network has Poisson multidegree distributions and averages given by In this case, the messages S m, n ( s j ) only depend on Therefore the fraction of replica nodes in the RMCGC S is determined by the variables The equations that these variables need to satisfy can be described in a symbolic way by suitable diagrams (see the Supplementary Material for details in how to read these diagrams). Diagrams that describe the equations to determine the value of all the variables x (ξ) µ,ν are presented in Fig. 8. These equations read as We note that, in absence of overlap, i.e., for z 1 = z, z 2 = 0 and z 3 = 0, Eqs. (SM38) reduce Eqs. (17). By defining a suitable order of the variables x (ξ) µ,ν , it is possible to introduce a vector x whose elements are the variables x (ξ) µ,ν , and rewrite the Eqs. (SM38) in a matrix form as The hybrid discontinuous phase transition (characterized by a square root singularity) can be found by imposing that the system of Eqs. (35) is satisfied together with the condition that the determinant of the Jacobian J of G(x) equals to zero, that is Simulation results of the percolation process for multiplex networks in this ensemble are presented in Fig.9, and they provide clear evidence of a perfect agreement to the theoretical prediction.

VI. CONCLUSIONS
In this paper, we introduced and fully characterized an alternative percolation model for multiplex networks. The model serves to quantify the robustness of networks with redundant interdependencies. According to the model, interdependencies make a system more fragile than it would be by considering each layer independently. This fact is consistent with the original model used to study percolation in multiplex networks [4,21,23], and it is apparent from the fact that the transition is abrupt for any number of network layers considered in the interdependent model. On the other hand, redundancy of interdependencies across multiple layers favors system robustness, as the height of the discontinuous jump and the location of the transition point decrease as the number of layers increases. This is a fundamental difference with respect to the model currently adopted to study the robustness of multiplex networks and multilayer network in general, where instead increasing the number of layers generates more and more fragile networks [23,[26][27][28].
To characterize the model, we deployed a comprehensive theoretical approach based on message passing. Our theory is valid for arbitrary multiplex network topologies as long as they are locally treelike. The theory is further developed in the context of ensembles of synthetic network models to analyze properties of the percolation transition emerging in the new model, and to perform systematic comparisons, and emphasize fundamental differences, with the percolation model introduced in Ref. [4].
We remark that ours is not a definitive model, but it represents only a good starting point towards a more real-istic description of real multiplex systems. For instance, we do not expect to observe a qualitative change of the results if the redundancy is weaker and let us say a node is damaged only if less than two interdependent nodes are functioning. However, our approach is purely structural and in many situations, forgetting about intrinsic dynamical features may be not appropriate. The primary role of our model is to emphasize the importance of redundancy or complementarity in multilayer networks, an obvious-yet-neglected feature of many real systems. In several realistic settings in fact, system robustness is augmented by the addition of new layers of interactions, as these layers are indeed created to provide backup options. For example, adding a new mode of transportation in a preexisting multimodal transportation system should make the system more reliable against eventual failures. Similarly in a living organism, the development of new types of interactions among constituents should increase the stability of the same organism against possible mutations. In the current setting, the model assumes that the functioning of individual nodes requires that nodes are correctly operating on at least two interdependent layers. The model can be, however, generalized to deal with a variable number of minimal functioning layers to describe more realistic scenarios in specific situations of interest.

Message passing for given multilayer network and given initial damage
Let us consider a given multilayer network G with M layers. Each layer α = 1, 2, . . . , M of the multilayer network has adjacency matrix a [α] . In this multilayer network, each pair of nodes i and j is connected by a multilink m ij = (a [1] ij , a [2] ij . . . , a Any two nodes i and j are connected by a nontrivial multilink is m ij = 0 implying that at least one link between the two nodes is present across the M layers. We assume that the initial damage configuration is known and that it is given by the set of variables {s iα } where s iα indicates if a replica (i, α) is initially damaged (s iα = 1) or not (s iα = 0). The message-passing algorithm given in Sec. III of the main text allows us to determine for any given initial damage configuration, if any replica node (i, α) is in the RMCGC (σ iα = 1) or not (σ iα = 0) as long as the multilayer network is locally treelike. Specifically the variables σ iα are determined in terms the set of messages n i→j = (n [1] i→j , n [2] i→j , . . . , n going from any node i to any node j joined by a nontrivial multilink m ij = 0. The messages n i→j are determined according to the following recursive equation where N α (i) indicates the set of nodes that are neighbor of node i in layer α and where θ(x) is the step function with values θ(v i→j , 2) = 1 for v i→j ≥ 2 and θ(v i→j , 2) = 0 for v i→j = 0, 1. Here the variable v i→j indicates in how many layers node i is connected to the RMCGC assuming that node j also belongs to the RMCGC, Finally the variables σ iα are expressed in terms of the messages n i→j and are given by In many situations, however, the initial configuration of the damaged {s iα } is not known, and instead it is only known the probability distributionP({s iα }) of the initial damage configuration.
In this case, one aims to know the probabilityσ iα = σ iα that a replica node (i, α) is in the RMCGC for a random configuration of the initial damage. The value ofσ iα , on a locally treelike multilayer network is determined by a distinct message-passing algorithm that can be derived from the message-passing algorithm valid for single realization of the initial damage, by performing a suitable average of the messages. Particular care should be taken when one aims to perform this average. In fact σ iα depends on all the messages n i→j sent by node i to node j in all the layers α. These messages are correlated and therefore they cannot be averaged independently.
An alternative formulation of the Eqs. (SM 3) − (SM 11) provides the necessary framework for deriving in few steps the message-passing algorithm to predictσ iα . This alternative formulation is written terms of the variables σ m, n i→j indicating whether (σ m, n i→j = 1) or not (σ m, n i→j = 0) node i send to node j a message n = n i→j given that node j is connected to node i by a multilink m = m ij . Using where n i→j is determined in terms of the messages σ m, n i→j as Finally a replica node (i, α) is in the RMCGC (σ iα = 1) or not (σ iα = 0) depending on the messages it receives from its neighbors, i.e. (SM10)

Message-passing algorithm for random damage
By averaging Eqs. (SM 6)−(SM 7)−(SM 10) we can derive the message-passing algorithm predicting the probabilitŷ σ iα that a replica node (i, α) is in the RMCGC when the initial damage {s i } is randomly drawn for the probability distributionP({s iα }). Assuming that each replica node is damaged independently the probability distribution P({s iα }) is given byP The message-passing algorithm valid for a random distribution of the initial disorder, is written in terms of the messagesσ m, n i→j ( s). The messagesσ m, n i→j ( s) take real values between zero and one. They indicate the probability that node i send to node j a message n = n i→j given that node j is connected to node i by a multilink m = m ij and that node j has initial damage configuration s = s j , i.e. (s 1 , s 2 , . . . , s α , . . . , s M ) = (s j1 , s j2 , . . . , s jM ).
Let us indicate withP ( s) the probability of a local initial damage configuration given bŷ and let us indicate with r the vector r = (r [1] , r [2] , . . . , r ( s i ) where Finally the probabilityσ iα that a replica node (i, α) is in the RMCGC is given bŷ

Average over multilayer ensemble with give multidegree sequence
In order to derive the phase diagram of the percolation transition in presence of redundant interdependencies over given multilayer network ensembles, it is useful to consider a further average of the messagesσ m, n i→j . To this end we consider the multilayer network ensemble that preserves the multidegree sequence {k m i }. Every multilayer network G in this ensemble has probabilityP whereZ is a normalization constant equal to the number of multilayer networks with the given multidegree sequence.
In this multilayer network ensemble the average messages S m, n ( s) = σ mij , n i→j | m = m ij indicate the probability that a message n is sent toward a node with initial damage configuration s over a multilink m. These average messages can be found by solving the following recursive equations: Finally the probability S α that a replica node in layer α is in the RMCGC in the multilayer network ensemble is given by 4. Derivation of Eq. (SM 14) In this section, we will discuss in detail the derivation of Eq. (SM 14) from Eq.(SM 6). A similar derivation (that we omit here) can be performed to derive Eqs. (SM 16)/(SM 18) from Eqs. (SM 7)/(SM 10).
We start from Eq. (SM6) written for the messages σ m, n i→j sent from a node i to a node j with n satisfying ν = M α=1 n [α] > 1 and m = (a [1] ij , a [2] ij . . . , a [M ] ij ). This equation is rewritten here for convenience, We given a set of variables p [α] = 0, 1 we can use the following identity where in the last expression we perform a sum over all the M -dimensional vectors r r = (r [1] , r [2] , . C m, n ( s i , s, r) where C m, n ( s i , s, r) is given by Eq. (SM 15). By using the fact that the messages σ m, n i→j take only values zero or one, that that out of all the messages σ m, n i→j from node i to node j only one is actually equal to one, and all the others are zero, we can rewrite Eq. (SM 28) as Finally, averaging over the probability distributionP ( s i ) of the configuration s i of the initial damage of node i, in the locally treelike approximation we obtain for the messagesσ m, n i→j ( s j ) the Eq. (SM 14) that we rewrite here for convenience,σ m, n i→j ( s j ) = (SM30)

Ensembles of multilayer networks link overlap and M = 2 layers
In the multilayer network case with M = 2 and link overlap, every replica node is in the RMCGC if and only if also its interdependent node in the other layer is in the RMCGC. Therefore, the only messages that are different from zero are the messages S m, n ( s j = (1, 1)) sent to nodes j in state s j = (1, 1). Specifically, we consider the case of a random multilayer network with Poisson multidegree distributions characterized by the averages The messages S m, n ( s j = (1, 1)) only depend on the multiplicity of overlap of the multilinks m and n given respectively by The fraction S of replica nodes in the RMCGC is determined by the variables x µ,ν = S m, n ( s j = (1, 1)). (SM33) The value of x 2,2 indicates the probability that node i to sends a message n = (1, 1) to its neighbor j with s j = (1, 1) connected by a multilink m = (1, 1). This fact occurs if and only if node i has both replica nodes that are not initially damaged (which occurs with probability p 2 ) and if at least one positive message in each layer α is reaching node i from neighbors different from j. The value of x 1,1 indicates the probability that node i sends a message n = (1, 0) to its neighbor j with s j = (1, 1) connected by a multilink m = (1, 0) or equivalently sends a message n = (0, 1) to its neighbor j with s j = (1, 1) connected by a multilink m = (0, 1). This fact occurs if and only if node i has both replica nodes that are not initially damaged (which occurs with probability p 2 ) and if at least one positive message in each layer α is reaching node i from neighboring nodes different from j. The latter is a necessary condition to have v i→j = 2. The value x 2,1 indicates the probability that node i is sending a message n = (1, 0) to its neighbor j in state s j = (1, 1) and connected by a multilink m = (1, 1) or equivalently sends a message n = (0, 1) to its neighbor j in state s j = (1, 1) and connected by a multilink m = (1, 1). This fact occurs if only if node i has both replica nodes that are not initially damaged (which occurs with probability p 2 ) and if at least one positive message is reaching node i in the layer for which n [α] = 1 from neighbors different from j and no positive message is reaching node i in the layer where n [α] = 1 from neighboring nodes different from node j. Finally, S is the probability that a replica node (i, α) is in the RMCGC which implies that (i) it is not initially damaged, (ii) its replica node in the other layer is not initially damaged, and (iii) at least one positive message reaches node i in both layers. The values of the variables x µ,ν and S are therefore determined by the following set of equations These equations are the same equations as those that determine the value of the MCGC as long as the fact that the damage in each replica node is independent is taken into account, which can be done by making the substitution p 2 → p [26,37].
The fraction of replica nodes in the RMCGC S is determined by the variables x (ξ) µ,ν = S m, n ( s j ).
Let us explicitly describe the equations that one of these variables needs to satisfy, and introduce a symbolic way to describe the equations. Specifically, we consider x 3,2 as the probability S (1,1,1),(1,1,0) ((1, 1, 1)) that a node i, connected 3,2 = S (1,1,1), (1,1,0) [(1, 1, 1)] in a multilayer network with M = 3 layers and k (1,1,1) = z3, k (1,1,0) = k (1,0,1) = k (0,1,1) = z2 and k (1,0,0) = k (0,1,0) = k (0,0,1) = z1. Filled circles indicate initially undamaged replica nodes siα = 1, whereas empty circles indicate initially damaged replica replica nodes siα = 0. The message are sent along the direction indicated by the arrows. A solid line reaching node i in layer α indicates that at least one positive message is reaching node i from nodes different from node j in layer α. Dotted lines joining node i in layer α indicate that no positive message reaches node i from nodes different from node j in layer α. A solid (dotted) line between node i and node j in layer α indicates m [α] = 1 (m [α] = 0). to a node j by a multilink m = (1, 1, 1), sends to node j a message n = (1, 1, 0) provided that node j is in the state s j = (1, 1, 1) (see Fig. SM1). This probability is equal to the sum of (i) the probability that node i is in the state s i = (1, 1, 0) [which occurs with probability (1 − p)p 2 ] and it sends the message n = (1, 1, 0) to node j and (ii) the probability that node i is in the state s i = (1, 1, 1) (which occur with probability p 3 ) and sends the same message to node j. Node i sends the message n = (1, 1, 0) only if the following conditions are met: (i) if node i is in the state s i = (1, 1, 0), node i must receive at least one positive message from nodes different from node j in layers α = 1 and α = 2.
(ii) if node i is in the state s i = (1, 1, 1), node i must receive at least one positive message from nodes different from node j in layers α = 1 and α = 2 and must not receive any positive message from nodes different from node j in layer α = 3.
In the main text we compared the result of these equation with simulation results of the percolation process in this ensemble of multiplex networks averaged over different multiplex network realizations. Here for completeness we also show simulation results of the percolation process on a single instance of multiplex network and given initial damage in this ensemble with the results of the message-passing algorithm proposed in Sec. III of the main text (see figure  SM 3).