From a single network to a network of networks

Network science has attracted much attention in recent years due to its interdisciplinary applications. We witnessed the revolution of network science in 1998 and 1999 started with small-world and scale-free networks having now thousands of high-profile publications, and it seems that since 2010 studies of ‘network of networks’ (NON), sometimes called multilayer networks or multiplex, have attracted more and more attention. The analytic framework for NON yields a novel percolation law for n interdependent networks that shows that percolation theory of single networks studied extensively in physics and mathematics in the last 50 years is a specific limit of the rich and very different general case of n coupled networks. Since then, properties and dynamics of interdependent and interconnected networks have been studied extensively, and scientists are finding many interesting results and discovering many surprising phenomena. Because most natural and engineered systems are composed of multiple subsystems and layers of connectivity, it is important to consider these features in order to improve our understanding of such complex systems. Now the study of NON has become one of the important directions in network science. In this paper, we review recent studies on the new emerging area—NON. Due to the fast growth of this field, there are many definitions of different types of NON, such as interdependent networks, interconnected networks, multilayered networks, multiplex networks and many others. There exist many datasets that can be represented as NON, such as network of different transportation networks including flight networks, railway networks and road networks, network of ecological networks including species interacting networks and food webs, network of biological networks including gene regulation network, metabolic network and protein–protein interacting network, network of social networks and so on. Among them, many interdependent networks including critical infrastructures are embedded in space, introducing spatial constraints. Thus, we also review the progress on study of spatially embedded networks. As a result of spatial constraints, such interdependent networks exhibit extreme vulnerabilities compared with their non-embedded counterparts. Such studies help us to understand, realize and hopefully mitigate the increasing risk in NON.


