Extreme risk induced by communities in interdependent networks

Networks in nature not only depend on each other but also have internal community structures, such as infrastructure networks with links within and across geographic regions. The communities play an important role when the networks undergo localized failures in specific regions, for instance when natural disasters or economic sanctions disrupt a local community region and consequently influence the whole system. How a disruption in one community propagates through the entire system is a crucial, but still open, question. Here we find that the community structure embeds extreme risk: weakening the community strength could abruptly drive the system to a precarious state. Examining the business-flight network among cities as a proxy for the world economy, we find this real coupled system evolving towards the extreme vulnerable phase due to ongoing globalization. This shows the community risk indeed exists in real world networks and deserves more attention from the scientific community. The effect of localized disruption or failure in interdependent networks with internal community structure remains an open question. Adopting the generating function approach, the authors are able to uncover rich phase transition behaviours and the associated risks for such system, and by studying real networks under random failures within a community, they find that weakening the community strength could rapidly drive the system to a precarious state.

I nterdependent networks are important since many real work systems have interactions at different levels, and such complexity often leads to new and rich phase transition behaviors not present in single layer networks [1][2][3] . Although much research has focused on the robustness to failure in multilayer interdependent networks 1,[4][5][6][7][8][9][10][11][12] , these studies have assumed that networks are unstructured. This is in sharp contrast to real-world networks that have an internal community structure 6,[13][14][15][16] . For example, transportation networks have more connections within urban regions than between urban regions. A group of countries usually have more economic ties within the group than with countries outside the group. A problem of particular interest is determining how these complex systems with a rich community structure behave under localized disruption. Such community structure has been studied on single layer networks 17 to understand the impact of removing intercommunity links, and rich phase transition has been found when the number of communities changes. Random attack on intercommunity links on interdependent networks have also been found 18,19 to have even more complex transitions and scaling relations. Yet, certain attacks are localized in specific communities, like natural disasters and economic sanctions. The underlying mechanism controlling how such localized disruption in one community of a complex system disseminates throughout the entire system has not been understood, as well as the conditions that the network becomes vulnerable to abrupt collapse. Analytically, it is challenging to study since the network system with both interdependency and communities structure is highly complex.
Here, we present a generalized framework of interdependent multilayer network with community structure 1,6 based on generating functions, and study cascading processes that occur across the entire network initiated by random damage disruption in a single community. Such framework can be generalized to any layers of interdependent networks with arbitrary number of communities. Our analysis reveals that such system has rich phase transition behaviors that are much more complicated than interdependent networks without community structures. In particular, we find both theoretically and empirically that network robustness changes abruptly from safe to vulnerable as the strength of the community changes. The safe region of the system is characterized by the lack of phase transition phenomenon, i.e. the system does not disintegrate even with one whole community fully removed from the system. The vulnerable region is when phase transition is present, i.e. removing enough fraction of one community will disintegrate the whole system. Employing the business-flight network among cities of North America, Asia, and Europe as an example with strong community structure, we confirmed the presence of such risk in this interdependent system. More strikingly, this system is evolving towards the transition point from the safe to the vulnerable phase due to ongoing globalization.

