Introduction

Studying complex systems includes analyzing how the different components of a given system interact with each other and how this interaction affects the system's global collective behavior. In recent years complex network research has been a powerful tool for examining these systems, and the initial research on isolated networks has yielded interesting results1,2,3.

A network is a graph composed of nodes that represent interacting individuals, companies, or elements of an infrastructure. Node interactions are represented by links or edges. Real-world systems rarely work in isolation and often crucially depend on one another4,5,6,7,8,9,10. Thus single-network models have been extended to more general models of interacting coupled networks, the study of which has greatly expanded our understanding of real-world complex systems. One intensive study of these "networks of networks" has focused on the propagation of failure among closely-related systems11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26. The great blackout of Italy in 2003 and the earthquake of Japan in 2011 were catastrophic events that demonstrated that breakdowns in power grids strongly impact other systems such as communication and transport networks, and that the failure of these networks in turn accelerates the failure of the power grid. The propagation of these "failure cascades" has received wide study in recent years11,12,13,14,15,16,17,18,19,20,21,27,28.

The simplest model of these systems consists of two interdependent networks in which nodes in one network are connected by a single bidirectional edge to nodes in a second network11. In this model a node is functional (i) if it belongs to the largest connected component (the "giant component") in its own network (the internal rule of functionality) and (ii) if its counterpart in the other network is also functional (the external rule of functionality). This original model has been extended to include localized and targeted attacks15,29,30,31,32 and mitigation13,25,26,33,34 and recovery strategies27,35. Recently it was found that the giant component membership requirement can be replaced by a weaker requirement of belonging to a cluster of a size larger than or equal to a threshold h * 28. Alternatively, a heterogeneous k-core condition can be applied as an internal functionality condition in which node i is functional when at least \({k}_{i}^{\ast }\) nodes among its k i immediate neighbors remain functional36,37,38,39,40. In this model the random failure of a critical fraction of nodes in an isolated network leads to an abrupt collapse of this network.

Although the original interdependent network model expanded our understanding of different coupled systems, the single-dependency relationship between nodes in different networks does not accurately represent what happens in real-world structures. A cascading failure model of a network of networks with multiple dependency edges has been applied to a scenario in which nodes fail only when they lose all their support nodes in the other network14,17, but nodes in complex real-world systems can be so fragile that the loss of a single support link can cause them to shut down. More generally, each node may require a certain minimal number of supply links connected to the nodes in the other network to remain functional. In the world-wide economic system, for example, banks and financial firms lend money to non-financial companies who must pay the amount back with interest after a stated period of time. If a single non-financial company becomes insolvent, the bank that lent money to this company will likely not fail, but if the number of companies that cannot pay back their loans is sufficiently large, the possibility of bank failure becomes real. This resembles the k-core process in a single network described above.

Here we model the process of cascading failure in a system of two interdependent networks A and B in which nodes have multiple connections or supply-demand links between networks. In the following, network X means either network A or B. Each node i in network X has k sX,i supply nodes in the other network that are connected to node i by supply links. This node remains functional at a certain stage of the cascade of failures if the number of its functional supply nodes in the other network remains greater or equal to its supply threshold k sX,i * ≤k sX,i . We call this the external functionality condition. We assume that a supply threshold is predefined for each node.

In principle, this model is non-trivial even if the survival of a node in network X does not directly depend on the internal connectivity of network X. In this case our model is equivalent to cascading failures in a bipartite network composed of two sets of nodes A and B connected only by supply-demand links, i.e., these networks only have external functionality. For generality, we add to the external functionality condition an internal functionality condition that can be one of the following: a node is functional (i) when it belongs to the giant component of its network ("giant component rule"), (ii) when it belongs to a finite component of size h that survives with probability 1 −q(h) ("mass rule"), and (iii) when a node i with internal connectivity k i has a number of functional neighbors greater than or equal to k i * ("k-core rule").

