Cascading Failures in Interdependent Networks with Multiple Supply-Demand Links and Functionality Thresholds

Various social, financial, biological and technological systems can be modeled by interdependent networks. It has been assumed that in order to remain functional, nodes in one network must receive the support from nodes belonging to different networks. So far these models have been limited to the case in which the failure propagates across networks only if the nodes lose all their supply nodes. In this paper we develop a more realistic model for two interdependent networks in which each node has its own supply threshold, i.e., they need the support of a minimum number of supply nodes to remain functional. In addition, we analyze different conditions of internal node failure due to disconnection from nodes within its own network. We show that several local internal failure conditions lead to similar nontrivial results. When there are no internal failures the model is equivalent to a bipartite system, which can be useful to model a financial market. We explore the rich behaviors of these models that include discontinuous and continuous phase transitions. Using the generating functions formalism, we analytically solve all the models in the limit of infinitely large networks and find an excellent agreement with the stochastic simulations.


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 functional 13 . 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: 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 ⁎ . We call k sX i , ⁎ 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 Scientific RepoRts | 7: 15059 | DOI:10.1038/s41598-017-14384-y 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 ⁎ = 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 = 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 state of the nodes varies according to their color: functional nodes ( ), nodes that do not fulfill the internal rule of functionality ( ) and nodes that fail due to the initial damage ( ). In addition, we have nodes that externally fail because they do not get enough supply from the other network ( );. In panel (a) we show the three internal rules of functionality for a node i (marked with the blue arrow) to be functional: (I) it must be connected to the GC (represented by the ∞ symbol), (II) it must belong to a component of size h which survives with probability 1 − q(h) (in this case k * ≡ k i * ), or (III) it must have a number of neighbors equal to or greater than k * ≡ k i * (we show k * = 2). In panel (b) we show the external rule of functionality for k s * = 1, and k s * = 2 in panel (c). In these cases P sA (k) = P sB (k) = δ k,3 , however, not all supplies are shown, nor are the internal connectivity links.
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 , ⁎ or more nodes to remain functional, similar to k-core percolation.
External functionality failure is similar to heterogeneous k-core percolation 37 . 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, 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 model 42 . 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 Thus the recursion relations for the stages n > 1 are 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 µ µ = = .
y g y y g y 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 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  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 model 11 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 transition 28 . Here q X (1) = 1 and   5 and P sA (k) = P sB (k) = δ k, 5 and system size N = 10 5 for different values of required supplies, ⁎ k sX = 1 (○), k sX ⁎ = 2 ( ), k sX ⁎ = 3 ( ), ⁎ k sX = 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 = 1 ( ). The symbols are the results of the stochastic simulations and the lines are the iterated values obtained by equations (3)(4)(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.
Scientific RepoRts | 7: 15059 | DOI:10.1038/s41598-017-14384-y 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.
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 . 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 = 2, the transition is continuous, but it is discontinuous for  . 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.
Scientific RepoRts | 7: 15059 | DOI:10.1038/s41598-017-14384-y k sX ⁎ > 2, and there is no transition for p > 0 if k sX ⁎ = 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 c A c 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 k sB 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 k sB sB , 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 percolation 43,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.
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,

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 networks 11,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 models 11,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 models 11,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 when γ → 0, and thus corresponding curves appear finite even for very small γ > 0. 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 = 10 6 to compute the steady state and N = 10 8 for the temporal evolution close to the critical threshold (See Fig. 2). We use the Molloy-Reed Algorithm 45 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.