Results
Theoretical analysis. We consider an interdependent multilayer network with a community structure. In a simple network consisting of two network layers A and B, each layer has the same number of communities m in which A i and B i are of size N i , i = 1,2, …...m. Every node in A i has exactly one interdependent node in B i , and vice versa [see Fig. 1a]. We define the generating function of A i as where p iA (k) is the degree distribution of nodes in community i of layer A. We define the community structure using the distribution of both intra-and intercommunity links 14,16 . For a community A i within network layer A and has average degree 〈K iA 〉, we use 〈k ijA 〉 to denote the average number of intercommunity links per node in community A i that connects to community A j . We further let α ij; = 〈k ijA 〉/〈K iA 〉 to be the fraction cross community links in A i . These parameters can also be defined in terms of stub that is a link with one end from a node and the other end not connected yet. Then in the network formation process one make α ii fraction of stubs in community A i to be connected to other stubs in community A i , and connect α ij fraction of stubs in community A i to stubs in community A j . β ii and β ij are defined similarly for community B. The generating functions for communities A i and B i are G iA ðξ 1 ; ξ 2 ; :::; ξ m Þ ¼ P P iA ðk 1 ; k 2 ; :::; k m Þξ k 1 1 :::ξ k m m and G iB ðζ 1 ; ζ 2 ; :::; ζ m Þ ¼ P P iB ðk 1 ; k 2 ; :::; k m Þζ k 1 1 :::ζ k m m , respectively 20 . P iA (k 1 , k 2 , ..., k m ) is the probability of finding a node in community A i with k 1 links connecting to nodes in community 1, k 2 links connecting to nodes in community 2,…, and k m links connecting nodes to community m. P iB (k 1 , k 2 , ..., k m ) is defined analogously for network layer B. These generating functions can be expressed by substituting x in Eq. (1) with (ξ 1 , ξ 2 , ..., ξ m ) (derivation provided in Methods), G iA ðξ 1 ; ξ 2 ; :::; ξ m Þ ¼ Q iA P m j¼1 α ij ξ j ! ; G iB ðζ 1 ; ζ 2 ; :::; We next attack community A in one of the network layers (say A i ) by initially removing a fraction of 1−p i nodes in community A i . In interdependent networks, it is usually assumed, based on percolation theory, that nodes become non-functional if they become disconnected to the network giant component 1,5,[21][22][23][24][25] . That is all nodes outside the giant component in network A are removed, and so are their interdependent nodes in network B. This cascading process continues until no more nodes can be removed from the system. We finally obtain the mutually connected giant component of the remaining functional nodes.
As we remove nodes and links during the cascading process, the entire network breaks down into the components connected through inter-and intra-community links. The components in community A i comprises the set of A i nodes belonging to a single percolation component of the entire network. The largest of these components is the giant component of A i , and only nodes of the giant component continue to function. Figure 1b-d shows the various types of giant component in multicommunity networks. Note that some nodes are not linked to nodes in the same community but to nodes in the other community.
At the end of the cascading process, the remaining size of A i and B i is where :::; f 1 m Þ are the giant component sizes of community i in network layers A and B, respectively. Here f 1 i is the probability that a node in A i along a randomly selected link is non-functional (i.e. not connected to the giant component at the steady state). It satisfies the self-consistent equation (see Methods for details).