We develop a theoretical model that is solved using the formalism of generating functions. We present numerical solutions and compare them with stochastic simulations. We find that for all internal rules of functionality, increasing the \({k}_{sX}^{\ast }\) value increases system vulnerability and often causes a discontinuous transition. For the mass rule of internal functionality we find a continuous transition for some parameter values. We also study the asymptotic limit of a large number of supply links, and we find a relation between the critical threshold of initial failure and the ratio \({k}_{sX}^{\ast }\)/k sX .

Model

We assume that the system consists of two networks A and B with internal degree distributions P A (k) and P B (k), respectively, where k is the degree of a node within its own network. Each node i in network A is supplied by k sA,i supply links from nodes in network B, and each node j in network B has k sB,j demand links that act as supply links for nodes in network A. For simplicity we assume that the demand links in network A serve as supply links for nodes in network B, and that supply links in network A serve as demand links for nodes in network B. Thus each supply-demand link is a bidirectional link that connects a node in network A with a node in network B. If the internal degree of all nodes in networks A and B is zero, our model is equivalent to a bipartite network. We assume that the degree distribution of supply-demand links in network A is P sA (k) and the degree distribution of supply-demand links in network B is P sB (k). In principle, some nodes may not have supply links and still remain functional13. If this is the case, P sA (0) > 0.

The functionality of the nodes in both networks is related to their connections within their own network, which we call the internal rule of functionality. In addition, the state of the nodes also depends on the supply demand links that connect both networks, which we call the external rule of functionality.

We study three different internal rules of functionality:

  1. (I)

    Model I (The "giant component" rule): nodes that belong to the giant component in their own network are functional.

  2. (II)

    Model II (The "finite component" or "mass" rule): a finite component of size h remains functional with a probability 1 −q(h). If it fails, all of its nodes fail. If it survives, all of its nodes remain functional.

  3. (III)

    Model III (The "k-core" rule): a node i with internal connectivity k i remains active if the number of its functional neighbors is greater than or equal to \({k}_{i}^{\ast }\).

The external rule of functionality states that nodes in network X must be connected with the other network through a number of functional supply links greater than or equal to \({k}_{sX}^{\ast }\).

We call \({k}_{sX,i}^{\ast }\) the supply-demand functionality threshold of node i, since in principle the threshold may be different for different nodes. For conceptual simplicity, we assume that the supply thresholds are predefined for each node by random selection from a cumulative probability distribution r sX (j, k) =P(\({k}_{sX}^{\ast }\) ≤j|k sX  =k), where P(|) is the conditional probability. Alternatively, function r sX (j, k) can be understood as a probability that a node with k supply links remains functional if j of its k supply nodes in the other network remains functional.

For example, in the case of a uniform supply threshold \({k}_{sX}^{\ast }\) =m where m is a constant, the distribution r sX is a step function, i.e., r sX (j, k) = 0 for j <m and r sX (j, k) = 1 for j ≥m. Another option is linear: r sX (j, k) =j/k. For autonomous nodes that can survive without any functional supply nodes in the other network, \({k}_{sX}^{\ast }\) = 0. This case is included in the general scheme if we assume that r sX (0, k sX ) > 0.

Figure 1a shows a schematic of the internal rules of functionality, and Fig. 1b,c show a schematic of the external rules of functionality. In each network green nodes are functional, i.e., they satisfy both internal and external conditions of functionality. Red nodes are affected by the initial failure, blue nodes do not satisfy internal conditions of functionality and pink nodes do not satisfy external conditions of functionality. Internal links are black, and supply links are orange. Here we use P sA (k) =P sB (k) =δ k,3, but for simplicity in Fig. 1b,c we omit the internal links and some of the supply-demand links in network A. For example, in Fig. 1b node A3 has two additional supply nodes from network B that are not shown. Figure 1b shows the case k s,i * = 1 for all i. Note that since all nodes in network B receive supplies from functional node A0 they are unaffected when other nodes in network A fail. On the other hand, Fig. 1c shows that when k s,i * = 2 all nodes must have two functional supply nodes from the other network to remain functional. Nodes B2 and B3 are connected to A0, receive supplies from functioning nodes A3 and A6, respectively, and remain active. On the other hand, because node B1 is only supported by node A0, it fails, as indicated by the pink color.

Theoretical approach