INTRODUCTION
The interdisciplinary field of networks has attracted enormous attention in the past 15 years due to its widespread appeal and applicability to many different complex systems in a large variety of disciplines [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Many real-world complex networks are not homogeneous and may obey a power-law degree distribution, and such networks are called scale-free (SF) networks [2,6]. There are plenty of real-world examples of SF networks, such as the Internet [20] representing the network of computers, the WWW [3] representing the network of hypertext documents, social networks [21][22][23][24][25][26] representing the relations between individuals, infrastructure networks [27], and biological networks [28,29]. The discovery of SF networks had a large impact on the basic properties of networks, which exhibit drastically different properties than the classical Erdős-Rényi (ER) networks [30]. For example, SF networks are significantly more robust than ER networks to random failures but less robust to targeted attack [3][4][5]. Percolation theory, based on REVIEW Gao et al. 347 statistical physics, provides us the method to analyse the robustness of a network in terms of node and link failures. Using percolation theory, we can study the robustness of a network and predict the critical percolation threshold, i.e. the fraction of removed nodes (or links) that leads to the collapse of the network [4]. Furthermore, using percolation theory one can address some other issues, such as efficient attacks or immunization [5,8,31,32], for obtaining optimal path [33] as well as for designing robust networks [34]. When networks are embedded in space, one can characterize them by a dimension which depends on the link length distribution [12,35,36]. More details can be found in the section entitled 'DIMENSION OF SPATIALLY EMBEDDED NETWORKS'. However, almost all network research has been focused on the properties of a single network that does not interact or depend on other networks. In reality, diverse critical infrastructures are coupled together and depend on each other, including systems such as water and food supply, communications, fuel, financial transactions and power stations [37][38][39][40][41]. For example, the electric power network provides power for pumping and control systems of water network, the water network provides water for cooling and emissions reduction of power network, the fuel network provides fuel for generators for electric power network and the electric power network provides power to pump oil for fuel network, etc. (see Fig. 1). In interdependent networks, the failure of nodes in one network leads to the failure of dependent nodes in other networks, which in turn may cause further damage to the first network, leading to cascading failures and possible catastrophic consequences. For example, blackouts in various countries have been the result of cascading failures between interdependent communication network and power grid [42,43]. In 2010, Buldyrev et al. [44] developed an analytical framework to study percolation of two interdependent networks subject to cascading failures based on a generating function formalism [45,46]. They found a first-order discontinuous phase transition, which is dramatically different from the second-order continuous phase transition found in isolated networks. Parshani et al. [47] studied two partial interdependent networks and found that as the coupling strength decreases below a critical threshold, the percolation transition becomes second order. Because in real-world scenarios, the initial failure of important nodes ('hubs') may not be random but targeted, Huang et al. [48] proposed a mathematical framework for understanding the robustness of interdependent networks under target attack, which was later extended by Dong et al. [49] to targeted attack on partially interdependent networks. They [48,49] developed a general technique that uses the random-attack solution to map the targeted-attack problem in interdependent networks. Also in real-world scenarios, each node in one network may depend on a number of nodes in the other network, while Buldyrev et al. [44] assume that a node in one network depends on a node in the other network. Therefore, Shao et al. [50] proposed a theoretical framework for understanding the robustness of interdependent networks with a random number of support and dependence relationships. In many cases, when high-degree nodes in one network depend with a high probability on low-degree nodes of another network, the coupled system becomes vulnerable. In the case of interdependent flight network and railway network, a city has high degree in the flight network and with high probability has high degree in the railway network. To better understand this phenomenon, Parshani et al. [51] proposed an 'intersimilarity' measure between the interdependent nodes. They studied a system composed of the interdependent world-wide seaport network and the world-wide airport network, and found that as the interdependent networks become more intersimilar, the system becomes more robust. The case in which all pairs of interdependent nodes in both networks have the same degree was solved analytically by Buldyrev et al. [52]. The effect on percolation of intersimilarity between the coupled networks was studied analytically by Cellai et al. [53] and Hu et al. [54]. In interdependent networks [44], there are two types of links: (i) connectivity links within each network and (ii) dependence links between networks. The connectivity links connect the components within a network to perform its function and the dependence links represent the fact that a given node in one network may depend on a node in the other network [55]. Hu et al. [56] studied a more realistic coupled network system with both dependence and connectivity links between the coupled networks, and found a mixed first-order and secondorder hybrid transition. A percolation approach for a single network in the presence of dependence links was also developed [57][58][59]. More recently, Gao et al. developed an analytical framework to study the percolation of a tree-like network formed by n interdependent networks [60][61][62], which was later extended by Dong et al. [63] to the study of targeted attacks on high-degree nodes. Gao et al. [60] found that while for n = 1 the percolation transition is of second order, for any n > 1, cascading failures occur and the network collapses in a first-order transition. Dong et al. [64] found that the cascading failures during the abrupt collapse involve a spontaneous second-order percolation. Baxter et al. [65] studied the avalanche collapse of interdependent networks,  [69]. These complex relationships are characterized by multiple connections between infrastructures, feedback and feed-forward paths, and intricate, branching topologies. The connections create an intricate web that, depending on the characteristics of its linkages, can transmit shocks throughout broad swaths of an economy and across multiple infrastructures. It is clearly impossible to adequately analyse or understand the behavior of a given infrastructure in isolation from the environment or other infrastructures. Rather, one must consider multiple interconnected infrastructures and their interdependences. For example, the reliable operation of modern infrastructures depends on computerized control systems, from supervisory control and data acquisition (SCADA) systems that control electric power grids to computerized systems that manage the flow of railcars and goods in the rail industry. In these cases, the infrastructures require information transmitted and delivered by the communication infrastructure [69].
which is called multiplex networks. The multiplex networks can be considered as n interdependent networks without feedback condition, i.e. tree-like network of networks (NON), where equation (2) in [65] is another description of equation (30) in [61]. And there is another extension which can be considered as the percolation of tree-like NON in [66]. More recently, Gao et al. also developed a general framework to study the percolation of any 'network of networks' [67], which can be used to study the robustness of any NON. More details can be found in the section entitled 'FRAMEWORK FOR A NET-WORK OF INTERDEPENDENT NETWORKS'.
Cascading failures have caused blackouts in interdependent communication and power grid systems spanning several countries [42]. To be able to design resilient infrastructures or improve existing infrastructures, one needs to understand how vulnerability is affected by such interdependences . Due to the broad range of applications of NON, there are many studies in this field [44,62,70,71], such as the robustness of network of ecological networks [72], robustness of multiplex networks [65], percolation of two antagonistic networks [73], transportation on two interconnected networks [74], eigenvalue and eigenvector REVIEW Gao et al. 349 analysis of multilayered networks [75][76][77][78] and many more. Furthermore, the idea of NON has been applied to some other models such as evolutionary game [79], synchronization [80] and epidemic models [81]. It was found recently that interdependent spatially embedded networks are extremely vulnerable to various types of attacks [82][83][84], see also the section entitled 'DIMENSION OF SPATIALLY EMBEDDED NETWORKS'.
In this paper, we discuss the percolation of a network formed by interdependent networks in the section entitled 'Framework for a network of interdependent networks'. We review the new phenomena found in interdependent networks and show that the new percolation laws of interdependent networks generalize the known results for single networks. Since spatially embedded networks are common in our world, we also discuss results on properties of spatially embedded networks in the section entitled 'Dimension of spatially embedded networks'. We review results on the extreme vulnerability of spatially embedded networks in the section entitled 'Extreme vulnerability of interdependent spatially embedded networks'.

FRAMEWORK OF TWO PARTIALLY INTERDEPENDENT NETWORKS
Consider a system composed by two networks A and B with the number of nodes N A and N B , respectively. The nodes within network A (B) are randomly connected by A (B) edges with degree distribution P A (k) (P B (k)), and a fraction q A (q B ) of network A (B) nodes depends on the nodes in network B (A) with no-feedback condition [62]. The no-feedback condition means that a node from one network depends on no more than one node from the other network, and if node A i depends on node B j , and B j depends on node A k , then k = i. We suppose that the initial random removal of nodes from network A is a fraction 1 − p. This removal causes further nodes in A that become disconnected to its giant component to also fail. Thus, the failed nodes in network A cause the dependent nodes in network B to fail and the nodes in network B that do not belong to the giant component also to fail, which causes further failures to network A. This iterates back and forth (cascading failure) until there is no further damage from one network on the other. At the steady state we have [57] where x and y satisfy and g A (x) and g B (x) are the generating functions of networks A and B, respectively [44,45,62,85].

FRAMEWORK FOR A NETWORK OF INTERDEPENDENT NETWORKS
In many real systems, there are more than two interdependent networks, for example, the biological networks, infrastructures and transportations. Understanding the way system robustness is affected by such interdependences is one of the major challenges when designing resilient infrastructures.
Here, we review the generalization of the theory of two interdependent networks [44,47] to a system of n interdependent networks [60][61][62], which can be regarded as an NON.

General framework
Consider a system (NON) composed of n nodes, and each node in the NON is a network itself and each link represents a fully or partially dependent pair of networks. We assume that each network i (i = 1, 2, ... , n) of the NON consists of N i nodes linked together by connectivity links and is characterized by a degree distribution of links P i (k). Two networks i and j form a partially dependent pair if a certain fraction q ji > 0 of nodes of network i directly depends on nodes of network j, i.e. they cannot function if the nodes in network j on which they depend do not function. Dependent pairs are connected by unidirectional dependence links pointing from network j to network i. We assume that after an attack or failure, only a fraction of nodes p i in network i will remain. We also assume that only nodes that belong to a giant connected component of each network i will remain functional. This assumption helps to explain the cascade of failures: nodes in network i that do not belong to its  The results obtained using Equation (6) for ER networks and Equation (7) for RR networks agree well with simulations. After [61]. where x i satisfy the system of n equations, and the product is taken over the Ki networks interlinked with network i (i.e. Ki neighboring networks of network i) by the partial dependence links (nonzero q ji ). The term has the meaning of the fraction of nodes in network j which survive after the damage from all the networks connected to network j except that of network i is taken into account. Equation (5) represents the nofeedback condition, and for the case of the absence of no-feedback condition, y ij = x j .

Percolation laws
In this section, we summarize the percolation laws regarding the NON, which has been studied in recent works [60][61][62]67].
(a) No-feedback condition (i) For a tree-like NON formed by n ER [30] networks with the same average degreesk, p 1 = p, p i = 1 for i = 1 and q ij = 1 (fully interdependent), the size of the mutual giant component [60] for all networks follows Note that for n = 1, Equation (6) reduces to the known single ER network result [30], with a secondorder phase transition, while for n > 1, the system undergoes a first-order phase transition, as seen in Fig. 2a.
(ii) For a tree-like NON formed by n random regular (RR) networks [61] with the same degree k of each network, the size of the mutual giant component for all networks follows Here also for n = 1, Equation (7) reduces to the single RR result, as seen in Fig. 2b.
(iii) For an NON composed of n ER networks with the same average degreek, where each network depends on exactly m other networks (RR of ER networks), and each network depends on a fraction q of nodes in its dependent network, the size of the giant component [60] for each network follows Here, P ∞ does not depend on n (see Fig. 3), but when q = 0 or m = 0, we obtain P ∞ = p(1 − ek P ∞ ) which is the single ER result [30]. We assume that the dependency loops in the NON form a mismatch and therefore, for q = 1, removal of one node may cause a complete collapse of all networks [61]. If the dependency loops do not form a mismatch but form periodical circles, the percolation results are identical to those of tree-like NON (see Equation (6) The absence of first-order regime in NON formed of ER networks is due to the fact that at the initial stage, nodes in each network are interdependent of isolated nodes (or clusters) in the other network. However, if only nodes in the giant components of both networks are interdependent, all three regimes, second order, first order and collapse, will occur, like in the case of RR NON formed of RR networks. After [67].
for q = 1), since when a node is removed, identical nodes will fail in both cases [61].
(b) Absence of no-feedback condition (i) For an NON composed of n ER networks with the same average degreek, where each network depends on exactly m other networks (RR of ER networks), the size of the giant component [60] for all networks follows (ii) For an NON composed of n RR networks with the same degree k, where each network depends on exactly m other networks (RR of RR networks), the size of the giant component [60] for all networks follows Here, again when m = 0 or q = 0, Equations (9) and (10) reduce to the single network result and the numerical solutions of (9) and (10) as well as simulation results are shown in Figs 4 and 5. This rich generalization shows that the percolation on a single network studied for more than 50 years is a limited case of the more general case of NON.

