Interdependent resistor networks with process-based dependency

Studies of resilience of interdependent networks have focused on structural dependencies between pairs of nodes across networks but have not included the effects of dynamic processes taking place on the networks. Here we study the effect of dynamic process-based dependencies on a system of interdependent resistor networks. We describe a new class of dependency in which a node’s functionality is determined by whether or not it is actually carrying current and not just by its structural connectivity to a spanning component. This criterion determines its functionality within its own network as well as its ability to provide support-but not electrical current-to nodes in another network. We present the effects of this new type of dependency on the critical properties of σ and B ∞ ?> , the overall conductivity of the system and the fraction of nodes which carry current, respectively. Because the conductance of current has direct physical effects (e.g. heat, magnetic induction), the development of a theory of process-based dependency can lead to innovative technology. As an example, we describe how the theory presented here could be used to develop a new kind of highly sensitive thermal or gas sensor.


Introduction
Resistor networks have been a central topic of percolation research for decades. The seminal work of Kirkpatrick [1][2][3] and others demonstrated that significant information about general transport phenomena in disordered media could be deduced from d-dimensional lattices of resistors from which a fraction − p 1 of the sites or bonds are removed at random [4][5][6]. One of the important results which emerged from the study of resistor networks is the distinction between current carrying nodes, (the 'backbone') and nodes which are connected to the current carrying nodes but do not carry any current (the 'dead ends'), see figure 1. As shown in figure 1(c), in two-dimensions near criticality, the mass of dead ends completely overwhelms the backbone mass in the thermodynamic limit [7][8][9]. Though the value of the percolation threshold p c is the same for the conductivity and the largest connected component, their behavior at criticality is different, see figure 4(a) and reference [10]. As we will show below, this difference leads to differences in p c for conductivity and the spanning cluster in a system of coupled resistor networks, see figure 4(b).
When utilizing percolation theory of networks to describe system robustness, connection to a giant connected component (GCC) which is of order N is treated as a proxy for functionality [5,6,11,12]. As long as a path exists between a node and a finite fraction of other nodes, the node is considered functional. This is reasonable in many systems, like the internet [11,12], but may not be the case in the spread of information or disease on social networks, movement of people on road networks or, as we discuss here, in the case of general transport as exemplified by electric flow in resistor networks. In these cases, 'functionality' is more readily identified with 'activity' as determined by a specific process and not just with 'connectivity'.
In interdependent networks [13][14][15][16], node functionality depends on two factors: connection to the GCC in its network and support from a specific node in another network. The supporting node also requires connectivity and support to function. When the same two networks support and depend on each other, they are interdependent and failures are amplified as they are transferred back and forth between the two networks, creating cascading failures.
Here, we break with previous work on interdependent networks and consider functionality to be determined by a dynamic process and not just the potential for such a process to occur, as determined by the existence of links connecting the node to the GCC. We call this dependency functional or process-based to distinguish it from previous research on structural or topological dependency. In this work, we focus on current flow as the process which determines the process-based dependency. A node can function only if it has current flowing through it (it is part of the backbone) and there is current flowing through the node that it depends on in another network. Note however, that there is no electrical connection or current flow between the networks, only dependency links.
We present here theoretical predictions and extensive simulations of this system. We find that process-based dependency leads to significantly increased vulnerability and first and second-order transitions in the bulk conductivity at higher values of p c than structural dependency. The first-order phase transitions occur over a wider range of parameters compared to structural dependency. We find that the second-order transition has critical scaling in coupled networkswhich is indistinguishable from single networks, even when p c changes dramatically. After presenting our theoretical findings, we discuss possible laboratory realizations of this system.
(a) Model. We study a system composed of two interdependent two-dimensional square lattices with nearest neighbor links of unit resistance. We take open boundary conditions for the sides and connect the top and bottom of each lattice to uniform voltage and ground, respectively ( figure 1(a)). There is no electrical connection between the two networks and therefore the only interaction between the systems is via the dependency links. Each node in network A depends on exactly one node in network B and the dependency is mutual, in line with the 'no feedback condition' commonly assumed in interdependent networks [17]. This oneto-one correspondence means that the number of functional nodes in each network is identical. We note, though, that even if only a fraction of the nodes were interdependent (as in [18,19]), the size of the backbone and giant component would be essentially the same in both systems, due to symmetry. The dependency links are assigned uniformly at random with the restriction that they connect nodes with horizontal distance not longer Nodes (squares) are connected to their nearest neighbors with unit resistance. Nodes cease to function in three ways: (i) direct attack (black), (ii) structural separation from the spanning component (blue) or (iii) dynamic loss of electrical current (green). In our model, only the nodes in the currentcarrying backbone remain functional (white). As the system approaches criticality, the green nodes (dead ends) far outnumber the white nodes (backbone) [5,6]. than r in each direction, in line with the models used in [19][20][21], figure 1(a). This means that a node at location (i, j) in lattice A can depend on the continued functioning of a node in lattice B located at ′ ′ i j ( , )only if | − ′| ⩽ i i r and | − ′| ⩽ j j r. Of the + r (2 1) 2 nodes that fulfill this condition for any given node, one is chosen uniformly at random. If = ≡ ∞ r L , we recover the case of purely random dependency links as described in [18,20]. The unconstrained case is an important limiting case for which we can develop more complete models and theoretical predictions due to the uniform distribution of damage from one network to the other.
In contrast to previous work on interdependent networks, we consider a node to be functional only if there is current flowing through it and also through the node that it depends on, see figure 1(a). However, as mentioned above, no current flows between the layers. As such, only part of the current-carrying backbone in each network can be functional at any moment. The entire backbone will generally not be functional because some nodes will depend on non-current carrying nodes in the other network (cf the green nodes in the upper right and lower left of figure 1(a)). In each network, we determine the backbone via direct electrification [3,22]. We then disable any nodes which are not in the backbone (the blue or green squares in figure 1) or are dependent on nodes which are not in the backbone of the other network. The removal of nodes due to dependencies causes a dramatic decrease in the backbone size from iteration to iteration (see figure 6) and the process is continued until it reaches a steady state, (cf the red nodes in figure 1(a)). At this point the bulk conductivity σ and backbone fraction ∞ B are measured.