We construct a system of two randomly connected networks in which connectivity links within each network follow degree distributions P A (k) and P B (k) and supply-demand links between the networks follow distributions P sA (k) and P sB (k). For this system we achieve a theoretical solution within the limit of a large number of nodes, N A and N B , where N A and N B are the number of nodes in networks A and B, respectively. The bidirectionality of the supply-demand links requires that relation N A 〈k〉 sA  =N B 〈k〉 sB is satisfied, where 〈k〉 sA and 〈k〉 sB are the average degrees of the supply links in networks A and B respectively.

When we randomly remove a fraction 1 −y X of nodes from network X, the remaining fraction of active nodes μ X for an isolated network X is determined by which internal functionality rule is followed. It can be expressed in the closed-form expression μ X  =y X g X (y X ), where g X (y X ) ≤ 1 is an exacerbation factor that takes into account additional node failures triggered by the random removal of a fraction of 1 −y X nodes. The explicit form of this factor is determined by the internal functionality rules of the model. The Supplementary Information presents equations for g X for Rules I, II, and III (see Supplementary Information: section Explicit form of the functionality rules). For example, for a bipartite network g X (y X ) = 1.

The cascading process begins with a random failure in network A. This failure causes an additional loss of nodes determined by the exacerbation factor. This event triggers a cascade in which failure is transmitted back and forth between networks A and B through the supply-demand links, and this further decreases the fraction of functional nodes. The external functionality rule states that node i with k s,i supply-demand links must have \({k}_{s,i}^{\ast }\) or more nodes to remain functional, similar to k-core percolation.

External functionality failure is similar to heterogeneous k-core percolation37. To describe this failure due to a lack of supply between networks A and B, we introduce the functions W sA (x),W sB (x) and Z sA (x),Z sB (x), which are the k-core generating functions of the degree distribution and the excess degree distribution of the supply-demand links in networks A and B, respectively. These functions depend on the degree distributions P sA and P sB of supply-demand links and the distribution of the thresholds r sA (j,k) and r sB (j,k) of the supply-demand links in networks A and B,

$${W}_{sX}(\beta )=\sum _{k=0}^{\infty }{P}_{sX}(k)\sum _{j=0}^{k}(\begin{array}{c}k\\ j\end{array}){r}_{sX}(j,k){\beta }^{j}{\mathrm{(1}-\beta )}^{k-j}$$

(1)

and

$${Z}_{sX}(\beta )=\sum _{k=0}^{\infty }\frac{k{P}_{sX}(k)}{{\langle {k}_{s}\rangle }_{X}}\sum _{j=0}^{k-1}(\begin{array}{c}k-1\\ j\end{array}){r}_{sX}(j+\mathrm{1,}\,k){\beta }^{j}{\mathrm{(1}-\beta )}^{k-j-1},$$

(2)

where 〈k s 〉 X is the average number of supply links per node in network X. In this context β is the probability that a functional node will be selected. Similar formulas were derived in ref.41 for a variant of the Watts opinion model42.

We next examine a theoretical approach to the temporal evolution of the cascading process. As explained above, initially a randomly selected fraction 1 −p of nodes fails in network A. Then the surviving fraction of nodes in network A in this first stage of the cascade is μ A,1 =pg A (p). We introduce a new parameter f B , which is the probability of randomly choosing a supply link that is connected to a functional node in the other network. When a node fails, all its demand links also fail. Thus f B,1 =μ A,1.

After applying the external functionality rule to network B, the fraction of nodes that fulfill the conditions is given by y B,1 =W sB (f B,1). Because there are additional disconnected nodes in network B given by the exacerbation factor g B , the number of functional nodes in network B at the first stage of the cascade is μ B,1 =y B,1 g B (y B,1). In the second stage of the cascade we cannot apply the same rules to obtain μ A,2, because f A,2 ≠ μ B,1. If, for example, a supply-demand link connects nodes i in A and j in B, then the probability that this link is active depends on how many other links belonging to nodes i or j are active. Thus the fraction of surviving links at this step is f A,2 =Z sB (f B,1)g B (y B,1).

Thus the recursion relations for the stages n > 1 are