DIMENSION OF SPATIALLY EMBEDDED NETWORKS
While topological properties of non-embedded networks have been extensively studied, many networks including transportation networks, power grid and even social networks are embedded in Euclidean space with spatial constraints. Furthermore, many realistic spatial networks including the Internet and airline network are found to have a broad distribution of link length, P(r) ∼ r −δ [86].
To characterize the structure and basic physical properties of these spatial networks, the fundamental concept of dimension and its measurement algorithm have been proposed by Li et al. [8,12]. The results indicate that networks characterized by a broad distribution of link lengths have a dimension d defined as M ∼ r d , where r is the Euclidean distance, which is higher than that of the embedding space and the dimension d depends on δ. In systems embedded in 2D space, for δ < 2, the system has an infinite dimension d, while for δ > 4, it has dimension d = 2. For 2 < δ < 4, the dimension d varies continuously from d = ∞ to d = 2 [12].
To illustrate the origin of dimension of spatial networks, we analyse how both the mass M and the Euclidean distance r increase with the topological distance l (minimum number of links between two sites in the network), from which we find the relation M ∼ r d that determines the dimension d [8,35]. Numerical results, based on extensive simulations, suggest that for networks embedded in 2D lattice, in the intermediate regime 2 < δ < 4, M(l) and r(l) increase with l as a stretched exponential. However, Figure 5. The giant component for an RR NON formed of RR networks with feedback condition, P ∞ , as a function of p for RR of degree k = 6 and m = 3, for two different values of q. The curves are obtained using Equation (10), which shows a first-order phase transition when q is large but a secondorder phase transition when q is small. After [67].  Figure 6. The size of the giant component P ∞ at steady state after random failure of a fraction 1 − p of the nodes of (a) two interdependent lattice networks with periodic boundary conditions and (b) two RR networks. All networks are of size 16 × 10 6 nodes and the same degree distribution P(k) = δ k, 4 . The coupling between the lattices and between the RR networks changes from q = 0 to 0.8 with step 0.1 (from left to right). In the case of interdependent lattices, only for q = 0 (no coupling, i.e. a single lattice) the transition is the conventional second-order percolation, while for any q > 0 the collapse is abrupt in the form of first-order transition. This is in marked contrast to the case of interdependent RR networks, where only for q > q c ∼ =0.43, the transition is abrupt, while for q < q c , the transition is continuous. After [83].
interestingly, the relation between M and r is found to be a power law which determines the dimension d of the embedded system.
The dimension found in this way is a key concept to understand and characterize not only the spatial network itself, but also dynamical processes on spatial networks, such as diffusion reaction [36] and critical phenomena including percolation [87]. For δ > 4, network dimension is similar to dimension of embedding lattice, and percolation belongs to the universality class of percolation in regular lattices. For δ < 2, the percolation transition belongs to the universality class of percolation in ER networks (mean field). In the intermediate range 2 < δ < 4, the percolation properties show new behavior different from mean field, with critical exponents that depend on network dimension. The problem of percolation of interdependent networks, where each (or one of them) has a power-law distribution of link length, is still an open question.
Besides percolation, dimension of spatial network is also found to influence the diffusion processes with applications ranging from epidemics and transport properties of materials to financial activities. Also, chemical reaction processes on the spatially constrained networks are found to be controlled by the system dimension d [36].