Results
We find that, for every system, there exists a critical fraction p c of nodes above which a spanning cluster exists and current flows through both networks and below which it does not. The phase transition is second-order for < r r c , and p c increases monotonically with r. For ⩾ ≈ r r 4 c , the transition is first order and p c decreases asymptotically towards its value for = ∞ r , see figures 2(a) and 3(a). In contrast, fully interdependent lattices with structural dependency undergo a first-order phase transition only when > ≈ r r 8 c [19,20]. In all cases, the p c of process-based dependency is substantially higher (more vulnerable) than structural dependency, due to the fact that the current-carrying backbone is always smaller than the GCC, see figure 3(a). This also leads to a larger discontinuity in the order parameter compared to structural dependency, as shown in figure 3(b).
As in the case of structural dependency of interdependent spatially embedded networks [19][20][21][22], we find three regimes of critical behavior which are determined by the value of r. In the first regime, for < r r c , the phase transition is driven by the spontaneous emergence of many small holes, each triggering localized cascading failures in their vicinity. Because the size of the small holes increases continuously as → + p p c [5], ∞ B and σ decrease gradually and there is a continuous second-order transition. In the second regime, for ⩽ < ∞ r r c , the dynamics at criticality are very different. Random fluctuations generated by the original node removal lead to the emergence of a hole of a critical size. The same hole will appear in approximately the same place in both networks due to the restriction of dependency links to length r. This hole causes extensive damage to nodes up to a distance r from its edge. At criticality, this damage leads to the separation of small clusters around the edge of the hole and its gradual propagation. This continues until the networks no longer have a spanning cluster and the flow of current is entirely cut off. In this regime there is no critical exponent above p c , the number of nodes which fail per iteration grows linearly until the hole reaches the system's boundary, and the number of iterations at criticality grows linearly with system size.
In the third regime, for = ∞ r , there is no critical hole which propagates. Instead, the system undergoes a mixed first and second order transition, similar to what is found in interdependent random networks [14,26,27]. This transition is characterized by a brief initial cascade in which a large fraction of the nodes fail followed by a long plateau during which the number of failing nodes stays approximately constant from step to step. After the slowly accumulating damage during this phase reaches a sufficient amount, the damage rate increases exponentially and the entire system collapses. Thus the majority of the nodes collapse at the beginning or end of the cascade while the majority of the time is spent in the plateau phase. This mechanism is analyzed in detail for the case of interdependent random networks in [27] and results are compared with spatially embedded networks with structural dependency in [19].
When the dependency links are of unconstrained length ( = ∞ r ), we can accurately predict the backbone fraction ∞ B and conductivity σ for any value of p. This is because damage from one network causes damage in the other network that is distributed evenly. Since the requirement that the nodes be members of the current carrying backbone ∞ B is analogous to the requirement that they both be in the GCC ( ∞ P ), we can predict the evolution of ∞ B in the same manner that ∞ P is calculated in [20] based on [13]. Since the size of the backbone is the same in both networks, it is sufficient to calculate it for either network. The fundamental observation underpinning this theory is that-when damage is randomly distributed-each iteration can be treated as a percolation problem with the initial p determined by the network size of the other network in the previous iteration. We let p i represent the fraction of nodes expected to survive at iteration i before including the effect of percolation and ∞ B p ( ) i be the expected size of the backbone at iteration i. Throughout we use ∞ B to describe the actual backbone fraction and ∞ B x ( )to describe the functional relationship between the fraction of nodes removed ( − x 1 ) and the backbone fraction ∞ B in single networks, as plotted in figures 2(b) and 4(a). The size of the current carrying backbone is thus: i i i 1 1 We obtain the steady state solution to equation (2) by imposing = − p p i i 1 and defining the p i that satisfies this as x:  When r = 0, the dependency links have no effect and the two systems have the same p c as in single 2d lattices ( ≈ p 0.5927 c [24]). However, for all other values of r, the kind of dependency has a strong impact on the overall robustness. The peak in p c occurs at = ≈ r r 4 c for the process-based dependency compared to around = ≈ r r 8 c for the structural dependency [20]. This peak divides the second-order transitions (on the left) from the first-order transitions (on the right). (b) The order parameter ∞ P , ∞ B at p c as a function of r. The order parameter is ∞ P for structural-dependency and ∞ B for process-based dependency. Though the transition is second-order for < r r c , there is still a finite jump in the order parameter because of the finite system size. The reason that the two curves do not coincide for r = 0 is that, for every finite system at p c , < ⩽ ∞ ∞ B P 0 . For both panels L = 1000, Δ ⩽ × − p 5 10 6 and − 50 100 realizations were averaged for each data point. The error bars denote standard deviation.
Thus for any p, if x is a solution to equation (3), ∞ B x ( )is the fraction of nodes in the backbone, in each network. If the initial removal of nodes is performed on both networks, the same equation holds with p replaced by p 2 in equation (3).
To determine σ, we utilize the fact that, on a single network or on a pair of networks with = ∞ r , the conductivity has a one-to-one relationship with the size of the backbone σ ∞ B ( ), see inset of figure 4(a). It is then straightforward to predict the conductivity σ, see figure 4(b). We can also predict p c by finding the first nontrivial solution to equation (3) and the expected values of ∞ B or σ at each iteration, see appendix. When < < ∞ r 0 , we are unable to make direct predictions based on the percolation profile of a single network. This is due to the fact that damage is not distributed uniformly and the phase transition itself is characterized by the emergence of a critically sized hole which grows from iteration to iteration and eventually destroys the whole system [19,20]. Another consequence of the non-uniform distribution of damage in systems with finite r is that the relationship between ∞ B and σ is no longer one-to-one. Therefore even if we had a way to predict the size of ∞ B , we would not be able to precisely deduce the value of σ from it. Scaling properties of the backbone and conductivity above p c on a single lattice have been studied extensively [10]. For a single lattice, the conductivity above p c scales as in which ≈ t 1.31 [28,29] and the size of the backbone scales as in which ν = 4 3 is the correlation length critical exponent for two-dimensions [5,30,31], ψ ≈ 0.3568 and ψν ≈ 0.4758 [10,29,32].
In figure 5, we show the scaling above p c for the backbone and sigma respectively for interdependent resistor networks with low values of r, where the transition is still second-order. We find that, despite the dramatic change in p c as a function of r ( figure 3(a)) and the difference in σ and ∞ B outside the critical region ( figure 2(a)), there is no appreciable difference in the scaling exponents. However, we do see that as r increases, the extent of the scaling region decreases, and already for r = 3 the behavior has broken down (see figure 5). We believe that this is due to the fact that as r approaches r c , the dynamics of the transition becomes a hybrid of the single-lattice second-order transition and the hole spreading first-order transition. Instead of a single hole propagating regularly, multiple holes propagate irregularly. This mixture of critical behaviors leads to a blurring of the critical region. Nevertheless, we expect that for larger systems where p can reach closer to p c , the scaling region will increase. In the absence of dependencies, the value of p c is the same whether functionality is determined by membership in the GCC or the current carrying backbone, as represented by fractions ∞ P and ∞ B , respectively. The conductivity (σ) is non-zero only when a current carrying backbone exists and therefore p c is the same. The inset shows the relationship between ∞ B and σ which holds for both single networks and interdependent networks with unconstrained dependency links. (b) In the presence of dependencies, the difference between determining functionality by structure alone or based on the current flow process has a strong impact on the value of p c . Heuristically, we can consider the single network case as degenerate in p c with regard to backbone or spanning component and that this degeneracy is removed when the dependency interaction is introduced. ∞ P and ∞ B represent the fraction of nodes which are in the spanning component or backbone of their own network, respectively. Due to symmetry, ∞ P and ∞ B are the same in each network. The lines are calculated based on the theory described in the text, in equation (3). The symbols on the predicted lines are the result of over 130 simulations on systems of size L = 1000.

