Multiple Tipping Points and Optimal Repairing in Interacting Networks

Systems that comprise many interacting dynamical networks, such as the human body with its biological networks or the global economic network consisting of regional clusters, often exhibit complicated collective dynamics. To understand the collective behavior of such systems, we investigate a model of interacting networks exhibiting the fundamental processes of failure, damage spread, and recovery. We find a very rich phase diagram that becomes exponentially more complex as the number of networks is increased. In the simplest example of $n=2$ interacting networks we find two critical points, 4 triple points, 10 allowed transitions, and two"forbidden"transitions, as well as complex hysteresis loops. Remarkably, we find that triple points play the dominant role in constructing the optimal repairing strategy in damaged interacting systems. To support our model, we analyze an example of real interacting financial networks and find evidence of rapid dynamical transitions between well-defined states, in agreement with the predictions of our model.

M ost real networks are not isolated structures but interact with other network structures. As a result, much research has been focused recently on the dynamics of interdependent 1-8 and multilayer [9][10][11] networks. Recent studies on network repair [12][13][14] have shown the importance of recovery of nodes as a process that leads to reverse transitions, hysteresis effects and such phenomena as spontaneous recovery 12,15 .
The cardiovascular and nervous systems in the human body are examples of two dynamically interacting physiological networks 16 . Diseases often result from complex pathological conditions that involve a dynamical interaction with positive or negative feedback between different functional subsystems in the body. Similarly, in the global economy there is a hierarchy of clustered and tightly connected countries, often grouped geographically, which are further interconnected to one large global interacting economic and financial network [17][18][19] . To understand the behaviour of these systems using network science, we develop a model of interacting networks with nodes that can recover from failure and we examine the resulting phase diagram.
Our model of a generic system consisting of interacting dynamical networks captures the important events found in real-world interacting networks, that is, node failure [20][21][22][23] , systemic damage propagation 24 and node recovery 12,15,25 . We first analiticaly solve the model in the mean field theory (MFT) approximation and confirm the results in numerical simulations. The phase diagram of this system is characterized by a number of phases, phase transition lines and tipping points [26][27][28][29] . We show that triple points play a critical role in devising an optimal strategy for efficient repairing of interconnected network systems. By analysing the network of credit default swaps (CDSs) of sovereign countries, we demonstrate the application of our model to a real system and show that all of its parameters are experimentally accessible.