G iB
1 ðζ 1 ; ζ 2 ; :::; ζ m Þ is defined analogously for community B i . For example, for a pair of Erdös-Renyi networks 26 , Eq. (3) reduces to a simple form, the derivation of which is provided in Supplementary Note 1. The internal structure of our network model is a generalization of a model 1 in which networks have no internal community structure.
To demonstrate critical phenomena in our network model, we consider a simple case of two equally sized communities that are symmetrical in each layer of the network 27 , i.e., m = 2. Without loss of generality, we attack community 1 at the initial stage, i.e., p 1 < 1 and p 2 = 1. We also set α 11 = α 22 , Here the set of parameters α 11 and β 11 is sufficient to describe the community structure. A key parameter that quantifies the robustness of the system is the size of the critical value p 1c that describes the threshold below which the entire system disintegrates with no remaining functional giant component. Thus the smaller the p 1c value, the less vulnerable the network, implying that when p 1c = 0 the network is perfectly robust.
Our first main finding is that a stronger community structure does not always increase the robustness of the interdependent networks, a phenomenon significantly different from the one found in single networks. In single networks, the stronger community structure always increases robustness as shown in Fig. 2a (see Supplementary Note 2 for derivation). Differently, Fig. 2b shows the critical point p 1c against the community strength α 11 of network layer A when the community strength β 11 of network layer B is fixed. The average degree is fixed at 〈K A 〉 = 〈K B 〉= 4 in both network layers. When the community structure is strong in layer B-which occurs when β 11 > 0.436-there is a monotonic dependence of p 1c on α 11 , but when the community structure in B is relatively weak (i.e., when β 11 < 0.436) the behavior of p 1c is non-monotonic.
In order to better understand this, we draw the contour of p 1c with respect to the changes of community structures strength in Fig. 2c. Figure 2c shows that when β 11 < 0.436 the contour of p 1c is the bulging of the equipotential line, specified by the constant value of p 1c , cuts the horizontal line characterized by constant β 11 value twice, resulting in non-monotonous changes as the parameter α 11 changes. Besides this, we find that that system robustness falls into one of the three regions (see Fig. 2c): "vulnerable" p 1c >p r 1c À Á , "robust" p 1c <p r 1c À Á , and "safe" (p 1c = 0, which means that even after the removal of all nodes in a community, there still exist a giant component of the system) regions with the increase of community structure strength in both two network layers A and B. Here p r 1c is the critical point of the corresponding interdependent networks without any community structures (α 11 = β 11 = 0.5).
We also note that the contour lines intersect at two interesting symmetrical points: (α 11 , β 11 ) equals to (0.436, 1) and (1, 0.436) (see Fig. 2c and Supplementary Note 3 for deviation), which implies p 1c changes abruptly at these two unusual points. α 11 = 1 or β 11 = 1 means network A or B has two disconnected (localized) communities. This is a good approximation of real-world network segmentations which could be either geographical or political/economic imposed by embargo/sanctions. For instance, the world-wide business-flight network, which we examine below, has β 11 = 0.98. We find that when one network layer has disconnected communities, as the community strength in another network layer weakens, the critical point suddenly jumps from 0 to a finite number, which is our second main finding. This abrupt jump of the critical point value is a first-order phase transition in which a small change in community strength dramatically increases structural risk (see Fig. 2d and Supplementary Note 4 for derivation).
Empirical implications. An example of a multilayered community structure network is the network of global cities. We examine the data from the system of North American, European, and Asian cities-three different communities-in which transportation and business connections among them define two layers of network in the system. As expected, there are more connections among cities located on the same continent than that among cities located on different continents 28 . We collect business and flight data for 145 North American cities, 158 Asian cities, and 334 European cities and for companies across 21 major industrial sectors in 2010 and 2013 (see Supplementary Note 10 for data description) 29,30 . We use the data to construct an interdependent network of business and flight connections among these cities, where cities are nodes and business and flight connections are links between the cities. In the business network layer (Fig. 3a), a connectivity link between two cities is formed when at least 10 pairs of companies have business connections with each other. In the flight network layer (Fig. 3b), a link is formed when there are at least 200,000 passenger trips between the two cities annually. Simulations demonstrate that the results are not sensitive to these two threshold values (see Supplementary Fig. 5).
The business and transportation network layers are interdependent because businessmen must travel to conduct their business, and airport of one city also dependents on the companies of this city. Usually, it is not easy to obtain the interdependency relationships on important infrastructure networks, thus we begin by assuming a single interdependency link 1 between business and transportation in the same city. We then assume the typical failure mechanism for interdependent networks in which the failure of one node leads to the failure of its interdependent node.
In numerical simulations, we first construct three interdependent networks comprising only cities in North America vs Asia, Europe vs Asia, and Europe vs North America. Secondly, altering the business network A but fixing the flight network B (β 11 is close to 1), we tune down community strength α 11 of business network by rewiring some of its links while keeping their node degrees unchanged (see Supplementary Note 5-8 for the technical details). This approach enables us to assess the system robustness, i.e. how p 1c depends on α 11 , and for what range of parameter values does the system reside in the safe phase or the vulnerable phase. As shown in Fig. 3c-e, we see a phenomena similar to Fig. 2d that p 1c changes abruptly from 0 to a finite value, signaling the emergence of the community structure risk in the system. As the trend of globalization continues, we expect the community structure in the business network to become weaker. Indeed this was the case when comparing the value of community strength α 11 value between year 2010 and 2013 in Fig. 3 for each pair of communities North America-Asia, Europe-Asia, and Europe-North America. Our results from Fig. 3c-e show that as globalization drives global economic system to be more integrated, the North American and European economies exhibit larger robustness than the Asian economy. Finally, we construct an interdependent network of business and flight connections comprising all cities in North America, Europe, and Asia and we find that the overall system is still shifting towards the abrupt change through removing nodes within the Europe community (Fig. 3f). Our results and analysis are based on communities divided according to different continents, and use them to serve as proxies to represent the coarse structure of the world economy. Using our mathematical framework, we can estimate the risk towards the unstable region given the trend from 2010 to 2013. In essence, the world economy, according to this mechanism, had been likely heading towards the unstable phase. In Supplementary Note 9 we discuss some additional weaker forms of interdependency between two network layers characterized by probabilistic interdependency links.