Discussion
In this study, we have analyzed the first model of interdependent networks in which nodes can only function and provide support when they are in a dynamically active state. To realize this scenario, we have revisited the canonical study of percolation on resistor networks as a classic example of general transport phenomena. We find that process-based dependency leads to substantially increased vulnerability (figure 3(a)) including higher p c values and a first-order transition emerging for lower values of r, but that the scaling exponents at p c of ∞ B and of σ are indistinguishable from the exponents of single networks ( figure 5). This isconsistent with the hypothesis that coupled lattices with < r r c and single lattices are in the same universality class and is particularly striking in light of the pronounced difference in the σ p ( )curves outside the critical region ( figure 2(a)). We find that in the absence of any interaction between the networks, there is a kind of degeneracy: the percolation threshold p c is the same whether functionality is determined by connectivity or current flow. However, once an interaction is introduced, this 'degeneracy' is removed and the two types of functionality have different percolation thresholds. This is shown theoretically for random links (figure 4) and numerically for links of limited length ( figure 3(a)).
Resistor networks are significant because they present a fundamental model for general transport phenomena. We believe that process-based dependency may be a more realistic and useful approach for dealing with many systems in which it is important to differentiate between nodes which are connected and nodes which are dynamically active. For instance, people in a social network or organisms in an ecological network may be linked according to some measure but those links may not reflect their actual functionality [33,34]. Furthermore, because the example of process-based dependency on interdependent resistor networks is based on a real physical process, we propose a physically realizable system of interdependent networks based on metal insulator transitions.
Many materials, especially transition-metal oxides, undergo a metal insulator transition at a critical temperature T c from a low-temperature semiconductor phase to a high-temperature metallic phase [35,36]. In a planar system composed of grains of a transition-metal oxide connected by metallic links, dependency can emerge in the form of heat transmission. If the system is cooled from an initial temperature above T c , even though the ambient temperature is too low to sustain conductivity, nodes can still conduct electricity if they receive heat from nearby conducting nodes-even from nodes with which they are not electrically connected. In such a system a warm backbone will emerge in each network which will shrink (as in figure 1) as nodes are disabled or as the background temperature decreases [37]. In a single network, bulk conductivity will continue for temperatures substantially below T c . With process-based dependency, however, since the red bonds in each network will generally not be co-located, this will cut off the conductivity sharply at a new ′ > T T c c . Similarly, if the nodes change conductive properties in the presence of certain gas molecules, the first-order transition can be triggered by the presence (or absence) of gas molecules. Whether driven by temperature or gas molecules, the sharp first-order transition described here can be used to construct a new kind of highly sensitive sensor.  figure 2(a), the scaling behavior shows little variation as r is varied.This is consistent with the hypothesis that interdependent lattices below r c are in the same universality class as single lattices. The case of r = 3 shows substantial deviation from the other cases. This is because it is not a pure second-order transition like = r 0, 1, 2 but rather combines first and second-order properties as discussed in the text and is visible in figure 3(b). Still, we see that in a small critical region, it has the same critical behavior. For = r 0, 1, 2, 3, we find t values of 1.309 ± 0.001, 1.309 ± 0.003, 1.311 ± 0.008 and 1.31 ± 0.046, respectively. Likewise, we find β values of 0.476 ± 0.001, 0.476 ± 0.002, 0.476 ± 0.004 and 0.480 ± 0.015, respectively. These results were obtained from over 100 realizations of systems with L = 1500 and fitted using orthogonal distance regression on all data points in the critical region to accurately measure errors [25]. Log binned points are shown for visibility only and are not used in the fitting process.
A number of technical challenges remain before such a physical experiment can be conducted.The nodes need to be prepared in such a way that the absence of heat from the resistive dissipation in one node triggers a complete loss of conductivity in another node, which may prove challenging. Alternatively, the model can be extended to describe a relationship in which a lack of current in a node causes a continuous decrease in the conductivity of the node it depended on. A greater challenge is to physically implement the dependency link. This is not trivial in the current model because we require ⩾ r r c for the first-order transition to take place. If the physical dependency is based on heat transmission, making dependency links of length ⩾ r r c would require heat screening between layers which would distribute the heat in a disordered but localized fashion. We believe that the best way to overcome these technical problems is to extend the model of resistor networks to other topologies for which dependency link creation is more straightforward. We believe that such extensions are feasible and would be an important topic for future research both for the purpose of the physical experiment described here as well as for providing deeper understanding of the differences between process-based and structural dependency.
For comparison, we also present the results for r = 5 where we do not have a theoretical prediction for ∞ B or σ at each iteration. Since > r 5 c we are in the spreading regime. The reason for the parabolic decrease in ∞ B (figure 6(c)) is that the critical hole spreads at a constant rate. Thus at each iteration the radius of the hole increases as t and its area (which is of order − ∞ B 1 ) as t 2 . We see similar results for the conductivity of the system ( figure 6(a)).