$$\begin{array}{rcl}{f}_{A,n} & = & {Z}_{sB}({f}_{B,n-1})\,{g}_{B}({y}_{B,n-1});\\ {f}_{B,n} & = & p\,{Z}_{sA}({f}_{A,n})\,{g}_{A}({y}_{A,n}),\end{array}$$

(3)

where

$$\begin{array}{rcl}{y}_{A,n} & = & p\,{W}_{sA}({f}_{A,n});\\ {y}_{B,n} & = & {W}_{sB}({f}_{B,n})\end{array}$$

(4)

are the fractions of nodes that satisfy the external rule of functionality, i.e., randomly removing a fraction of 1 −y X,n nodes leaves the same number of functional nodes as in stage n of the cascade. The fractions of functional nodes at stage n of the cascade are

$$\begin{array}{rcl}{\mu }_{A,n} & = & {y}_{A,n}\,{g}_{A}({y}_{A,n});\\ {\mu }_{B,n} & = & {y}_{B,n}\,{g}_{B}({y}_{B,n}\mathrm{).}\end{array}$$

(5)

The process begins with f A,1 = 1 and y A,1 =p, which is equivalent to an initial random failure on network A.

Results

We next present these theoretical results using several simple examples and verifying them with stochastic simulations.

To test the validity of the equations, Fig. 2 shows the temporal evolution of the order parameter of networks A and B close to the critical threshold p c , computed using the equations and stochastic simulations when the giant component functionality rule is applied (see Supplementary Information: subsections Giant Component and Numerical Solution for the threshold p c ). Note that the plots show the simulation results are in total agreement with the theoretical results.

Figure 2
figure 2

Temporal evolution, close to the critical threshold, of the giant component μ A (n) and μ B (n) of networks A and B, when both are random regular (RR) networks with delta degree distribution P X (k) =δ k,5, with X =A,B. The degree distributions of supply links are also delta-distributions with P sA (k) =P sB (k) =δ k,5 and k s * = 2. The critical threshold for this system is p c  = 0.381. (a) p = 0.38, (b) p = 0.381. Network A (○, ), Network B (, ). The dashed lines are the results from the equations and the symbols are the results from the stochastic simulations.

Full size image

Figure 3 shows a plot of μ A and μ B in the steady state as a function of the initial fraction of surviving nodes p when the giant component rule is applied. The results for the k-core rule are shown in the Supplementary Information. We use two random regular (RR) networks with a degree distribution P X (k) =δ k,5, with X =A, B, and where the distribution of supplies is also RR with P s,A (k) =P s,B (k) =δ k,5. For the external rule of functionality we use r sX (j, k) = 0 if j <m and r sX (j, k) = 1 if j ≥m for all m from m = 1 to m = 4. The results obtained from the equations (dashed lines) agree with the results of the simulations (symbols). In addition we compare the results of the present model with the results of the original model of cascading failures11 shown as a dashed-dotted line in which P X (k) =δ k,5, but P s,A (k) =P s,B (k) =δ k,1 and m = 1.

Figure 3
figure 3

Two random regular (RR) networks with P A (k) =P B (k) =δ k,5 and P sA (k) =P sB (k) =δ k,5 and system size N = 105 for different values of required supplies, \({k}_{sX}^{\ast }\) = 1 (○), \({k}_{sX}^{\ast }\) = 2 (), \({k}_{sX}^{\ast }\) = 3 (), \({k}_{sX}^{\ast }\) = 4 (), as a function of the initial fraction of survived nodes p. Also shown two RR networks with P A (k) =P B (k) =δ k,5, but P sA (k) =P sB (k) =δ k,1, \({k}_{sX}^{\ast }\) = 1 (). The symbols are the results of the stochastic simulations and the lines are the iterated values obtained by equations (3–5). The dashed-dotted lines represent only the theoretical results since they have been obtained in Ref.11. In panels (a) and (b) we show the order parameter of network A and B, μ A and μ B , respectively for the giant component rule.

Full size image