Discussion
In summary, we have studied the robustness of multilayer interdependent networks under community attacks, characteristic of many real system constrained by geographic locations. Through our developed framework, we find that such a system attack is plotted against its community strength α 11 . As the positive community structure strengthens with increasing α 11 , p 1c decreases monotonically before reaching 0, after which point the two communities both can have a giant component even when the other community is removed from the network. b Critical transition point p 1c of interdependent networks subjecting to community attack is plotted against the community strength α 11 of layer A for different β 11 values of layer B. For large β 11 values (0.5 and 0.8), the behavior is similar to single networks with monotonic decrease of critical threshold p 1c value. Whereas for small β 11 (0 or 0.4), the change of p 1c is non-monotonic and does not reach the safe state of p 1c = 0. The solid red circle represents random failures of an unstructured random network system which is equivalent to α 11 = β 11 = 0.5. c Phase diagram of p 1c with respect to the changes in community structures described by α 11 and β 11 in interdependent networks. The upper right region is the safe region, in which even if we remove a whole community, there still exist a giant component in the whole system. The robust region between boundary p 1c = 0 and p 1c ¼ p r 1c represents a phase that is more robust than unstructured networks without community structures. The vulnerable region below the boundary p 1c ¼ p r 1c represents a phase that is more vulnerable than unstructured networks. In b and c, the two interdependent networks are Erdös-Rényi (ER) networks and both have two communities with equal size and equal average degree 4. d When one layer of the networks is completely localized (β 11 = 1), extreme community structure risk exists when the community strengths α 11 decreases, as p 1c changes discontinuously from 0 to 0.34. In the system, K = 4 for two equally sized communities. The network has a total of N ; = 10,000 nodes, with each community having 5000 nodes. Each data point is an average of 10 simulations and the standard deviation is presented by the error bar exhibits much more complex behaviors than the ones without clear community structures. Depending on the respective community strengths of the different network layers, the whole system may reside in one of the three different regions of phase diagram, with different phase transition behaviors. In particular, we find that the system embeds extreme risk when one of the layers has strong community structure. In contrary to a single network which is always more robust with enhanced community structure, interdependent networks showing strong community structure sometimes makes the system more vulnerable. In cases of extremely strong community structure with very sparse intercommunity in one layer of network, a small change in the other layer of network's community strength could induce abrupt changes in total systemic resilience if one community is under attack. This new finding adds to the growing knowledge of resilience of interdependent networks with community structure, in particular attacks on specific communities that is representative of realistic events like natural disaster and economic embargo. From the global business-flight networks taken as a proxy for the world economy, we observe that, as globalization weakens the community structure, the entire network is approaching a state in which there is a potential for abrupt disruptions in certain communities. This also leads to the finding that Asia strongly depends on North America and Europe economies but not the other way round. Our result is indicative for a broad range of systems which have community structures defined by geographic location and physical infrastructure. While we mainly focused on the analytical results of two communities on two-layer interdependent networks, our adopted mathematical framework allows for more complex network structures, and their phase transition behaviors could be very interesting for further studies. On the empirical side, our methods can be similarly extended to study other real interdependent networks with community structures, to give better understanding of their dynamical  Table 1 and Supplementary Tables 1 and 2. a Firm and b flight networks between Asia and North America. The red nodes represent cities in Asia, and the blue nodes represent the cities in the North America. The links represent air travel activities in the flight network and business activities in the firm network. We consider the Asian cities as one community and the North America cities as another community in both networks. c Empirical study of the structural risks shows the same abrupt change. The two communities are cities in North America and Asia, and the two layers of networks are business interaction networks and flight networks. The critical transition point p NA c (removing nodes within the North America community) is plotted against community structure strength α 11 of the business layer for the business-flight networks. The original value of α 11 is the empirical value of the system in 2010 and is marked on the horizontal axis (green dot), while the new value of α 11 is measured in 2013. Each data point is an average of 1000 simulations and the shaded region is the bounded by the standard deviation. The result shows that the whole system is shifting towards the abrupt change, where there will be extreme systemic risk. The decrease in α 11 value is possibly due to the effect of globalization, which is weakening the community structure. d, e Similar to c but the results is on the two communities of Europe vs. Asia and North America. From 2010 (green dot) to 2013 (purple dot), the system has shifted from the safe phase to the critical point p E c (removing nodes within the Europe community) close to the vulnerable state. f Similar trend for the systemic risks can be observed in an interdependent network involving all of three communities, and the node removal is done on the Europe community. For all of the simulations in c-f, we only show the attack on one community, as the attack on the other communities does not have the abrupt change on the systemic resilience, i.e. no p c is found.
behaviors for risk assessment, which then leads to the development of efficient risk mitigation strategies in those systems.

