Published on

09-Apr-2017View

213Download

1

Transcript

PHYSICAL REVIEW E 87, 052804 (2013)Robustness of network of networks under targeted attackGaogao Dong,1,2,* Jianxi Gao,2,3 Ruijin Du,1,2,4 Lixin Tian,1 H. Eugene Stanley,2 and Shlomo Havlin51Nonlinear Scientific Research Center, Faculty of Science, Jiangsu University, Zhenjiang 212013, China2Center for Polymer Studies and Department of Physics, Boston University, Boston, Massachusetts 02215, USA3Department of Automation, Shanghai Jiao Tong University, Shanghai 200240, China4College of Mathematics Science, Chongqing Normal University, Chongqing 401331, China5Department of Physics, Bar-Ilan University, 52900 Ramat-Gan, Israel(Received 12 February 2013; published 16 May 2013)The robustness of a network of networks (NON) under random attack has been studied recently [Gao et al.,Phys. Rev. Lett. 107, 195701 (2011)]. Understanding how robust a NON is to targeted attacks is a major challengewhen designing resilient infrastructures. We address here the question how the robustness of a NON is affectedby targeted attack on high- or low-degree nodes. We introduce a targeted attack probability function that isdependent upon node degree and study the robustness of two types of NON under targeted attack: (i) a tree of nfully interdependent Erdos-Renyi or scale-free networks and (ii) a starlike network of n partially interdependentErdos-Renyi networks. For any tree of n fully interdependent Erdos-Renyi networks and scale-free networksunder targeted attack, we find that the network becomes significantly more vulnerable when nodes of higherdegree have higher probability to fail. When the probability that a node will fail is proportional to its degree,for a NON composed of Erdos-Renyi networks we find analytical solutions for the mutual giant component Pas a function of p, where 1 p is the initial fraction of failed nodes in each network. We also find analyticalsolutions for the critical fraction pc, which causes the fragmentation of the n interdependent networks, and for theminimum average degree kmin below which the NON will collapse even if only a single node fails. For a starlikeNON of n partially interdependent Erdos-Renyi networks under targeted attack, we find the critical couplingstrength qc for different n. When q > qc, the attacked system undergoes an abrupt first order type transition.When q qc, the system displays a smooth second order percolation transition. We also evaluate how the centralnetwork becomes more vulnerable as the number of networks with the same coupling strength q increases. Thelimit of q = 0 represents no dependency, and the results are consistent with the classical percolation theory of asingle network under targeted attack.DOI: 10.1103/PhysRevE.87.052804 PACS number(s): 89.75.Hc, 64.60.ah, 89.75.FbI. INTRODUCTIONMost studies of complex networks over the last decadehave focused on single networks that do not interact withother networks [121]. Most real-world networks, however,are actually networks of networks and are made up of two ormore networks that interact and are interdependent. Examplesinclude various infrastructures such as water and food supplynetworks, communication networks, fuel networks, financialtransaction networks, and power-station networks [2226].Recently, Buldyrev et al. [27] developed a theoretical frame-work for studying the robustness of two fully interdependentnetworks subject to random attack. They found that due tothe dependency coupling between networks, they becomeextremely vulnerable to random failures and system collapse inan abrupt first order transition. Parshani et al. [28] generalizedthis framework and investigated networks that are only par-tially interdependent. They found that reducing the couplingstrength could lead to a change from first order (abrupt) to asecond order (continuous) percolation transition due to randomattack. Recently, Gao et al. [29] generalized this theory toinclude the network of networks (NON) case, a system of ninteracting networks. They developed an analytical frameworkfor studying percolation under random failures and presentedexact analytical solutions for an interdependent NON system.*dfocus.gao@gmail.comThe results for a NON suggest that the classical percolationtheory extensively studied in physics and mathematics is alimiting case of n = 1 of the general case of percolation in aNON. For more recent studies of a NON, see Refs. [16,30].In real-world scenarios, however, failures may not be randombut may be the result of a targeted attack on specific centralnodes [3135]. Huang et al. [31] studied the robustness oftwo interdependent networks when high- or low-degree nodesare under targeted attack by mapping the targeted-attackproblem in interdependent networks to the random-attackproblem. They found that interdependent scale-free networksare difficult to defend using strategies such as protecting thehigh-degree nodes. In this paper, by introducing a probabilityfunction of targeted attack, we study the robustness of a NONsystem of any tree of n fully interdependent networks and of astarlike NON structure of n partially interdependent networks(see Fig. 1).II. THE MODELGallos et al. [35] proposed a probability function fora targeted attack on an isolated single network, W(ki) =ki /Ni=1 ki , which represents the probability that a nodewith degree ki will be removed. For the case < 0, it wasassumed that ki = 0. This function was used to study therobustness of a single scale-free network to different attackstrategies [35], but in a system comprised of interdependent052804-11539-3755/2013/87(5)/052804(11) 2013 American Physical Societyhttp://dx.doi.org/10.1103/PhysRevLett.107.195701http://dx.doi.org/10.1103/PhysRevE.87.052804DONG, GAO, DU, TIAN, STANLEY, AND HAVLIN PHYSICAL REVIEW E 87, 052804 (2013)n-1(a)q12q21q41qq1ii1q13q14q3123 4i(b)1qn1 q1nnFIG. 1. Each circle represents one network. Two networks iand j form a partially interdependent pair of networks if a certainfraction qji > 0 of nodes of network i directly depends on nodes ofnetwork j ; i.e., they cannot function if the nodes in network j onwhich they depend do not function. Two networks i and j are fullyinterdependent if qji = qji = 1, thereby establishing a one-to-onecorrespondence. (a) An example of a treelike NON composed offive fully interdependent networks. (b) A starlike NON composedof n partially dependent networks. The arrows pointing on q1i andqi1(i = 1,2,3, . . . ,n) represent the partial dependency of the centralnetwork, 1, on a neighboring network i, and q1i represents thedependency of the neighboring network i on the central network.networks, zero-degree nodes strongly influence overall systemrobustness. In order to understand the robustness of a NONunder targeted attack and to solve the problem analytically wemust avoid singularities for ki = 0 in W(ki). We thus proposea probability function for targeted attack that represents theprobability that a node i with degree ki has been removed,W(ki) = (ki + 1)Ni=1(ki + 1), < < +. (1)When < 0, nodes with higher ki are protected, and nodeswith lower ki have a higher probability of failure. In orderto avoid singularities for ki = 0 and < 0, the term ki + 1is taken. When > 0, nodes with higher ki have a higherprobability of failure and are therefore vulnerable. The case corresponds to the targeted attack case where nodesare removed in strict order, from high degree to low degree.When = 0, all nodes in the networks have the sameprobability of failure, W0 = 1N , and the nodes are removed atrandom. The important special case equal to= 1 correspondsto the acquaintance immunization strategy [13,31,32,35] andcan be solved analytically.We assume that in an attack, a fraction of nodes 1 piaccording to Eq. (1) in each network i is removed and onlynodes which belong to the giant connected component of eachnetwork i remain functional. The failure of these nodes causesfurther failures in other networks and vice versa [2729,36,37],resulting in a cascade of failures. We now describe thegenerating function and the giant component of one networki in the NON after a fraction of nodes is removed with aprobability W(ki), which can be applied to all networksin the NON system. The generating function of the degreedistribution of network i is defined asGi0(x) kPi(k)xk, (2)where Pi(k) is the degree distribution of network i. Thegenerating function of the associated branching process isGi1(x) = Gi0(x)Gi0(1)[2729,31,32,3841]. The average degree ofnetwork i is ki =k Pi(k)k. We follow an approach similarto that of attacking a pair of coupled networks [31].Applying Eq. (1), after removing a fraction of 1 pi ofnodes of network i but keeping the edges of the remainingnodes, the degree distribution Ppi (k) of the remaining nodes isPpi (k) =Api (k)piN, (3)where Api (k) is the number of nodes with degree k in networki. When one more node is removed, Api (k) changes toA(pi1/N)(k) = Api (k) Ppi (k)(k + 1)k Ppi (k)(k + 1). (4)In the limit N , the derivative of Api (k) with respect topi is obtained from Eq. (4),dApidpi= N Ppi (k)(k + 1)k Ppi (k)(k + 1). (5)By differentiating Eq. (3) with respect pi and using Eq. (5),we obtain pi dPpi (k)dpi= Ppi (k) Ppi (k)(k + 1)k Ppi (k)(k + 1). (6)In order to solve Eq. (6), we define G(x) k Pi(k)x(k+1)and introduce a new variable, ui G1 (pi) [31,32,42]. Wefind the solution of Eq. (6) to bePpi (k) =1piPi(k)u(k+1)i , (7)kPpi (k)(k + 1) =uiG(ui)G(ui), (8)which can be shown to satisfy Eq. (6). After removing afraction 1 pi of nodes from network i according to Eq. (1)and their links, the generating function of the nodes left innetwork i isGib(x) kPpi (k)xk = 1pikPi(k)u(k+1)i xk. (9)In this step, after removing the links between removed nodesand remaining nodes, the generating function of network icomposed of the remaining nodes is [41,43]Gic(x) Gib(1 pi + pix), (10)where pi is the ratio between the number of links ofthe remaining nodes and the number of original links innetwork i,pi = piNik(pi)Nik=k Pi(k)ku(k+1)ik Pi(k)k, (11)where k(pi) =k Ppi (k)k is the average degree of theremaining nodes in network i.If we find the generating function Gi0(x) of network iafter randomly removing the 1 pi fraction of nodes, thenew generating function of the remaining nodes is the same asGic(x). The targeted-attack problem on network i is mapped to052804-2ROBUSTNESS OF NETWORK OF NETWORKS UNDER . . . PHYSICAL REVIEW E 87, 052804 (2013)a random-attack problem on network i by solving the equationGi0(1 pi + pix) = Gic(x) [31],Gi0(x) = Gib(1 pipi+ pipix). (12)The generating function of the associated branching processis [41,43]Gi1(x) = Gi0(x)Gi0(1). (13)When a fraction of 1 pi of nodes in the ith network isremoved with the probability given by Eq. (1), the fraction ofnodes that belongs to the giant component is given bygi(pi) = 1 Gi0[1 pi(1 fi)], (14)where fi fi(pi) satisfies a transcendental equationfi = Gi1[1 pi(1 fi)]. (15)In particular, for the case = 1, from Eqs. (9), (10),and (12), Gib(x),Gic(x), and Gi0(x) becomeGib(x) = uipikPi(k)uki xk = uipiGi0(uix), (16)Gic(x) = uipiGi0[ui(1 pi + pix)], (17)Gi0(x) = uipiGi0[ui(1 pipi+ pipix)], (18)where ui satisfiesui G1i0 (pi), (19)pi u2i Gi0(ui)Gi0(1). (20)Next, we apply the framework of the n equations developedin Ref. [36] to the no-feedback condition. The n unknowns xt,irepresent the fraction of nodes that survived in network i atstep t of the cascading failures after removing the failed nodesand the nodes dependent upon the failed nodes in the othernetworks,xt,i = piKj=1[1 qji + qjiyt,j igj (xt,j )], (21)where, for t = 1, xt,i = pi andyt,j i = xt,j1 qij + qij yt,ij gi(xt,i) . (22)The terms yt,j i represent the fraction of nodes in network jthat survive after the damage from all the networks connectedto j is taken into account (except from network i at stage t ofthe cascading failure) and yt=1,j i = pj . The giant componentof each network i at stage t of cascading failures can also befound from the equation Pt,i = xt,igi(xt,i), where gi(xt,i) isthe fraction of the remaining nodes of network i at stage t ofcascading failures that are part of its giant component.When t , the cascading failures have ended, and thefinal giant component of network i is P,i = xigi(xi), wherexi = piKj=1[1 qji + qjiyjigj (xj )], (23)yji = xj1 qij + qij yij gi(xi) , (24)where xi is calculated over the K networks interlinked withnetwork i by partial dependency links qji . The terms yjirepresent the fraction of nodes in network j surviving afterthe damage from all the connected networks to network jexcept from network i is taken into account [29,36,39]. Whenthe dependency has a feedback between two interdependentnetworks, Eq. (23) becomes simpler as yji = xj [39]. We nextaddress two cases of n coupled networks: (i) a treelike NONwith fully interdependent networks and (ii) a starlike NONwith partially interdependent networks.III. TREELIKE FULLY INTERDEPENDENTNETWORK OF NETWORKSWe now analyze the specific classes of a treelike NONformed of Erdos-Renyi (ER) or scale-free (SF) networks[2,3,10,11,27,31]. By using the dynamics of cascading failuresshown in Fig. 2, we first compare our theoretical results ofEqs. (16)(20) with simulation results shown in Fig. 3 for aNON formed of five ER networks and SF networks. The theoryis in excellent agreement with simulations.t=3,5,7,... t=4,6,8,...t=2t=12 34 52 314 52 314 5 4 52 311FIG. 2. The dynamic of cascading failures of a treelike NONsystem composed of five fully interdependent networks. Each noderepresents an ER network or SF network. The arrow (on the link)illustrates the direction of the damage spreading from one network toanother within the treelike NON at step t . At t = 1, arrows point tofive networks, which means that a fraction 1 pi of nodes in the fivenetworks is removed according to Eq. (1). Furthermore, we removenonfunctional nodes, which become disconnected from the largestcomponent of each network. At t = 2, arrows from network 1 pointto networks 2 and 3, illustrating that all nodes that depend on theremoved nodes in network 1 are removed and nonfunctional nodesare also removed in networks 2 and 3. In this step (t = 2), network1 spreads the damage to networks 2 and 3. At t = 3,5,7, . . . , arrowsfrom networks 2 and 3 point to networks 1, 4, and 5. This reflects thatnetworks 2 and 3 spread damage to their neighboring networks, 1, 4,and 5. At t = 4,6,8, . . . , networks 1, 4, and 5 spread damage to theirneighboring networks, 2 and 3. Thus, the damage is spread betweennetworks, back and forth, until they arrive at a mutually stable giantcomponent or collapse.052804-3DONG, GAO, DU, TIAN, STANLEY, AND HAVLIN PHYSICAL REVIEW E 87, 052804 (2013)10 20 30 40 50 6000.20.40.60.81tPt,1theorysimulation(a) = 1k = 4n = 5p = 0.9900ER10 20 30 40 500.750.80.850.90.95Pt,1ttheorysimulation(b) = 1k = 4n = 5p = 0.994ER10 20 30 40 5000.20.40.60.81tPt,1theorysimulation(c) = 1 = 2.7p = 0.990SF10 20 30 40 500.920.940.960.98tPt,1theorysimulation(d) = 1 = 2.7p = 0.9920SFFIG. 3. For the case of = 1, which corresponds to the solvable acquaintance immunization strategy, we present the giant component(theory and simulation) of network 1, Pt,1, during the t cascading failures for fully interdependent treelike NON composed of five (a) and (b) ERand (c) and (d) SF networks. (a) For each ER network in the NON, we chose ki = k = 4, Ni = N = 104, and pi = p = 0.990 < pc = 0.9920.(b) The giant component of network 1 (Fig. 2) with the same parameters as in (a) but for pi = p = 0.9940 > pc = 0.9920. (c) For each SFnetwork in the NON, we chose = 1, = 2.7, Ni = N = 104, and pi = p = 0.9900 < pc = 0.9910. (d) The giant component with the sameparameters as in (c) but for pi = p = 0.9920 > pc = 0.9910. The results of the simulations are obtained by averaging over 50 realizations.For the case of a treelike NON formed of n ER networks,the generating function of the ER network isGi0(x) = eki (x1).From Eqs. (13) and (18)(20), we obtain Gi0(x) = Gi1(x) =ekiu2i (x1) for = 1, where ui can be expressed in terms of theLambert function ui = W (pi ki eki )ki . From Eqs. (14) and (15),we can get gi(xi) = 1 Gi0[1 xi(1 fi)], fi = Gi1[1 xi(1 fi)], and gi(xi) = 1 fi . For a treelike NON with fullyinterdependent qij = qji = 1, we getxi = pinj=1yji(1 fj ), (25)yji = xjyij (1 fi) , (26)P = P,i =ni=1pigi(xi) =ni=1pi(1 fi), (27)xi = P,igi(xi)=nj=1 pj (1 fj )1 fi , (28)where fi = Gi1[1 xi(1 fi)] = ekiu2inj=1 pj (1fj ),i = 1,2, . . . ,n.When pi = p and ki = k, we have ui = u, and fi = fsatisfiesf = eku2pn(1f )n . (29)Then, from Eq. (27), we obtainP =ni=1pigi(xi) = pn(1 f )n. (30)From Eqs. (29) and (30), we can getf = eku2P . (31)The size of the mutual giant component for all values of p, k,u and n becomesP = pn(1 eku2P )n, (32)whereu = W (pkek)k. (33)Let r = epnu2(1f )n ; then r = f 1k ,f = rk , and Eq. (29)becomes r = epnu2(1rk )n . Then the critical case correspondsto the tangential condition, 1 = ddrepnu2(1rk )n , and the criticalvalue of p = pc can be found,pnc =1nf (1 f )n1ku2 . (34)052804-4ROBUSTNESS OF NETWORK OF NETWORKS UNDER . . . PHYSICAL REVIEW E 87, 052804 (2013)0.94 0.96 0.98 100.20.40.60.81pP,itheorysimulation = 1n = 5k = 4.5 k = 4(a)k = 50.94 0.96 0.98 100.20.40.60.81pP,1theorysimulationk = 5 k = 4 = 1n = 5(b)k = 4.50.94 0.96 0.98 10100200300400p(c) = 1n = 5 k = 5k = 4.5k = 40.94 0.96 0.98 1050100150200250300p(d) = 1n = 5 k = 4.5k = 4k = 5FIG. 4. (a) and (b) Comparison between theory and simulations of the giant component of network i, P,i , of a treelike NON systemformed of five ER networks as a function of p. (c) and (d) Simulations of the number of cascading stages as a function of p for different k.In the simulations, Ni = N = 104, and the results are obtained by averaging over 50 realizations. The sharp peaks represent the location of pc.Substituting Eq. (34) into Eq. (29), we can get the critical valueof f = fc,fc = efc1nfc , (35)where the solution can be expressed in terms of the Lambertfunction W (x), fc = [nW ( 1ne1n )]1.We support our analytical results, Eqs. (32), (33), (9)(15),(23), and (24), by simulating the giant component of thenetworks as a function of p, as shown in Figs. 4(a) and 4(b).The results of the theory agree well with the simulations. Theaverage convergence stage as a function of p for several kis simulated and shown in Figs. 4(c) and 4(d). From Figs. 4(c)and 4(d) the first order transition threshold can be easilyidentified by the sharp peak, which represents the divergenceof the number of cascading stages [44].As an important characteristic of percolation theory, thecritical threshold pc is analyzed to study the robustness ofa fully interdependent treelike NON in Figs. 5(a)5(f). Thethreshold pc as a function of obtained from Eqs. (9)(15),(23), and (24) is shown in Fig. 5(a). Figure 5(a) shows that thecritical fractionpc continuously increases with for a differentnumber n of networks, which means that the NON becomesmore vulnerable when the nodes with higher k have a higherprobability of failure. Figures 5(b) and 5(c) show that the NONbecomes less robust with increasing n. Moreover, for fixedn, pc decreases when k increases, as seen in Figs. 5(b)5(f).Furthermore, for a fixed n, when k is smaller than the minimumaverage degree kmin(n), pc = 1. In this case [k kmin(n)] theNON will collapse even if a single node fails. Using Eq. (34),for pc = 1, we findkmin(n) = 1nfc(1 fc)n1 , (36)which is independent of (see Ref. [36] for = 0). In fact,when all n interdependent ER networks have a minimumaverage degree kmin, the ER networks are at critical condition,and the NON collapses even if a single node fails irrespectiveof whether it is high degree or low degree. This is an inherentproperty of interdependent networks and is not influenced bythe probability that nodes are attacked.For the case of a treelike NON formed of n SF networks,Pi(k) = (k+1)1k1(M+1)1m1 is the degree distribution of nodes ofnetwork i, where k,M , andm are the degree, maximum degree,and minimum degree of network i, respectively [4,27,31]. Thecritical thresholds pc of the NON are obtained from Eqs. (9)(15), (23), and (24) and are presented in Fig. 6. Figure 6(a)shows that pc for interdependent SF networks is nonzero forthe entire range of , which differs from that of a single SFnetwork for whichpc = 0 [4,5,31,35]. For fixedn,pc increaseswhen increases since the NON becomes less robust andthe high-degree nodes are removed with a higher probability.Also, with the same probability of removing nodes (fixed )the NON system becomes more vulnerable when the numberof networks increases, as seen in Fig. 6. When lower-degreenodes are removed with a higher probability, Fig. 6(b) showsthat the NON becomes less robust with increasing . This isalso illustrated for = 0 and 1 in Figs. 6(c) and 6(d), wherepc becomes larger with larger .052804-5DONG, GAO, DU, TIAN, STANLEY, AND HAVLIN PHYSICAL REVIEW E 87, 052804 (2013)2 0 20.70.750.80.850.90.951pcn=2n=3n=4(a)k=41010.60.70.80.91npc = 1k = 4k = 6k = 5k = 7(b)1010.60.70.80.91npc = 0k = 4k = 5k = 7k = 6(c)3 4 5 60.650.750.850.9510.55kpcn=2n=3n=4n=7(d) = 13 4 5 60.650.750.850.9510.55kpcn=2n=3n=4n=7 = 0(e)3 4 5 60.550.650.750.850.951kpcn=2n=3n=4n=7(f) = 1FIG. 5. (a) The critical threshold pc as a function of for different n. The critical threshold pc as a function of n for (b) = 1 and (c) = 0 for several values of k. (d)(f) For = 1,0, 1, pc as a function of k for different n. Note that for k kmin, pc = 1; i.e., the networkcollapses even if a single node is removed.IV. THE PARTIALLY INTERDEPENDENT STARLIKE NONWe now study the robustness of a starlike partially in-terdependent NON (see Fig. 7) using the analytical frame-work provided in Eqs. (23) and (24). Using the dynamicsof cascading failures shown in Fig. 7, we compare thetheoretical results of Eqs. (9)(15), (21), and (22) with thesimulation results in Fig. 8 for a starlike NON formed of fiveER networks.For ki = k, qi1 = q1i = q, p1 = p2 = = pn = p, and = 1, by using the generating function of ER networks andEqs. (16)(20), (23), and (24), we obtainx1 = p[1 q + pqg2(x2)]n1 = p[1 q + pq(1 f2)]n1,(37)x2 = p{1 q[1 y12g1(x1)]}, (38)wherey12 = x11 q[1 y21g2(x2)] . (39)From Eqs. (37)(39), we obtainx2 =p{1 q +pq(1 f1)[1 q +pq(1 f2)]n2}, (40)wheref1 = eku2(f11)x1 , (41)f2 = eku2(f21)x2 , (42)ui = u = W (pkek)k. (43)052804-6ROBUSTNESS OF NETWORK OF NETWORKS UNDER . . . PHYSICAL REVIEW E 87, 052804 (2013)2 1 0 1 2 30.650.70.750.80.850.90.95pc = 2.7(a)n = 2n = 31010.50.60.70.80.91npc = 1 = 3 = 2.7 = 2.3(b)2.5 3 3.50.50.60.70.80.91pcn = 2n = 4n = 3 = 0(c)2.5 3 3.50.90.920.940.960.981pc = 1n = 4n = 3n = 2(d)FIG. 6. (a) The critical threshold pc as a function of for a NON composed of n = 2 and 3 SF networks with = 2.7. (b) The criticalthreshold pc as a function of n for = 1 and different . (c) and (d) pc as a function of for = 0,1 and for different n, respectively. Thelowest degree of the SF networks is taken to be m = 2.Substituting Eqs. (37) and (40) into Eqs. (41) and (42), weobtainf1 = ekpu2(f11)[1q+pq(1f2)]n1 , (44)f2 = ekpu2(f21){1q+pq(1f1)[1q+pq(1f2)]n2}. (45)t=15 13 33 14545214t=2,4,6,8,10,... t=3,5,7,9,11,...22FIG. 7. The dynamics of cascading failures for a starlike NONwith five networks. Each node represents a network, and the arrow (onthe link) illustrates the damage spreading from one network to anothernetwork in the NON. At t = 1, we remove a fraction 1 pi of nodesin the five networks according to Eq. (1). Furthermore, unconnectednodes in the networks are also removed. At t = 2,4,6, . . . , the centralnetwork, 1, spreads damage to the neighboring networks, 2, 3, 4, and5. At t = 3,5,7, . . . , networks 2, 3, 4, and 5 spread damage to centralnetwork 1. Finally, the mutual giant component of the NON becomesstable, and no further removal of nodes occurs.Furthermore, the giant components of networks areP,1 = x1g1(x1) = p(1 f1)[1 q + pq(1 f2)]n1= ln f1ku2, (46)P,2 = x2g2(x2) = p(1 f2){1 q + pq(1 f1) [1 q + pq(1 f2)]n2} = ln f2ku2, (47)where f1 and f2 are obtained from Eqs. (44) and (45).Figures 9(a) and 9(b) show excellent agreement betweenthe simulations of the giant components and the results inEqs. (44) and (46). Furthermore, the central network, 1, showsa first and a second order percolation transition by regulatingthe coupling strength q, as shown in Figs. 9(a) and 9(b).Figures 9(c) and 9(d) present simulation results of as afunction of p. The first order transition point pIc can be easilyidentified by the sharp peaks of .From Eqs. (44) and (45), we obtainf1 = 1 ln f2kpu2(f21) + q 1pq[1 q + pq(1 f2)]n2 , (48)f2 = 1 1pq[(ln f1kpu2(f1 1)) 1n1+ q 1], (49)where u = W (pkek )k.From Eqs. (46), (48), and (49), the critical threshold pcas a function of n for different q is obtained and shown inFig. 10(a). We see that for = 1 the central network becomes052804-7DONG, GAO, DU, TIAN, STANLEY, AND HAVLIN PHYSICAL REVIEW E 87, 052804 (2013)20 40 6000.20.40.60.8tPt,1theorysimulation(a) = 1,k = 4,n = 5q = 0.8,p = 0.958510 20 30 4000.10.20.30.40.50.6tPt,itheorysimulation(b)i = 1i = 2, 3, 4, 5 = 1,k = 4,n = 5q = 0.4,p = 0.7430FIG. 8. For the case of = 1, which corresponds to the solvable acquaintance immunization strategy, we present the theory, andsimulated results of the giant component of networks Pt,i , i = 1,2,3,4,5 during the t cascading failure for a partially interdependent starlikeNON composed of five ER networks are shown. (a) For each network in the NON, ki = k = 4, Ni = N = 104, q1i = qi1 = 0.8, and pi =p = 0.9585 < pc = 0.9600. (b) Theory and simulations of the giant component Pt,i of networks as function of t with the same parameters asin (a) but forpi = p = 0.7430 > pc = 0.7420 and q1i = qi1 = 0.4. The results of the simulations are obtained by averaging over 50 realizations.significantly more vulnerable with increasing n for the samecoupling strength q. In addition, for the same n the centralnetwork becomes more robust when the coupling strengthq is reduced. Figure 10(b) shows that using Eqs. (9)(15),(23), and (24) these conclusions can also be demonstratedfor = 1. Figures 10(a) and 10(b) also show that thecentral network becomes more robust for the same couplingstrength q and n when nodes of lower degree have a higherprobability of failure. We also analyze pc as a function of theaverage network degree for different , as shown in Figs. 10(c)and 10(d). Figure 10(c) shows that the central networkbecomes more fragile and finally collapses as the averagedegree is decreased for the same coupling strength. ComparingFigs. 10(c) and 10(d), we see that when nodes of lowerdegree have a higher failure probability ( = 1), the centralnetwork is more robust for the same average degree and fixedcoupling strength. Therefore, in the partially starlike NON, ifwe decrease the number of networks and the coupling strength0.7 0.8 0.9 100.20.40.60.8pP,1theorysimulation = 1k = 4n = 5q = 0.6q = 0.4(a)q = 0.80.4 0.6 0.8 100.20.40.60.81P,1ptheorysimulation(b) = 2k = 4n = 5q = 0.6 q = 0.8q = 0.40.7 0.8 0.9 1050100150200250300p = 1k = 4n = 5(c)q = 0.6q = 0.80.4 0.6 0.8 1050100150200250300350p(d) = 2k = 4n = 5q = 0.6q = 0.8FIG. 9. (a) Comparison between simulation and theory for the giant component P,1 as a function of p for different q for the starlike NON.For weak coupling strength q = 0.4 and = 1, the central network shows a second order phase transition, and for strong coupling strengthq = 0.6,0.8, the central network shows a first order phase transition. (b) Comparison between simulation and theory for the giant componentP,1 with the same parameters in (a) but for = 2. (c) and (d) The average time at the convergence stage as a function of p for different and q. The results are obtained by averaging over 50 realizations for N = 104.052804-8ROBUSTNESS OF NETWORK OF NETWORKS UNDER . . . PHYSICAL REVIEW E 87, 052804 (2013)1010.50.60.70.80.91npcq = 0.3q = 0.4q = 0.5q = 0.7 = 1k = 4(a)1010.20.40.60.81pcn(b)q = 0.4q = 0.3q = 0.5q = 0.7 = 1k = 42 3 4 5 60.60.810.30.4kpc = 1n = 5q = 0.5q = 0.3q = 0.7q = 0.9(c)2 3 4 5 60.40.60.810.3pckq = 0.5q = 0.7q = 0.9 = 1n = 5q = 0.3(d)FIG. 10. Critical threshold pc as a function of n for several q values and average degree k. (a) and (b) pc as a function of n for = 1, 1and for different coupling strengths q. (c) and (d) pc as a function of k for different coupling strengths for = 1, 1.and increase the average network degree, it will become morerobust.From Eqs. (48) and (49), the second order transitionthreshold pIIc can be obtained by taking the limit f1 1,ln f2f2 1 = kpIIc u2(1 q), (50)f2 = 1 1pIIc q[(1kpIIc u2) 1n1+ q 1]. (51)From Eq. (50), we getf2 = 1kpIIc u2(1 q)W(kpIIc u2(1 q)e(kpIIc u2(1q))), (52)and by substituting Eq. (52) into Eq. (51), we obtain pIIc .Figures 9(a) and 9(b) show that changing the couplingstrength can change the phase transition from first order tosecond order. Thus we examine the relation betweenpc and thefraction q of interdependent nodes when, using Eq. (1), 1 pnodes are removed from the central network. Figure 11 showsthat for fixed average degree k and , there exists a criticalthreshold qc of coupling strength that changes for differentn. When the coupling is strong (q > qc), the NON exhibitsa first order phase transition at a transition threshold pIc atwhich the giant component abruptly approaches zero. Whenthe coupling is weak (q < qc), the NON shows a second orderphase transition at transition threshold pIIc at which the giantcomponent smoothly changes from a finite value to zero. Notethat the phase transition lines of both pIc and pIIc converge atq = 0 for different n, which is the limit of a single network,i.e., the classical percolation of a single network under targetedattack [35]. Figure 11 also shows that when q = 0, the centralnetwork becomes more vulnerable when the higher-k nodeshave a higher failure probability.In summary, by introducing a targeted attack probabilityfunction which avoids singularities for degree ki = 0 and < 0, we address the question of robustness of two types ofNON under targeted attack: (i) a treelike network comprised ofn fully interdependent Erdos-Renyi or scale-free networks and(ii) a starlike network comprised of n partially interdependentErdos-Renyi networks. For the treelike NON of n fullyinterdependent ER networks under targeted attack, we findan exact solution for the giant component P when = 1as a function of the 1 p initial fraction of removed nodesfrom each of the n networks. Analytical solutions are foundfor the critical fraction pc that leads to the fragmentation ofn interdependent networks, and the minimum average kminat which the NON becomes so vulnerable that the failureof a single node will bring down the entire system. Fordifferent values of , we present the numerical solutionfor the giant component and the critical fraction pc as afunction of k and n for ER and SF networks. When a treelikenetwork of n fully interdependent ER networks is undertargeted attack, increasing the average degree of each of thenetworks, increasing the probability of targeted attacks onlower degree nodes, and decreasing the number of networkswill make the tree more robust. For a starlike network of npartially interdependent ER networks we can use these samestrategies, but we can also decrease the coupling strength052804-9DONG, GAO, DU, TIAN, STANLEY, AND HAVLIN PHYSICAL REVIEW E 87, 052804 (2013)0 0.2 0.4 0.6 0.800.20.40.60.811q1 pcSecond orderCritical lineFirst order(a)n = 2n = 3n = 5n = 11n = 7k = 5 = 10 0.2 0.4 0.6 0.800.20.40.60.811 pc1qSecond orderCritical lineFirst order(b)n = 2n = 3n = 5n = 7n = 11k = 5 = 00 0.2 0.4 0.6 0.800.20.40.60.811 pc1qSecond orderCritical lineFirst order(c)n = 2n = 3n = 5n = 11k = 5 = 1n = 7FIG. 11. The relation between 1 q and 1 pc for different n for (a) = 1, (b) = 0, and (c) = 1. The curves connecting the circlesshow the critical lines of qc, below which the system shows a first order phase transition and above which a second order phase transition.to make the starlike NON more robust. We also find thecritical threshold qc for different values of below whichthe system exhibits a second order transition and above whichthe system shows a first order transition. When the couplingstrength q = 0, the NON system becomes n single andindependent networks. When is increasing or decreasing, thenetwork becomes more vulnerable or robust, which coincideswith the classic percolation of a single network to targetednetwork.ACKNOWLEDGMENTSWe thank X. Huang for discussions. L.T. acknowledgesthe support from the National Natural Science Foundation ofChina (Grants No. 11171135, No. 71073071, No. 71073072,and No. 51276081), Major Program of the National SocialScience Foundation of China (Grant No. 12&ZD062). G.D.acknowledges the support from the China Scholarship Fund(Grant No. 2011832326). R.D. acknowledges the supportfrom the Graduate innovative Foundation of Jiangsu ProvinceCX10B_272Z and the Youth Foundation of ChongqingNormal University (Grant No. 10XLQ001). J.G. thanks theDoctoral visiting scholar program of SJTU, the ShanghaiKey Basic Research Project (Grant No. 09JC1408000) andthe National Natural Science Foundation of China (Grant No.61004088) for support. S.H. thanks the European EPIWORKand MULTIPLEX (EU-FET project 317532) projects, theDeutsche Forschungsgemeinschaft (DFG), the Israel ScienceFoundation, and the LINC ITN Project for financial support.We wish to thank ONR (Grant N00014-09-1-0380, GrantN00014-12-1-0548), DTRA (Grant HDTRA-1-10-1-0014,Grant HDTRA-1-09-1-0035), NSF (Grant CMMI 1125290),CONGAS (Grant FP7-ICT-2011-8-317672), the Next Gener-ation Infrastructure (Bsik).[1] D. J. Watts and S. H. Strogatz, Nature (London) 393, 440(1998).[2] A.-L. Barabasi and R. Albert, Science 286, 509 (1999).[3] R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, 47(2002).[4] R. Cohen, K. Erez, D. ben-Avraham, and S. Havlin, Phys. Rev.Lett. 85, 4626 (2000); 86, 3682 (2001).[5] D. S. Callaway, M. E. J. Newman, S. H. Strogatz, and D. J.Watts, Phys. Rev. Lett. 85, 5468 (2000).[6] S. N. Dorogovtsev and J. F. F. Mendes, Evolution of Networks:From Biological Nets to the Internet and WWW (OxfordUniversity Press, New York, 2003).[7] R. P. Satorras and A. Vespignani, Evolution and Structure of theInternet: A Statistical Physics Approach (Cambridge UniversityPress, Cambridge, 2006).[8] A. Bashan, R. P. Bartsch, J. W. Kantelhardt, S.Havlin, and P. C. Ivanov, Nat. Commun. 3, 702(2012).052804-10http://dx.doi.org/10.1038/30918http://dx.doi.org/10.1038/30918http://dx.doi.org/10.1126/science.286.5439.509http://dx.doi.org/10.1103/RevModPhys.74.47http://dx.doi.org/10.1103/RevModPhys.74.47http://dx.doi.org/10.1103/PhysRevLett.85.4626http://dx.doi.org/10.1103/PhysRevLett.85.4626http://dx.doi.org/10.1103/PhysRevLett.86.3682http://dx.doi.org/10.1103/PhysRevLett.85.5468http://dx.doi.org/10.1038/ncomms1705http://dx.doi.org/10.1038/ncomms1705ROBUSTNESS OF NETWORK OF NETWORKS UNDER . . . PHYSICAL REVIEW E 87, 052804 (2013)[9] C. Song, S. Havlin, and H. A. Makse, Nature (London) 433, 392(2005); Nat. Phys. 2, 275 (2006).[10] R. Cohen and S. Havlin, Complex Networks: Structure, Robust-ness and Function (Cambridge University Press, Cambridge,2010).[11] G. Caldarelli and A. Vespignani, Large-Scale Structure andDynamics of Complex Webs (World Scientific, Singapore, 2007).[12] M. E. J. Newman, Networks: An Introduction (Oxford UniversityPress, New York, 2010).[13] R. Cohen, S. Havlin, and D. ben-Avraham, Phys. Rev. Lett. 91,247901 (2003); D. Zhou, H. E. Stanley, G. DAgostino, andA. Scala, Phys. Rev. E 86, 066103 (2012).[14] S. Son, G. Bizhani, C. Christensen, P. Grassberger, andM. Paczuski, Europhys. Lett. 84, 16006 (2012).[15] Z. Wang, A. Szolnoki, and M. Perc, Europhys. Lett. 97, 48001(2012).[16] G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, and J. F. F.Mendes, Phys. Rev. Lett. 109, 248701 (2012).[17] Z. Wang, A. Szolnoki, and M. Perc, Sci. Rep. 3, 1183 (2013).[18] S. Amina, G. A. Schwartzb, and S. S. Sastryb, Automatica 49,1 (2013).[19] J. F. Castet and J. H. Saleh, PLoS One 8, e60402 (2013).[20] C. D. Brummitta, R. M. DSouzab, and E. A. Leichtf, Proc. Natl.Acad. Sci. USA 109, E680 (2012).[21] C. M. Schneider, A. A. Moreira, J. S. Andrade Jr., S. Havlin,and H. J. Herrmann, Proc. Natl. Acad. Sci. USA 108, 3838(2011); C. M. Schneider, N. Yazdani, N. A. M. Araujo, S. Havlin,and H. J. Herrmann, arXiv:1106.3234.[22] J. C. Laprie, K. Kanoun, and M. Kaniche, in Computer Safety,Reliability, and Security: 26th International Conference, SAFE-COMP 2007, Nuremberg, Germany, September 1821, 2007.Proceedings, Lecture Notes in Computer Science, Vol. 4680(2007), p. 54.[23] J. Goldenberg, Y. Shavitt, E. Shir, and S. Solomon, Nat. Phys.184, 1 (2005).[24] P. Pederson, D. Dudenhoeffer, S. Hartley, and M. Permann,Technical Report INL/EXT-06-11464, Idaho National Labora-tory, 2006.[25] V. Rosato, Int. J. Crit. Infrastruct. 4, 63 (2008).[26] S. Rinaldi, J. Peerenboom, and T. Kelly, IEEE Contr. Syst. Mag.21, 11 (2001).[27] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin,Nature (London) 464, 1025 (2010).[28] R. Parshani, S. V. Buldyrev, and S. Havlin, Phys. Rev. Lett. 105,048701 (2010).[29] J. Gao, S. V. Buldyrev, S. Havlin, and H. E. Stanley, Phys. Rev.Lett. 107, 195701 (2011).[30] K. Zhao and G. Bianconi, J. Stat. Mech. (2013) P05005.[31] X. Huang, J. Gao, S. V. Buldyrev, S. Havlin, and H. E. Stanley,Phys. Rev. E 83, 065101(R) (2011).[32] G. Dong, J. Gao, L. Tian, R. Du, and Y. He, Phys. Rev. E 85,016112 (2012).[33] I. Simonsen, L. Buzna, K. Peters, S. Bornholdt, and D. Helbing,Phys. Rev. Lett. 100, 218701 (2008).[34] A. E. Motter and Y.-C. Lai, Phys. Rev. E 66, 065102(2002).[35] L. K. Gallos, R. Cohen, P. Argyrakis, A. Bunde, and S. Havlin,Phys. Rev. Lett. 94, 188701 (2005).[36] J. Gao, S. V. Buldyrev, S. Havlin, and H. E. Stanley, Phys. Rev.E 85, 066134 (2012); T. Kalisky, R. Cohen, D. ben-Avraham,and S. Havlin, Europhys. Lett. 38, 269 (2004).[37] S. Havlin, N. A. M. Araujo, S. V. Buldyrev, C. S. Dias,R. Parshani, G. Paul, and H. E. Stanley, arXiv:1012.0206v1;D. Zhou, J. Gao, H. E. Stanley, and S. Havlin,arXiv:1206.2427v2.[38] Y. Hu, B. Ksherim, R. Cohen, and S. Havlin, Phys. Rev. E 84,066116 (2011).[39] J. Gao, S. V. Buldyrev, S. Havlin, and H. E. Stanley, Nat. Phys.8, 40 (2012).[40] Q. Li, L. A. Braunstein, S. Havlin, and H. E. Stanley, Phys. Rev.E 84, 066101 (2011).[41] M. E. J. Newman, Phys. Rev. E 66, 016128 (2002).[42] J. Shao, S. V. Buldyrev, S. Havlin, and H. E. Stanley, Phys. Rev.E 83, 036116 (2011).[43] J. Shao, S. V. Buldyrev, L. A. Braunstein, S. Havlin, and H. E.Stanley, Phys. Rev. E 80, 036105 (2009).[44] R. Parshani, S. V. Buldyrev, and S. Havlin, Proc. Natl. Acad.Sci. USA 108, 1007 (2011).052804-11http://dx.doi.org/10.1038/nature03248http://dx.doi.org/10.1038/nature03248http://dx.doi.org/10.1038/nphys266http://dx.doi.org/10.1103/PhysRevLett.91.247901http://dx.doi.org/10.1103/PhysRevLett.91.247901http://dx.doi.org/10.1103/PhysRevE.86.066103http://dx.doi.org/10.1209/0295-5075/97/16006http://dx.doi.org/10.1209/0295-5075/97/48001http://dx.doi.org/10.1209/0295-5075/97/48001http://dx.doi.org/10.1103/PhysRevLett.109.248701http://dx.doi.org/10.1016/j.automatica.2012.09.001http://dx.doi.org/10.1016/j.automatica.2012.09.001http://dx.doi.org/10.1371/journal.pone.0060402http://dx.doi.org/10.1073/pnas.1110586109http://dx.doi.org/10.1073/pnas.1110586109http://dx.doi.org/10.1073/pnas.1009440108http://dx.doi.org/10.1073/pnas.1009440108http://arXiv.org/abs/arXiv:1106.3234http://dx.doi.org/10.1504/IJCIS.2008.016092http://dx.doi.org/10.1109/37.969131http://dx.doi.org/10.1109/37.969131http://dx.doi.org/10.1038/nature08932http://dx.doi.org/10.1103/PhysRevLett.105.048701http://dx.doi.org/10.1103/PhysRevLett.105.048701http://dx.doi.org/10.1103/PhysRevLett.107.195701http://dx.doi.org/10.1103/PhysRevLett.107.195701http://dx.doi.org/10.1088/1742-5468/2013/05/P05005http://dx.doi.org/10.1103/PhysRevE.83.065101http://dx.doi.org/10.1103/PhysRevE.85.016112http://dx.doi.org/10.1103/PhysRevE.85.016112http://dx.doi.org/10.1103/PhysRevLett.100.218701http://dx.doi.org/10.1103/PhysRevE.66.065102http://dx.doi.org/10.1103/PhysRevE.66.065102http://dx.doi.org/10.1103/PhysRevLett.94.188701http://dx.doi.org/10.1103/PhysRevE.85.066134http://dx.doi.org/10.1103/PhysRevE.85.066134http://arXiv.org/abs/arXiv:1012.0206v1http://arXiv.org/abs/arXiv:1206.2427v2http://dx.doi.org/10.1103/PhysRevE.84.066116http://dx.doi.org/10.1103/PhysRevE.84.066116http://dx.doi.org/10.1038/nphys2180http://dx.doi.org/10.1038/nphys2180http://dx.doi.org/10.1103/PhysRevE.84.066101http://dx.doi.org/10.1103/PhysRevE.84.066101http://dx.doi.org/10.1103/PhysRevE.66.016128http://dx.doi.org/10.1103/PhysRevE.83.036116http://dx.doi.org/10.1103/PhysRevE.83.036116http://dx.doi.org/10.1103/PhysRevE.80.036105http://dx.doi.org/10.1073/pnas.1008404108http://dx.doi.org/10.1073/pnas.1008404108