EXTREME VULNERABILITY OF INTERDEPENDENT SPATIALLY EMBEDDED NETWORKS
There exists a critical coupling of interdependent nodes in interdependent non-embedded networks, above which single failure can invoke cascading failures that may abruptly break down the system, while below this critical coupling, the transition is continuous and small failures lead only to small system damage but not to a collapse. As a result of this critical coupling, resilience of these interdependent networks is threatened by the risky abrupt collapse phenomenon. Therefore, when the coupling of interdependent nodes in NON is below this critical coupling, the system can be considered in a safe region free from the risk of abrupt collapse.
Many critical infrastructure networks are embedded in space. These include the transmitters of the telecommunications networks which are functioning with electricity from power grid, which is in turn controlled by the telecommunications network. Interdependent spatially embedded networks, modeled by coupled lattices, have been found significantly more vulnerable compared to non-embedded networks [83] as seen in Fig. 6. In contrast to nonembedded networks, there is no critical coupling of interdependence, but any small coupling of interdependent nodes will lead to an abrupt collapse of first-order transition. In such systems there is no safe region. Furthermore, when the length of interdependence links between sites of two lattices is constrained by distance r, as seen in Fig. 7, Li et al. [82] suggest that for full coupling (q = 1), the transition changes from first order to second order at r = r max = 8: percolation for small r, r < r max , is a secondorder transition, and for large r, r > r max , it is a Figure 7. Two square lattices A and B, where in each lattice every node has two types of links: connectivity links and dependence links. Every node is initially connected to its four nearest neighbors within the same lattice via connectivity links. Also, each node A i in lattice A depends on one and only one node B j in lattice B via a dependence link (and vice versa), with the only constraint that |x i − x j | ≤ r and |y i − y j | ≤ r. If node A i fails, then node B j fails. If node B j fails, then node A i fails. Network A is shifted vertically for clarity. After [82]. Gao et al. 353 first-order transition. When the coupling is reduced (q < 1), r max increases as shown by Danziger et al. [88]. Spatial constraints not only manifest themselves in network vulnerability to random failures, but also make the system even more susceptible to spatially localized failures caused by natural disasters or malicious attacks. Local failures (like the Tohoku earthquake and tsunami in 2011 or the recent typhoon in the Philippines or chemical/biological attack) can cause total failures within a given radius from the perturbation, which may propagate the entire interdependent network system.