Note that in network A the order parameter for all values of k s * is proportional to p until it begins to drop and become close to the critical threshold p c . This means that the depletion of the supply from network B does not significantly impact network A until it reaches the collapse threshold at which the system breaks down with a discontinuous transition. We calculate this critical value numerically using the generating functions (see Supplementary Information: section Numerical solution for the threshold p c ). Note also that, as expected, the behavior of network B is different. Because there is no initial random failure in network B, it remains more intact than network A. When network A crumbles, however, both networks collapse. Thus despite its damage being minor the transition in network B is more abrupt, more unexpected, and, therefore, more dangerous. This is the key difference between the present mode and the original model11 in which the behaviors of network A and B are identical. In addition, note that the system is more resilient when k s * is smaller, i.e., when the supply level decreases. We also observe that the interdependent system with only one supply-demand link (the dashed-dotted line) is more resilient than a system with more connections between the two networks, but with large functionality thresholds m ≥ 3.

If instead of the giant component we apply the k-core as an internal functionality rule we get the same qualitative results. For different values of k * and k s * the order parameters also undergo a discontinuous transition, and the system becomes more vulnerable when the threshold of internal links and the threshold of supply links increases (see Supplementary Information: section k-core Percolation).

When applying the "mass" rule, finite components of size h in network X survive with a probability 1 −q X (h). When all nodes have a single supply-demand link, i.e., when k s  = 1 and k s * = 1, and all finite components of size greater than or equal to h = 2 are preserved, the system undergoes a continuous transition28. Here q X (1) = 1 and q X (h) = 0 for h ≥ 2. If the number of supply links increases and the threshold k s * = 1 is fixed, the system becomes more resilient and the transition remains continuous. In contrast, if all the components of size h = 2 are removed [q X (2) = 1] the transition becomes discontinuous irrespective of the number of supply-demand links connecting the networks. Nevertheless, not all the components of size h = 2 need to survive to have a continuous transition. Figure 4 shows the order parameters for q(2) = 0.3 and q(2) = 0.85 when q A (h) =q B (h) =q(h). Note that when q(2) = 0.3 the transition is continuous even when some of the components of size h = 2 are deleted. When q(2) = 0.85 the number of surviving h = 2 components is insufficient to prevent an abrupt transition.

Figure 4
figure 4

Order parameters for the "mass rule", for a system of networks with internal distribution P A (k) =P B (k) =δ k,5, supply distributions P sA (k) =P sB (k) =δ k,2 and thresholds k A * =k B * = 1. All the components of size h = 1 are deleted (q(1) = 1), and all the components of size h ≥ 3 are preserved (q(3) =q(4) = ... =q(h max ) = 0 where h max is the maximum value of h). The curves represent the case q(2) = 0.3 (●, ), for which there is a continuous transition, and q(2) = 0.85 (, ), which leads to an abrupt breakdown of the order parameter. The dashed lines represent the theoretical results and the symbols the stochastic simulations. (a) Network A, (b) Network B.

Full size image

Thus when k s * = 1 there is a critical value of q(2) =q c (2) that separates the zone of continuous transition from the zone of discontinuous transition. Figure 5 shows a phase diagram for a system of networks following the "mass" rule with an internal distribution P A (k) =P B (k) =δ k,5 and supply distribution P sA (k) =P sB (k) =δ k,ks . Note that the behavior of the critical probability as a function of the number of supply-links k s between the networks delimits these two zones. As k s increases the system becomes more robust, and more components must fail to cause an abrupt transition. In the limiting case k s  → ∞ the curve reaches the value q c (2) = 1, but also p c  → 0. On the other hand, when k s * > 1 the transition is always discontinuous for any value of q(s) and sufficiently large k s .

Figure 5
figure 5

Phase diagram that shows the continuous and discontinuous transitions zones when the "mass rule" is applied. The curve represents the critical probability of failure of the components of size h = 2 as a function of the number of the supply-demand links. In this case P A (k) =P B (k) =δ k,5, P sA (k) =P sB (k) =δ k,ks and k sA * =k sB * = 1. For clarity, the k s axis is shown on a log scale.

Full size image