Results
Model. The structure of our system for the n ¼ 2 case is modelled as follows. We start with two isolated networks, network A and network B, and for simplicity we assume that both networks have the same number of nodes N and the same degree distribution f(k) (these assumptions can be relaxed but the results stay qualitatively similar). We assume that within each network the nodes are randomly connected. Now, to allow networks A and B to interact, we introduce interdependency links that connect nodes across the two networks 2 . This can be achieved in different ways and we use a simple one-to-one dependency: each node in network A is dependent on exactly one node in network B and vice versa. The pairs of nodes of both networks are chosen randomly.
The dynamic behaviour of our system is governed by two categories of event-failure and recovery-and we assume that every node is in either a failed or an active state. Node failure can result from internal failure or from the spread of damage from neighbour nodes in either the same network or the interdependent network. We thus assume that there are three ways a node can fail. The first way is the internally induced failure, when a node's internal integrity has been compromised, for example, an organ in the body can fail due to a malfunction within the organ or a company can fail due to bad management. The second type of failure is externally induced failure through failure propagation due to connections with failed nodes within the node's own network. Finally, there is a failure induced through the dependency link as a result of being dependent on a failed node from another (opposite) network. Apart of these three types of failures, we assume the existence of associated simple recovery processes for every type of failure. We specify quantitatively each of these processes below.
For internal failures (I), we assume that in both networks any node can fail due to internal problems, independent of other nodes. For each node in network A we assume that there is probability p A dt that the node will fail internally during any time period dt. The equivalent parameter in network B is p B .
Every node in network A and network B is connected by links to nearby nodes in its own network. These nodes constitute the node's neighbourhood. The number of links a node has within the network indicates its degree or connectivity, denoted by k. If a large number of nodes in a node's neighbourhood have failed, that is, if the neighbourhood is substantially damaged, we assume that the probability that the node itself will fail is increased. This is modelled by external failures (E). As in refs 12,30, we use a threshold rule to define a substantially damaged neighbourhood, which is a neighbourhood containing rm active nodes, where m is a fixed integer threshold. If node j has 4m active neighbours during time dt, we consider its neighbourhood to be 'healthy' and there is no risk of external failure. On the other hand, if j has rm active neighbours during time dt, there is a probability r A dt (for network A) or r B dt (for network B) that node j will externally fail. For certain systems it is more appropriate to define a fractional threshold 0rm frac r1 as in ref. 15. That is, the minimum number of active nodes as a requirement for a 'healthy' neighbourhood is replaced by a minimum fraction of active nodes in the neighbourhood. In the example of random regular network that we consider below, both are equivalent and related by m ¼ km frac .
In the case of two interdependent networks (A and B), we assume that each node in the first network is dependent on a node in the second network via an interdependent link and vice versa. We assume that if one node in the pair fails, there is a finite (but not 100%) probability, r d dt, that during time dt the other node in the pair will fail as well (dependency failure: type D). This represents the probability that the damage will spread through the interdependency link.
We also assume that there is a reversal process, a recovery from each of these three types of failure. A node recovers from an internal failure after a time period t 6 ¼ 0, it recovers from an external failure after time t 0 and from a dependency failure after time t 00 . In simulations, and without loss of generality, we use t¼100 and for simplicity we set t 0 ¼t 00 ¼1, to take into account the assumption that real-world systems usually require a longer time period to recover from internal problems (physical faults) then from a lack of environmental support. However, changing the numerical values does not introduce any qualitative difference. For the node activity notation, we assume that every node is in one of two states: active or failed. A node is considered active in the observed moment, if it is not experiencing internal (I), external (E) or dependency (D) failure.
Mean field theory. We characterize this system by studying the order parameters chosen naturally as the fraction of active nodes in network A and network B, z A and z B , respectively. To simplify the calculation, however, we first concentrate on the complementary and equally intuitive fraction of failed nodes a A and a B , in networks A and B, respectively Using the MFT presented in Methods, we obtain two coupled equations that connect a A and a B , which the system must satisfy in the equilibrium Here F k; x ð Þ¼ P m j¼0 ð k j Þx k À j 1 À x ð Þ j and we have also introduced simplifying parameters p Ã A 1 À e À p A t and p Ã B 1 À e À pBt , to make the equations more elegant and to reduce the number of parameters by replacing p A , p B and t that appear as a product. We find that the parameters p Ã A and p Ã B are very convenient to work with, because they correspond to the fraction of internally failed nodes in network A and network B, respectively.
Despite the seeming complexity of equations (1) and (2), it is noteworthy that there are only two unknown variables, a A and a B , and that all other parameters are fixed. These two equations define two curves in the (a A , a B ) plane. Figure 1a shows a graphical representation of the curves for a random regular 31 network (in which all the nodes have the same degree) with degree of k ¼ 16 and threshold m ¼ 8, for the symmetric parameter values p Ã A ¼p Ã B ¼0:16, r A ¼ r B ¼ 0.60 and r d ¼ 0.15. The size of each network is N ¼ 2 Â 10 4 . The blue curve is a graphical representation of equation (1) and the brown curve is defined by equation (2). The curves, similar to two 'ropes', create a 'knot' that can have up to nine intersections, representing mathematical solutions of the system of equations. However, not all of these solutions represent observable and stable physical states. To see that, observe one of the curves in Fig. 1a, for example, the blue curve described by equation (1). If we increase damage done to network B (that is, we increase a B ) and keep everything else constant, some damage will undoubtedly spread to network A. Thus, we expect that when a B is increased, a A must also increase (it would be very unusual if one network improves its activity as a result of damaging the other network). We conclude that the parts of the blue and brown curve that produce physical solutions are only those where a A and a B increase together or decrease together along the curve. This elimination leaves only four states in Fig. 1a that are stable (green circles), whereas the other five states are unstable (red crosses), for this particular choice of parameters. In simulated finite networks, when the network system evolves according to the rules of the model, at t ¼ 0 we have a freedom to set initial conditions for the activities. Systems initially prepared to have a pair of values (a A , a B ) corresponding to an unstable solution of equations (1) and (2) will be disturbed by a small fluctuation of a A or a B , owing to the system dynamics, and the values of a A or a B will rapidly change until one of the stable states is reached. Systems that are initially prepared to have values of a A or a B corresponding to a stable solution will fluctuate around these values, until perhaps a large finite fluctuation occurs and the system 'jumps' to another stable state. In general, for any choice of parameters, we have between one and four stable (physical) states. Figure 1b shows the scenario for the same network system when p Ã A ¼0:20, p Ã B ¼0:24, r A ¼ r B ¼ 0.60 and r d ¼ 0.15. In this case we have two stable states and one unstable state.
In the Methods, section 'Additional phase diagrams', we show diagrams for z A ¼ 1 À a A for a range of different values of p Ã A and all other parameters fixed. This MFT calculation agrees well with the states that we observe in our simulations, as we will demonstrate below.
The four stable solutions found above correspond to the following four scenarios: '11' or 'up-up', when there is high activity in both network A and network B; '12' or 'up-down' when there is high activity in network A and low activity in network B; '21' or 'down-up' when there is low activity in network A and high activity in network B; and '22' or 'down-down', when there is low activity in both network A and network B.
Depending on the parameters, we obtain between one and four stable states. Each of the states exists in a certain volume of the multi-dimensional space of parameters. Results of the MFT calculation for a particular set of parameters are presented in Fig. 2a-d as a phase diagram with four layers. Figure 2 shows the regions in which each of the four states exist in the (p Ã A , p Ã B ) parametric sub-space, when other parameters are fixed at values r A ¼ r B ¼ 0.60 and r d ¼ 0.15, with k and m remaining the same as before.
For example, in Fig The state in which the system is found depends on the initial conditions or the system's past. There are a total of 10 different transitions (11-12, 11-22, 11-21, 12-11, 12-22, 21- ARTICLE discovered in multiplex networks. In particular, Baxter et al. 32 introduced weak bootstrap percolation and weak pruning percolation in multiplex networks, which have potential applications in infrastructure recovery and information security, and can even provide a way to diagnose missing layers in a multiplex network. We next can examine the activity profile for various crosssections in the phase diagram. In Fig. 3 Figure 4a shows the activity measured in simulations of network A as we move along the black dashed line, changing both p Ã A and p Ã B , and preserving the relation p Ã B ¼0:1 þ 4=3p Ã A . We perform simulations for various initial conditions and find (Fig. 4a) three different states denoted by green, orange and blue colours, which we identify as 11, 12 and 22 states, respectively. We find four different transitions: 11-12, 12-22, 12-11 and 22-12. The solid lines show the MFT prediction (equations (1) and (2)) for the activity of network A.
The good agreement shows that the MFT correctly captures all the properties of the system. We note that qualitative agreement between the MFT and the simulations is better for higher values of k, because for higher k the fluctuations are smaller, which improves the accuracy of the MFT. Figure 4b shows the activity when moving along the red dashed line. Here we obtain four states and six different transitions.
The phase diagram of a system of n ¼ 2 interacting networks ( Fig. 3) is much richer than the phase diagram of a single network with damage and recovery 12 . The analytical results we presented here for n ¼ 2 can be generalized to n interacting networks in any topological configuration, although as n increases they become increasingly difficult to visualize. In general, a system with n interacting networks can have up to 2 n physical states.
The problem of optimal repairing. Knowing and understanding the phase diagram of interacting networks enable us to answer some fundamental and practical questions. A partially or completely collapsed system of nZ2 interacting networks in which some of them are in the low activity state is a scenario common in medicine, for example, when diseases or traumas affect the human body and a few organs are simultaneously damaged and need to be treated, and the interaction between the organs is critical. It is also common in economics, when two or more coupled sectors of the economy 18 experience simultaneous problems, or when a few geographical clusters of countries experience economic difficulties. The practical question that arises is: what is the most efficient strategy to repair such a system? Many approaches are possible if resources are unlimited, but this is usually not the case and we would like to minimize the resources that we spend in the repairing process. For simplicity, consider two interacting networks, both damaged (low activity). Is repairing both networks simultaneously the more efficient approach, or repairing them one after the other? What is the minimum amount of repair needed to make the system fully functional again? In other words, what is the minimum number of nodes we need to repair, to bring the system to the functional 11 ('up-up') state, and how do we allocate repairs between the two networks? An optimal repairing strategy is essential when resources needed for repairing are limited or very expensive, when the time to repair the system is limited, or when the damage is still progressing through the system, threatening further collapse, and a quick and efficient intervention is needed.
We show below that this problem is equivalent to finding the minimum Manhattan distance between the point in the phase diagram where the damaged system is currently situated and the recovery transition lines to the 11 region. The Manhattan distance between two points is defined as the sum of absolute horizontal and vertical components of the vector connecting the points, with defined vertical and horizontal directions. It is a driving distance between two points in a rectangular grid of streets and avenues. In our phase diagram, it is equal to Dp Ã . It turns out that two triple points of the phase diagram play a very important role in this fundamental problem. We find that these special points have a direct practical meaning and are not just a topological or thermodynamic curiosity.
To show this, we start by making some simplifying but reasonable assumptions. First, we assume that only internal failures can be repaired by human hands, as these failures are physical faults in nodes (any external and dependency failures and recoveries are 'environmental', and are a spontaneous recognition of the changing neighbourhood of a node). We mentioned above that the parameters p Ã A and p Ã B correspond to fractions of internally failed nodes in networks A and B, respectively. This implies that the number of internally failed nodes repaired in, say, network A, is directly proportional to the change of p Ã A . Hence, repairing nodes in networks A and B means decreasing p Ã A or p Ã B . We also assume that these repairs are done fast enough that there is only a small probability that the newly repaired nodes will internally fail again before the repair process is completed. The total number of repaired nodes is therefore and it is proportional to the Manhattan distance between the starting and final point in the phase diagram.
To optimize repairing we need to minimize this metric. Figure 5 shows the solution to the minimization problem and a detailed discussion is provided in the Methods section. The different colours in Fig. 5 correspond to the different optimal repair strategies, which depend on the failure state of the system. If the system is initially at point S 1 , both networks are in a low activity state, that is, they are non-functional. Our goal is to decrease p Ã A and p Ã B , and arrive to the region where the system is fully recovered (the green region) by performing a minimal number of repairs, that is, minimal N rep . We find that for any point in the red region there are actually two closest points in the green region, at an equal Manhattan distance away from the red region point. These two points are the triple points R1 and R2 shown in Fig. 5, which also correspond to the triple points in Fig. 2b. Although R1 may be closer to point A than R2 by Euclidian distance, the Manhattan distance is the same. Thus, two equally good repairing strategies are available. One involves allocating more node repairs to network A and the other allocating more repairs to network B. For the yellow regions (points S 2 and S 3 ), the closest points by Manhattan distance are R1 (for point S 2 ) or R2 (for point S 3 ). Here, only one triple point represents the optimal solution. It is noteworthy that the path samples in Fig. 5 are 'zig-zag' in shape (to highlight that we are minimizing Dp Ã A þ Dp Ã B ); however, even when a diagonal path (direct straight line) to a triple point is used, the Manhattan distance is the same. For the dark blue regions (points S 4 and S 7 ), the optimal strategy is to decrease p Ã B only, until the system is recovered. Similarly, for the light blue regions (points S 5 and S 6 ), the optimal strategy is to decrease only p Ã A . From our optimal repairing strategy analysis we find that the order of repair (the specific path taken between the initial point and final point) does not affect the final result. Minimizing the Manhattan distance only determines the optimal destination point. Therefore, there is actually a set of paths corresponding to equally optimal repairing processes.
States and transitions in real-world networks. In relatively small networks (NE10-1,000) fluctuations are very large. Thus, in small network systems exhibiting multistability it is possible to observe phase flipping 12,15,33 between different states. Figure 6a shows the fraction of active nodes for both networks, in time, for where the collapsed system is initially situated (S i ) to the nearest border of the green region where it becomes fully functional. For a system having the initial condition within the red section (for example, point S 1 ), there are two solutions: it is equally optimal to reach any of the two triple points R1 and R2 by decreasing p Ã A and p Ã B . For the systems starting in the yellow regions, it is optimal to reach only one triple point, R1, for the sector containing point S 2 , or R2 for the sector containing point S 3 . Starting in the dark blue regions it is optimal to decrease p Ã B only, that is, repairing only network B. Similarly, in the light blue regions it is optimal to decrease p Ã A only. Triple points play a crucial role when both networks are initially significantly damaged (red and yellow regions). switching between well-defined high-activity and low-activity states, as well as correlated collective behaviour of the two networks in interaction. We identify collective states 11, 22, 12 and 21, and mark them with connected black ovals. It is noteworthy that, as the CDS value grows with risk, a higher activity in a CDS network corresponds to bad economic news. a symmetric choice of parameters, p Ã A ¼p Ã B ¼0:21, r A ¼ r B ¼ 0.60 and r d ¼ 0.15, when each network has only N ¼ 100 nodes. Large fluctuations cause the system to jump between the different states allowed for this set of parameters. It is noteworthy that interdependent links cause the two networks to have partially dependent and correlated dynamics. Very often a transition in one network triggers a transition in the other. In Fig. 6a we can identify examples of all four global states: 22, 11, 21 and 12. For example, at time tE400 both networks are in the high activity state (11), whereas at tE620 network A is in the low activity and network B in the high activity state (21).
As many real-world interacting network systems have a small number of nodes, in those systems we can potentially uncover dynamics similar to what we observe in our model networks. As an example of a real system, we investigate the interacting sovereign 5-year CDS system, consisting of 25 European Union (EU) and Latin American countries (see Methods for the full list of countries) that began to issue the CDSs from 2005. We divide countries into two groups on a geographical basis: 8 countries belong to the Latin American group and 17 belong to the European Union group. Sovereign CDSs are financial instruments, for which the value reflects the probability that the reference country will default on its debts. Each country has a CDS value assigned and this value changes in time reflecting the economic news about this country and the perceived risk of default, which results in a time series that we can observe. CDSs are highly sensitive to important economic news, positive or negative. There is also a significant contagion and influence between the countries, especially between those with strong economic ties, which is reflected in the correlation between their CDSs. These characteristics make the CDS signals a candidate for modelling using our interacting network approach.
We can draw a parallel between the CDS system and our model network if we assume that each country (with its associated CDS signal) can be represented as a node, which has connections (links) to other countries within its own geographical region, as well as ties with countries from another continent. In this case, we might expect that random and independent bad (or good) economic news appearing in any given country have behaviour similar to random internal processes in nodes in our artificial model (random internal failures/recoveries). When economic problems in one country propagate to a neighbouring country within the same geographical region, the process resembles the external failures in our artificial model, whereas interaction between countries from different continents may be modelled by the interdependent links from our network model. For the CDS network system we also suppose that the fractional definition for the threshold (m frac ) is somewhat more natural then the absolute definition, as it is less dependent on the country size or its importance, that is, the number of links a country has to other countries.
We study the international CDS system during the period between June 2005, the earliest date when CDSs traded for all countries, and February 2014. We apply the network model to it as follows. We represent each country with one node that can have two states: active or failed. As the raw CDS values are continuous by nature and our model uses binary node states (up or down), we perform a trend mapping procedure to form a binary signal (0 or 1) for each country. In particular, for each time t, we consider the interval [t À 252, t] of 252 business days (the usual number of business days in a year). If the CDS value of a country has a net increase during that period, we consider the node of the country to be active at t (state ¼ 1). If it does not, it is inactive (state ¼ 0). Having individual binary signals for each country, we can calculate the average activity 0rz(t)r1 for both EU and Latin American networks. The resulting time series for EU and LA activities are shown in Fig. 6b. First, we note that the two geographical networks spend most of time having either a significantly high activity or significantly low activity (that is, there is an indication for two well-defined single-network states). We confirm this by measuring the frequency distribution of network activities (Fig. 7a,b), which exhibit a strong bimodality in z. The CDS network system in Fig. 6b shows rapid transitions between the high and low activity states, much similar to the artificial network system in Fig. 6a. Figure 7c shows the calculated correlations between binary signals of pairs of individual nodes. The correlation matrix reveals two strongly correlated blocks, which we identify as Latin American block (numbers 1-8) and EU block (numbers 9-25).
In Fig. 6b, we also observe that the two networks sometimes make transitions simultaneously, but not always. This behaviour also resembles the behaviour observed in the artificial networks in Fig. 6a.
Finally, we find that it is possible to estimate numerical values for all the model parameters of this real system (internal p Ã EU , p Ã LA ; external m frac,EU , m frac,LA , r EU , r LA ; and interdependent r d ) from the data. The basic idea is that for each parameter we identify an observation experiment in which this particular parameter dominates, enabling us to effectively isolate individual parameters from the noise of many others. For example, when both networks (EU and LA) are in the high activity phase, most of the failures are in fact internal failures. This allows us to almost directly estimate p Ã EU and p Ã LA from Fig. 6b, by observing p Ã EU ¼1 À z EU h i and p Ã LA ¼1 À z LA h i. External failures are most significant when a network is in a low activity state. Interdependent parameter r d can be estimated by studying the correlation between z EU (t) and z LA (t), as this is an increasing function of r d . Threshold parameters m frac,EU and m frac,LA can be estimated by exploiting the fact that they most significantly determine the fraction of time that each network spends in the high, or low, activity states. We describe in detail our procedures for numerically estimating these model parameters in the Supplementary Methods. The procedure for estimating some of these parameters is also illustrated in Supplementary Fig. 1. Numerical results for the parameter estimates are presented in Supplementary Table 1. Our dynamical network model also independently predicts that the typical fluctuation size of z(t) is not uniform for all values of z, but has a spike around z % 1 2 . We observe this phenomenon in both our simulations and the real network dynamics (Supplementary Fig. 2).