Methods
Multivariate generating function. In this subsection, we harness multivariate generating function to formalize the cascading process of interdependent networks with multiple communities 20,30,31 . We first construct the multivariate generating function of the interdependent networks with communities. For conciseness, we illustrate the generating function associated with the bi-community structure (a network with two communities). The result can be easily generalized to account for multicommunity systems. For bi-community system, the probability normalization requires that α 12 = 1 − α 11 and α 21 = 1 − α 22 and the inter-link identity requires that N 1 ⋅〈K 1A 〉⋅α 12 = N 2 ⋅〈K 2A 〉⋅α 21 . Therefore one can describe the structure of network A by one parameter α 11 , because the other three parameters can be calculated from the above three relations. By the definition of α 11 , the probability of finding a node in community 1 with intra-degree k 1 and inter-degree k 2 should be k 1 þ k 2 k 1 α k 1 11 ð1 À α 11 Þ k 2 . Such node with degree k 1 + k 2 has a prior probability p iA (k 1 + k 2 ). Thus one can write the generating function of community 1 in terms of ξ 1 and ξ 2 : The prefactor of ξ k 1 1 ξ k 2 2 represents the probability of finding a node in community 1 of network A with intra-degree k 1 and inter-degree k 2 . By a similar derivation, one can write out the generating function for community i in multicommunity network A and B as G iA ðξ 1 ; ξ 2 ; ::; Percolation equation. Here we show the percolation equation in the form of generating functions 32 . Initially a fraction (1 − p i ) of nodes from community i are randomly removed. The initial removal is followed by a cascading process of percolation failures of the rest nodes. Recall a single complex network with generating function G(x) and branching process G 1 (x) = G′(x)/G′(1) (ref. 33 ), it is known that after a random removal of 1 − p nodes, the generating function of the degree distribution of the remaining nodes can be written by G remain (x) = G(1 − p (1 − x)). At step n, say the size of the giant component of community i in layer A and B due to percolation failures are g n i and h n i . Analogously, we introduce u n i and v n i as the probability of finding a non-functional node by following a link of a randomly chosen node in community i of network A and B at step n after the initial removing of 1−p i fraction of nodes from each community, respectively. We define f n i by The physical meaning of f n i is the probability of finding a node along a randomly selected link in community i of the original network to be non-functional at step n. So the fraction of nodes that are still functional at step n + 1 is :::; f n m Þ; ð8Þ or equivalently g nþ1 i ¼ 1 À G iA ð1 À p 1 ð1 À u n 1 Þð1 À v n 1 Þ; 1 À p 2 ð1 À u n 2 Þð1 À v n 2 Þ; :::; 1 À p m ð1 À u n m Þð1 À v n m ÞÞ: ð9Þ u n i should obey a transcendental equation 3 u n i ¼ G iA 1 ð1 À p 1 ð1 À u n 1 Þð1 À v n 1 Þ; 1 À p 2 ð1 À u n 2 Þð1 À v n 2 Þ; :::; 1 À p m ð1 À u n m Þð1 À v n m ÞÞ; in which G iA 1 ðu n 1 ; u n 2 ; :::; u n m Þ is the generating function associated with the branching process: Similarly, the fraction of functional nodes h nþ1 i in community i in network B at step n ; + 1 is h nþ1 i ¼ 1 À G iB ð1 À p 1 ð1 À u n 1 Þð1 À v n 1 Þ; 1 À p 2 ð1 À u n 2 Þð1 À v n 2 Þ; :::; 1 À p m ð1 À u n m Þð1 À v n m ÞÞ; in which v n i is the probability of finding functional node by following a link of a randomly chosen node in community i of network B. It satisfies the following transcendental equation: v n i ¼ G iB 1 ð1 À p 1 ð1 À u n 1 Þð1 À v n 1 Þ; 1 À p 2 ð1 À u n 2 Þð1 À v n 2 Þ; :::; 1 À p m ð1 À u n m Þð1 À v n m ÞÞ; where G iB 1 ðv n 1 ; v n 2 ; :::; v n m Þ is the generating function related to the branching process: At equilibrium when no more nodes can be removed due to percolation and interdependency failures, the above important parameters u 1 i and v 1 i should obey the following self-consistent equations: À p 1 ð1 À u 1 1 Þð1 À v 1 1 Þ; 1 À p 2 ð1 À u 1 2 Þð1 À v 1 2 Þ; :::; 1 À p m ð1 À u 1 m Þð1 À v 1 m ÞÞ; v 1 i ¼ G iB 1 ð1 À p 1 ð1 À u 1 1 Þð1 À v 1 1 Þ; 1 À p 2 ð1 À u 1 2 Þð1 À v 1 2 Þ; :::; 1 À p m ð1 À u 1 m Þð1 À v 1 m ÞÞ; ð15Þ and the remaining size of i'th community μ 1 i becomes μ 1 i ¼ p i g 1 i h 1 i ¼ p i 1 À G iA 1 À p 1 ð1 À u 1 1 Þð1 À v 1 1 Þ; À Â 1 À p 2 ð1 À u 1 2 Þð1 À v 1 2 Þ; :::; 1 À p m ð1 À u 1 m Þð1 À v 1 m Þ ÁÃ Á 1 À G iB 1 À p 1 ð1 À u 1 1 Þð1 À v 1 1 Þ; À Â 1 À p 2 ð1 À u 1 2 Þð1 À v 1 2 Þ; :::; 1 À p m ð1 À u 1 m Þð1 À v 1 m Þ ÁÃ ;

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Received: 13 September 2018 Accepted: 11 March 2019 Community size, average degree, and community structure index are indicated α r ii ¼ K iA N i =ðK 1A N 1 þ K 2A N 2 Þ is the index if the given network has no community structure and α ii is the actual community structure index. β ii and β r ii are the same for network B