What happens if no internal functionality rule is applied? This could be the case in a bipartite system in which nodes within each network do not interact but use nodes in the other network as bridges to establish connections. Here the exacerbation factor is simply g X (y) = 1, which simplifies the equations. If we analyze this system for different functions r sX (j,k) (see Supplementary Information: section Examples of r sX (j,k) functions) we see that if r sX (j,k) is a step function with fixed threshold \({k}_{sX}^{\ast }\) = 2, the transition is continuous, but it is discontinuous for \({k}_{sX}^{\ast }\) > 2, and there is no transition for p > 0 if \({k}_{sX}^{\ast }\) = 1. Also if we choose a linear function, i.e., r sX (j, k sX ) =j/k sX , there is again no transition because here functions W s (β) and Z s (β) become linear functions of β. On the other hand, when the function r sX is by nonlinear, the behavior changes. Figure 6 shows the behavior of the order parameter of network A for a polynomial function r sX (j, k sX ) = 3(j/k sX )2 − 2(j/k sX )3 and for a supply-demand distribution P s,X (k) =δ k,ks . Note that for small values of k s the order parameter moves smoothly to zero but for k s  = 8 the system undergoes a discontinuous transition. The existence of these transitions can be explained studying Eqs (3) and (4) (see Supplementary Information: section Numerical solution for the threshold p c ).

Unlike the previous results, the transition here does not produce a total collapse of the system, and after the jump a small fraction of nodes remains functional for any p > 0. If a delta-distribution of supply links is replaced by the Poisson distribution with 〈k s 〉 X  =λ, we find a critical point on a (p, λ) plane λ c  = 7.58465, p c  = 0.728102 at which the first order phase transition emerges. For λ >λ c the transition is first order and for λ <λ c there is no phase transition for p > 0. At this point the system belongs to the mean-field universality class, such as the Ising model in infinite dimensions where p corresponds to the ordering field and λ to the thermal field.

We next analyze the limiting case of large k s values when all nodes in network B have a fixed threshold k sB *, and we find that the collapse threshold p c converges to a value determined by the ratio γ ≡k sB */k sB given by

$$\gamma ={p}_{c}{g}_{A}({p}_{c}),$$

(6)

which is valid for all of the internal failure rules.

The p c value depends on γ in this limit because when 〈k sX 〉 → ∞ the functions W sB (β) and Z sB (β) become step functions equal to 0 for β <γ and to 1, otherwise. Note that γ only relates to the external properties of network B, but that the value of p c depends solely on the topology of network A. This is because network B is intact above p c , but when p <p c all the supply-demand links maintaining the integrity of network B fail and the entire structure crumbles. Thus here the topology of network B does not affect the final state of the system. See Supplementary Information: section Asymptotic properties of the functions W s and Z s for the derivation of Eq. (6).

Figure 7 shows the behavior of Eq. (6) for each internal rule of functionality and for several values of internal connectivity z A in network A when it has an internal degree distribution P A (k) =δ k,zA . Note that all curves go to p c  = 1 when γ → 1, i.e., \({k}_{sB}^{\ast }\sim {k}_{sB}\), and thus even a small perturbation can cause a system breakdown. In contrast, curves with higher z A values have lower p c values because increased connectivity means increased resilience. In addition, when γ → 0 then \({k}_{sB}\gg {k}_{sB}^{\ast }\), rendering the influence of network B on network A insignificant. Here network A behaves as an isolated system. We see this in the giant component rule [see Fig. 7a] in which p c  → 1/(z A  − 1) as γ → 0, a value that corresponds to the critical threshold of node percolation43,44 in isolated RR networks. Similarly, for the "mass" rule we find that p c  → 0 when γ → 0 because when there is an initial attack 1 −p on an isolated network there are always components of varying masses in the thermodynamic limit (with an infinite number of nodes). Thus for any size h there are always surviving components when p > 0.

Figure 7
figure 7

Critical threshold p c as a function of γ =k sB */k sB for different values of z A , the internal connectivity of network A, where its internal degree distribution is RR. The curves represent different values of z A : z A  = 3 (─), z A  = 5 () and z A  = 10(). Panel (a) corresponds to the Giant Component rule. Panel (b) corresponds to the "mass rule", with q(h) = 1 for h = 1,2,3, and panel (c) to the k-core rule with k X * = 2 Note that in panel (b) p c  ~γ 1/4 when γ → 0, and thus corresponding curves appear finite even for very small γ > 0.