Discussion
Interacting networks appear across many disciplines, from medicine, physics, biology and ecology, finance and economics, to infrastructure (refs 34,35). We propose a generic model that captures some of the most common processes found in real interacting networks-node failure, systemic damage propagation and node recovery. We report several intriguing results. Our solution of the model produces a rich phase diagram with a number of tipping points (critical points, triple points and transition lines). Using the phase diagram, we solve a fundamental problem of the optimal repairing strategy for a damaged system consisting of interconnected networks. Solving this problem enables us to determine the minimum set of node repairs required to repair a failed interconnected network system; thus, it becomes fully functional again. Remarkably, we find that the triple points from our phase diagram play the dominant role in constructing the optimal repairing strategy in damaged interacting systems. This implies that triple points are not only a thermodynamic or topological curiosity, but they have a very direct practical meaning and application, specifically in constructing repairing strategies in interacting systems. The problem of functional repair using minimal intervention is relevant in medicine, economics and other disciplines, when repairs are invasive or expensive. Finally, we apply our model to a selected real system: interconnected networks of CDSs, for two interconnected groups of countries. We propose a methodology to measure all of the model parameters in a real system, by using observational experiments in which a particular parameter dominates the behaviour of the system (see Supplementary  Information).

Methods
Damage conductivity parameters. Parameters r A and r B are introduced, because they describe how easily the damage is spread through the network. When r ¼ 0 there is no damage spread between the nodes, and when r ¼ 1 there is perfect damage conduction. Assuming that external failures occur with certainty would mean fixing r to be equal to 1. In the case of a single network with recovery, it has been shown 12 that many important phenomena (for example, spontaneous recovery) are lost when r ¼ 1. The most interesting parts of the phase diagram are in fact where r is far from 1.
It is noteworthy that our choice of the value for r d is quite limited. If r d is too large, we find that the damage spreads through dependency links extremely efficiently and the only possible stable state is total system collapse. The extreme  (1) and (2), for a range of p Ã A values, in a system of two interdependent networks (k ¼ 16, m ¼ 8). Blue lines correspond to stable physical states and red dotted lines represent unstable solutions. In Fig. 8a, the same parameters as in Fig. 1a are used, except p Ã A , which is not fixed but varied. When p Ã A ¼0:16 (vertical black line), the corresponding values on the blue solid lines (green circles) match the graphical solutions in Fig. 1a (also green circles). (b) An analogous relationship holds between Figs 1b and 8b, in which case two stable states exist.
vulnerability of interdependent networks is well known 2,22 . As there is always at least one functional stable state in biological or man-made systems, total system collapse as the only stable state is not realistic. Thus, we need the r d parameter to 'soften' the dependency links 22 and allow a more realistic behaviour and its value should not be close to 1.
Mean field theory. Fractions a A and a B denote the fraction of nodes that are failed due to any of the three types of failures: internal (I), external (E) or dependency failure (D). We denote the probabilities that a node at a time of observation experiences a failure of I, E or D type as P(I), P(E) and P(D), respectively. As a first approximation, we assume that these failures are mutually independent events. Considering network A first, we write an expression for the probability a A,k that a node of degree k in network A is failed. The node can fail due to I, E or D events, or to a combination of them. Using the inclusion-exclusion principle for independent events, we write Next, we separately calculate P(I), P(E) and P(D). P(I) is also the average fraction of internally failed nodes in a network, as internal failures are independent events. This is a Poisson process on individual nodes 12,36 and therefore P I ð Þ¼e À pAt . As parameters p A and t come in this expression as a product, we can replace them with a single parameter, p Ã A e À pA t , which is bounded and also has the property 0 p Ã A 1 and thus P I ð Þ¼p Ã A for network A.
Next, we calculate P(E), the probability that a randomly chosen node with degree k has externally failed. Focusing once again on network A, without a loss of generality, we let F(k) be the probability that a node of degree k in network A is located in a critically damaged neighbourhood (where fewer than m þ 1 nodes are active). By definition, the time-averaged fraction of failed nodes (for any reason) in network A is 0ra A r1. In a mean-field approximation, this is also the average probability that a randomly chosen node in that network has failed. Using combinatorics, we obtain F k; a A ð Þ¼ P m j¼0 ð k j Þa 12). The probability that a node of degree k in network A has externally failed is then P(E) ¼ r A F(k, a A ). An analogous result is valid for network B.
Finally, we calculate P(D), the probability that a node has failed due to the failure of its dependent counterpart node in the other network. For network A, this probability is equal to the product of parameter r d and the probability that a counterpart node in B has failed: P(D) ¼ r d a B . In network B by analogy, this probability is equal to r d a A .
Writing equation (3) for both networks and inserting the results for P(I), P(E) and P(D) after summing over all k (and noting a A ¼ P k f k ð Þa A;k and a B ¼ P k f k ð Þa B;k ), we get a system of two coupled equations that describes the system of networks, Additional phase diagrams. Figure 8a shows the collection of stable solutions (solid blue lines) and unstable solutions (dashed red lines) for the activity z A ¼ 1 À a A of network A, with parameter values as used in Fig. 1a Fig. 1a. Green circles in this figure correspond to the stable states found in Fig. 1a and red crosses correspond to the unstable solutions for z A form Fig. 1a. Figure 8b shows an analogous phase diagram for the parameters with values as in Fig. 8a, again for a range of p Ã A .
Forbidden transitions. Transition lines for 12-21 and 21-12 do not appear in our phase diagram and it is quite easy to understand why. Let us assume that the transition line for 12-21 does exist. To obtain that transition, the idea would be to simultaneously increase p Ã A and decrease p Ã B (that is, increase the damage in one part of the system and decrease in another part). Suppose we are in phase 12 and infinitesimally close to the supposed transition line. Considering the local geometry of this line, we may be able to observe its angle with respect to the p Ã A axis. If a transition occurs when increasing p Ã A and decreasing p Ã B , the tangent on the supposed line would have an angle of y 2 0; p 2 Â Ã . From here, it follows that by increasing p Ã A only, while keeping p Ã B constant, we would also make a transition (cross the transition line). The only other possibility would be that we were moving along the transition line, but this is easy to disprove because it would imply that the transition does not depend on p Ã A . If increasing p Ã A only causes a transition, the transition must end in state 22 and not in 21. This is because if we only increase p Ã A , we increase damage to both network A (directly) and network B (indirectly, through the interdependent links).
Geometry of the Manhattan distance minimization problem. The optimal strategies shown in different colours in Fig. 5 are derived from the geometrical reasoning shown in Fig. 9. Figure 9a shows a plot of a series of curves consisting of points at identical Manhattan distances from point S 1 (equidistant curves). They produce a 'diamond' shape, and the minimal Manhattan distance between point S 1 and the green region translates into the task of 'fitting' the diamond so that it just touches the green region and its centre is at S 1 . The diamond in Fig. 9a touches the green region at two points-triple points, which are the solution to the minimization problem. Figure 9b shows the solution for point S 6 in the light blue region. Here the solution suggests a different strategy-decreasing only p Ã A .
Credit default swaps. Figure 6b shows an analysis of 5-year sovereign debt CDSs for a set of European countries: France, Germany, Italy, Spain, Portugal, Belgium, Austria, Denmark, Sweden, Greece, Ukraine, Hungary, Poland, Croatia, Slovenia, Romania, Bulgaria and Slovakia. This is the set of European countries that had a sovereign debt CDS in 2005. The set of Latin American countries we analyse consists of Brazil, Colombia, Argentina, Mexico, Venezuela, Chile, Peru and Panama. A CDS is typically used to transfer the credit exposure of fixed income products from one party to another. The buyer of the CDS is then obligated to make periodic payments to the seller of the CDS until the swap contract matures. In return, the seller of the CDS agrees to compensate (pay off) the seller who holds this third party debt if this (third party) defaults on the issued debt. A CDS is, in effect, an insurance against non-payment of a debt owed by a third party. The buyer of a CDS does not have to hold the debt of the third party but can speculate on the possibility that the third party will indeed default and the buyer can purchase the CDS for this speculative purpose. CDSs were developed in the 1990s and, given their simple structure and flexible conditions, they are now a major part of the credit derivative activity in the OTC market used to hedge credit risk. One of the most important aspects of a CDS is the definition of the 'credit event' that triggers the CDS. These events include bankruptcy, obligation acceleration, obligation default, failure to pay, repudiation (moratorium) and restructuring. In the case of the sovereign bond market, the last three are typically included in the contracts. CDSs are used by investors to hedge exposure to a fixed income instrument, to speculate on likelihood of a third party (reference asset) default, or to invest in foreign country credit without currency exposure.