REVIEW
Under random attack, a finite coupling of the system has to be removed to trigger the collapse of system. Surprisingly, it is found [84] that a localized failure in interdependent spatial networks can invoke substantially more damages to the networks than an equivalent random failure in many scenarios. Especially, if the localized failure size is larger than a critical size (but still a zero fraction), it will spread globally and fragment the entire system.

Evolutionary game
Traditional spatial evolution game on a single network assumes that players only interact with their direct neighbors, through ignoring real situations that networks are interdependent and interacting with each other. In fact, the interactions between networks have important influence on the results of evolution game. A recent article [89] studies the evolutionary dynamics of the prisoner's dilemma game and finds that the cooperative behaviors are enhanced by the interdependence between networks. Furthermore, the impact of the interdependence strength on the evolution game has been also analysed. Wang et al. [90] find that there exists an intermediate region for the fraction of interdependent links between networks leading to optimal result of the evolution game. Jiang and Perc [91] consider interdependent networks with different topology and strategy updating rule and reveal that there always exists a fraction of interdependent links leading to the optimal result. Furthermore, the impact of other properties of interdependences on the evolution game has been studied. An unbiased coupling has been proposed in [92] to investigate the condition of the spontaneous emergence of interdependent network reciprocity. Tang et al. [93] study the game of public goods on two interdependent networks and find that the more robust the links between cooperators, the more prevalent is the cooperation. The symmetry of cooperation has been stud-ied in [79], and a threshold value of interdependent fraction has been found.