Full size image

If there is a Poisson internal degree distribution in network A, i.e., P A (k) =exp[−〈k〉 A ]〈k〉 A k /k! where 〈k〉 A is the mean connectivity, we can write a closed-form expression for p c for the giant component rule,

$${p}_{c}=\frac{\gamma }{1-{\exp }[-\gamma \,{\langle k\rangle }_{A}]}\mathrm{.}$$

(7)

Note that p c does not depend on the internal degree distribution of network B. The derivation of Eq. (7) is supplied in the Supplementary Information: section Asymptotic properties of the functions W s and Z s . On the other hand, if the system is bipartite then from Eq. (6) the critical value is simply p c  =γ.

Discussion

We have analyzed the cascading failure process in a system of two interdependent networks in which nodes within each network have multiple connections, or supply-demand links, with nodes from their counterpart network. In this model each node must have at least a given number of supply-links leading to functional nodes in the other network to remain active. We call this number the supply threshold and we call this condition the external functionality rule. We have studied the process under three internal functionality rules, (I) nodes must belong to the giant component in their own network, (II) nodes that belong to a finite component survive with a probability determined by the mass of the component, and (III) an internal version of the external functionality rule, known as heterogeneous k-core percolation. In addition, we have studied a system in the absence of any internal functionality rule, which is equivalent to a bipartite network. Our system is a generalization of the models of interdependent networks11,13 that represent a particular case of our model with P sX (k) = 0 for k > 1 and a giant component rule of internal functionality. Our model shows a rich behavior for various parameter values that is characterized by the appearance of discontinuous first order transitions. In some cases, multiple first order transitions can be observed, a situation impossible in the original models11,13.

We have found that for all the internal functionality rules the system is more robust when the supply threshold is lower. Under internal rules I and III there is a discontinuous transition at a collapse threshold p =p c . The main difference between our model and the previously studied models11,13 is that in the case of multiple supply links the initial attack on network A does not immediately affect network B, and it remains more functional than network A for any p >p c . This makes the transition, when it occurs in network B, more abrupt than in network A. These sudden breakdowns can come without warning. In some catastrophic events, e.g., an earthquake of sub-threshold strength, the damage to network B may be minor and the development of precautions or recovery strategies thus deemed of minor importance. This becomes problematic when the strength of an earthquake exceeds a certain threshold and causes a total breakdown in network B. In contrast, in "mass" rule II for k s * = 1 the transition can be continuous depending on the probability that components of size h = 2 remain functional and on the number of supply-demand links. For each value of k s there is a critical probability q(2) below which the transition becomes discontinuous.

When the model is applied to a bipartite system, the behavior is determined by function r sX . In particular, when this function is polynomial there is no transition in k sX  ≤ 7, but when k s increases this curve breaks and becomes discontinuous.

Finally we have studied the asymptotic limit value of the number of supply-demand links, and find that when r sB is a step function there is an exact relationship between the ratio γ = k sB */k sB and the collapse threshold p c . We also find that in this limit the resilience of the interacting system is enhanced up to the point at which the critical threshold p c is solely dependent on the topology of network A.

Methods

For the stochastic simulations we use for both networks a system size of N = 106 to compute the steady state and N = 108 for the temporal evolution close to the critical threshold (See Fig. 2). We use the Molloy-Reed Algorithm45 for the construction of the networks. The simulation results are averaged over 1000 network realizations.

For model II, the "mass" rule, a finite component of size h survives with probability 1 −q(h). In the stochastic simulations if a finite component remains after the internal failure at a step of the cascade, then in the following steps of the cascade this component only can fail due to the external rule of functionality.

In our theoretical analysis, to calculate the values of the order parameters at the steady state we iterate the temporal evolution Eqs (3)–(5) until the condition μ A  ≡μ A,n  =μ A,n − 1 is satisfied. At this stage the magnitudes of all order parameters reach a steady state and no longer change.