Synchronization
Synchronization on interdependent networks is a recently developed research area. The impact of the coupling strength between interdependent networks on the synchronization of the system has been studied in [80]. For a single 1D network, no synchronized phase exists at any finite coupling strength, while for a single Watts-Strogatz small-world network, there is a mean-field transition. Um et al. [80] find that when these two kind of networks are coupled, the relation between synchronization and the coupling strength depends on the intranetwork coupling strength of 1D network. The synchronization dynamics as a function of the couplings and communication lag have been analysed in [94], and a new breathing synchronization regime has been found.

Epidemics on interconnected networks
Disease transmission is one of the typical dynamical processes in complex networks. Susceptibleinfected-susceptible (SIS) model and susceptibleinfected-recovered (SIR) model are usually used to describe the disease transmission. Recently, it has been found that the disease spreads not only in a single network, but also in the interconnected network. Saumell-Mendiola et al. [95] have developed a method to calculate the epidemic threshold of the SIS model in the interconnected networks, based on the heterogeneous mean-field approach. They also find that the global epidemic threshold of interconnected networks is smaller than that of each component network separately. Wang et al. [96] have developed another method to calculate the epidemic threshold of the SIS model, using the largest eigenvalue of the interacting matrix of the interconnected networks. The SIR model in interconnected networks has also been studied [81]. It has been shown that for strong interconnected networks, the endemic state occurs or dies out simultaneously in the component network, but in weak interconnected networks, the endemic state may occur only in one network, without presence in other networks.

SUMMARY AND REMAINING CHALLENGES
We have reviewed recent studies on the robustness of a system of interdependent non-embedded networks and spatially embedded networks. In interacting networks, when a node in one network fails, it REVIEW usually causes dependent nodes in other networks to fail which, in turn, may cause further damage in the first network and result in a cascade of failures with catastrophic consequences.
Thus far, only a very few real-world interdependent systems have been analysed using the percolation approach [51,97,98]. We expect the framework of interdependent NON to provide insights leading to further analysis of real data on interdependent networks. Further studies of interdependent networks should focus on (i) an analysis of real data from many different interdependent systems and (ii) the development of mathematical tools for studying real-world interdependent systems.
Many real networks are embedded in space, and the spatial constraints strongly affect their properties [12]. The embedded interdependent networks are extremely vulnerable [83,84]. It is important to develop algorithms for improving the robustness of embedded interdependent networks according to their topological structure and function. Currently, three methods have been proposed to improve the robustness of the non-embedded interdependent networks: (i) increase the fraction of autonomous nodes [47], especially nodes with high degree; (ii) design the dependence links such that they connect the nodes with similar degrees [51,52]; and (iii) protect the high-degree nodes against attack [99]. For a recent study on spontaneous recovery and breakdown of networks, see Majdandzic et al. [100].

FUNDING
This work was supported by the National Natural Science Foundation of China (61104144 and 61374160). Shlomo Havlin acknowledges the European LINC project (289447) funded by the EC's Marie Curie ITN program and MULTIPLEX (EU-FET project 317532) projects, the Deutsche Forschungsgemeinschaft (DFG), the Israel Science Foundation, BSF, ONR and DTRA for financial support.