Dynamical leaps due to microscopic changes in multiplex networks
DDynamical leaps due to microscopic changes in multiplex networks
Marina Diakonova,
1, 2
Jos´e J. Ramasco, and V´ıctor M. Egu´ıluz School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, United Kingdom Instituto de F´ısica Interdisciplinar y Sistemas Complejos IFISC (CSIC-UIB), Campus UIB, 07122 Palma de Mallorca, Spain (Dated: November 7, 2018)Recent developments of the multiplex paradigm included efforts to understand the role played bythe presence of several layers on the dynamics of processes running on these networks. The possibleexistence of new phenomena associated to the richer topology has been discussed and examples ofthese differences have been systematically searched. Here, we show that the interconnectivity ofthe layers may have an important impact on the speed of the dynamics run in the network andthat microscopic changes such as the addition of one single inter-layer link can notably affect thearrival at a global stationary state. As a practical verification, these results obtained with spectraltechniques are confirmed with a Kuramoto dynamics for which the synchronization consistentlydelays after the addition of single inter-layer links.
The ubiquity of processes naturally described by dy-namics on networks raises important issues on the roleplayed by network structure on the emergence of collec-tive phenomena. It has been shown that spectral analysisof the associated adjacency and Laplacian matrices canoffer insights into a variety of fundamental phenomenasuch as those relying upon spreading or diffusion mecha-nisms [1–3]. In particular, spectral methods are the basisto characterize synchronization and random walk diffu-sion in networks [4–7]. The second smallest eigenvalue λ is known to be related to the timescale to synchro-nization [8] and consensus [9], and is often interpretedas the ‘proper time’ of the system to relax [10]. Thisquantity, which depending on the literature is known as‘algebraic connectivity’ or ‘spectral gap’, is also indica-tive of the time of diffusion [11]. While the role of λ in“simple” graphs is well understood, more effort is neededto characterize its role and behavior in more realistic (andtherefore complex) contexts.One such novel framework is that of a multilayer net-work [12–14]. Its usefulness extends from finance [15–17]and mobility [18, 19], to epidemics [20, 21] and societaldynamics [22–28]. The multiplex scenario is ostensiblynon-trivial, in the sense that the phenomena observedon this system of interconnected networks cannot bestraightforwardly reduced to an aggregate network [29].The implication is that the multiplex structure plays afundamental role in diffusive processes, and that there-fore its effect of on λ is a pertinent and open issue.Diffusion on a multiplex, that is, a system with layersof N nodes with an association of corresponding nodesin each layer, were considered in [30, 31]. By varying theinter-layer diffusion constant, these works found bound-aries of λ for the multiplex in terms of the values forindividual layers (see also [32]). For a general weightedtwo-layer multiplex with a varying inter-layer link weight p , these results were rephrased in terms of a structuraltransition [33]: it was found that below some critical p c , λ ( p ) ∼ p , while for larger p the algebraic connectivityapproaches an asymptote given by L = λ ( L + L ), where L , are the Laplacian matrices of the two layers.The changing character of the functional growth of λ suggests existence of two regimes, an ‘underconnected’one with the two layers functioning autonomously, and a‘multiplex’ one where the inter-layer connectivity playsthe dominant role. It has been noted that the existenceof the point of inflection p c might follow from linear al-gebra arguments [34]. At this point it is also not cer-tain whether p c tends to zero in the limit of N → ∞ , orwhether it is indeed significant in describing the observedchanges in the dynamic timescales [35].This work extends the analysis of multiplex timescalesto cover the effect of the specific type of inter-layer cou-pling. Instead of gradually tuning up the intensity ofinter-layer links p , we now switch on the inter-layer linksat unit intensity one by one. Our strategy of consider-ing binary weight on the inter-layer links corresponds toa natural situation where these links are either absentor present, and where one can preferentially control theamount, and the placement, of such links. Spectral prop-erties resulting from such inter-layer link sequences, andtheir placement, were considered in [10]. The authorsstudied identical partially-coupled networks, and foundthat the trends of average algebraic connectivity showeda qualitative change at an (analytically derivable) mini-mum number of inter-layer links, depending on the place-ment of their end points. The authors of [36] elucidatedthe optimal relation between the similarity of layers andthe weight of inter-layer link, given some minimizationconstraint. Here we consider non-identical, general lay-ers (see [37] for study of a multiplex with widely differentinter-layer structure), and instead focus on the effect ofthe addition of inter-layer links to individual systems.We find that λ ( q ) of the fraction q of such an inter-layer connection sequence is characterized by leaps, im-plying sudden decreases in the timescales of the resultingmultiplex. We characterize the statistics of these jumpsand elucidate the way in which layer degree correlation[38, 39] affect them. Our results show that not only arehighly correlated multiplexes connected with high inten- a r X i v : . [ phy s i c s . s o c - ph ] A p r http://ifisc.uib-csic.es The Workflow ascending descending ordered
PC NC
Interlayer Correlation Sequence Type
Uncorrelated (
UnC ) Positively Correlated ( PC ) Negatively Correlated ( NC ) Random Random, Ascending, Descending Random, Ordered Figure 1. Some connection sequence schema for a positively-correlated and negatively-correlated multiplex. sities on average slower than anti-correlated multiplexes,but that the statistics of jumps are qualitatively differ-ent for distinct types of ordering. Finally, we validateour findings by running Kuramoto dynamics on the mul-tiplex, and verifying that λ does indeed inform on thescaling of the approach to the fully-synchronized state.The first step is to describe how the networks used inthe analysis are built. We consider multiplexes formedof statistically equivalent networks with the same num-ber of nodes N on each layer. For simplicity, the multi-plexes are formed of only two layers G and G , whichalong with the inter-layer connections constitute the mul-tilayer G . In each of the layers, the network is built usinga Molloy-Reed configurational algorithm with γ = 2 . q of the N nodes ofthe two layers with an intensity p . Three proceduresto draw inter-layer connections are considered dependingon the inlayer number of connection of the nodes (de-gree) (see Figure 1). The first one is simply uncorrelated(UnC) regardless of the nodes inlayer degrees. This isthe baseline scenario and has been profusely used in theliterature. The other two procedures include some corre-lation between the degrees of the nodes connected acrosslayers. The positively-correlated (PC) method preferen-tially connects high degree nodes in one layer G withtheir counterparts in the other layer G . Conversely, thenegatively-correlated (NC) procedure establishes linksbetween high degree nodes in G and low degree nodesin G and vice versa. Within the different categories ofcorrelations, the inter-layer links can be drawn at ran-dom or following some order based on the nodes inlayerdegrees. The possible sequence types corresponding toeach layer correlation is shown in Fig. 1. For a PC multi-plex, an ascending sequence connects node-pairs startingfrom those that have the lowest degree, the descendingsequence starts from the node-pairs with the highest de-gree. In a NC multiplex the ordered sequence establishesconnections between high-degree nodes on one layer andthe low-degree nodes in another first. Thus, to create a UnCPCNC λ p a) q = 1 Random p , q σ b) Figure 2. λ for a N = 100 scale-free multiplex. In a), meanand standard deviation over realizations (shaded region) of λ of M ( p, q = 1) as a function of p with positive (PC), nega-tive (NC), and neutral (UnC) inter-layer correlations, for 100realizations. In b), average λ for the PC case of a) (labelled‘intensity-based’) as a function of p , as well as the trends of λ as a function of q for multiplexes M ( p = 1 , q ) with PClayers and three different connection sequences, all for 10 re-alizations. Inset: standard deviation of the respective curvesacross the realizations. multiplex M ( p, q ) we must specify a) the size N of eachlayer, b) the layer correlation strategy (UnC, PC or NC),c) the strength of inter-layer connections p and the frac-tion of connected nodes q , and d) the sequence type thatgives the order in which layers become interconnected.For each realization of M ( p, q ), we compute the alge-braic connectivity λ of the supralaplacian L . The pro-cedure for constructing the supralaplacian of a multiplexwith q = 1 is detailed in [33], yielding L = (cid:18) L + p − p − p L + p (cid:19) , (1)where L k stands for the Laplacian of the layer G k , sep-arately, and is the unit N × N matrix. In order toimplement a situation in which some nodes are not con-nected across layers ( q (cid:54) = 1), we assign to such nodes,e.g. node i in G , a partner down in G (arbitrarily, thenode with the same index i ), but set the correspondinginter-layer link intensity equal to zero, p i = 0. Thus, amultiplex M ( p, q ) will have a supralaplacian where theelements in the matrices p corresponding to positions i are zero instead of p .The behavior of λ for different multiplex structures M ( p, q ) is shown in Figure 2 as a function of the inter-layer link intensity p . In Figure 2a, the multiplex is com-plete and we recover the results of [33]. The algebraicconnectivity displays an inflection point for a particularvalue of p , p x , around 0 . p < p x to an integrated network one for larger p . The correla-tions do play a role in the second integrated regime wherethey influence the relaxation time of the system (slowerfor PC and faster for NC with respect to the uncorrelatedbaseline). Increased system size N only magnifies thesedifferences. Interestingly, after the inflection point thefluctuations in λ between realizations of the multiplexincrease and remain constant with increasing p . The par-ticular position of p x displays a dependency on the sys-tem size in the range explored here numerically, but thepresence of this regime of high λ variability across real-izations always appears. It is worth noting as well thatin contrast to ordinary phase transitions the increase in σ is not constrained to the neighborhood of p x .If the intensity of the inter-layer connections is keptconstant, p = 1, and q is varied instead, these curves no-tably change and their shape depends on the particularorder implemented in the inter-layer node connection tobuilt the multiplex (see Figure 2b). λ is lower for a mul-tiplex with fewer inter-layer links (with heavier weightof unity) than for a multiplex with fully interconnectedlayers but at small intensities. That is to some extent in-tuitive, as the missing inter-layer connections might havebeen much more necessary for the ‘emergent’ multiplex,and having them present albeit at a small intensity wouldconnect the layer more strongly. In addition, the threepossible sequences behave differently depending on q . Forsmall q , until about q < q = p x , they increase linearlywith q . The descending sequence results in a consistentlygreater λ , implying a faster dynamics on the multiplex(this seems to be a generic result [41]). This is manifestlynot the case at high q where the situation reverses. Be-sides the average λ , σ displays a similar behavior as inthe case of q = 1: beyond p x it increases and remainshigh afterward, although with peculiarities due to thesequence and correlations of the inter-layer connections(see inset of Figure 2b).On single multiplex realizations instead of on averageproperties, we find what is the most important result ofthis work. The multiplex M ( p = 1 , q ) from the previ-ous section is constructed by adding inter-layer links ofweight unit one by one until their fraction is equal to q . This process accumulates microscopic changes to endup in a macroscopic configuration of the multilayer net-works. Naively, one could expect one of such microscopicchanges to be innocuous to the macroscopic picture and,consequently, to have a very minor impact on the dynam-ics of any process taking place on the multiplex. However,this expectation is wrong as can be seen in Figure 3a and3b on an example with a positively correlated multiplex.The value of λ experiences significant jumps after theintroduction of a single link across layers. The smoothbehavior characteristic of the absence of jumps in λ isonly visible in a sparsely interconnected multiplex (lower q ), where for some range (but depending on the sequence)the increase is linear, λ ( q ) ∝ q . Although the frequencyand location of λ ( q ) discontinuities depend on the spe-cific realization, on average they happen at some large q .Since λ is a global feature of the multiplex, the ‘offend-ing’ link is intrinsically connected to some global prop-erty of the layers, or to the sequences themselves. This . . . . . . λ . . . . . . q qPC ascending PC descendinga) b) sP(s) c) . . . . . . . . . . . . . . DescendingAscendingRandom Order P ε (q) . . . . . . q qd) e)ε = 0.02 ε = 0.001 DAR
P(s)P(s) s
Figure 3. Partial Connectivity Sequences in positively-correlated SF multiplex with N = 100, p = 1. In a) andb), λ as a function of q for the ascending and descendingdegree sequences. The plots contain ten realization in totalalthough the curves are displayed realization by realization.In c), probability density distribution of the differences be-tween consecutive values of λ as a new link is introducedfor the three difference sequence types D (descending), A (as-cending) and R (random), computed over 10 realizations at q = 0 . q = 0 . q = 0 . P (cid:15) to seejumps higher than (cid:15) , computed at each q value over 1000 re-alizations, for each of the ascending, descending and neutralsequences. reflects the precarious nature of adding final inter-layerconnections, implying that after some minimal number ofconnections has been set, it becomes very hard to predictthe exact effect of the addition of each inter-layer link.The distribution of jumps is included in Figure 3c,where beyond the sequence the width of the distribu-tion is controlled by the final value of q . Varying thisparameter, one can observe long tailed distributions inthe changes of λ spanning several orders of magnitudeas in the main plot or in the top-right inset, or more con- Interlayer connectivity q a) b)
Time -6.8-6.6-6.4-6.2-6 q = 0.3q = 0.47q = 0.49q = 0.77q = 0.78q = 0.84 c) λ (1/2) d ln(1-R)/dt l n ( - R ( t )) Figure 4. Validation of the effects of the λ jumps on the timescales of processes running on a PC multiplex M ( p = 1 , q ) with N = 100 generated with a sequence descending in degree. In a), λ and the numeric derivative of ln (1 − R ( t )) as a function of q .In b), a graphical representation of the multiplex M ( p = 1 , q ∗ ) in concentric circles, with nodes arranged counterclockwise andtheir size and location proportional to their degree, where q ∗ is from the set of lines in panel a) with the colors corresponding.In c), the temporal relaxation of the order parameter ln(1 − R ( t )) for each q ∗ . strained jump distributions as those of the left-bottominset for smaller values of q . For most of the values of q ,the dominating effect is the non trivial long-tailed distri-butions. To explore this further, we introduce a cutoff (cid:15) and measure the probability P (cid:15) ( q ) to have an increasein λ by at least (cid:15) at the addition of the subsequentlink. At intermediate (cid:15) values (Figure 3d), it identi-fies that the main difference between the ascending andthe descending sequence is not that the former has fewerjumps (which could have been expected from the lowerstandard deviation), but that the jumps are centered atsome q , with a well-defined mean. The descending se-quence and the random one, on the other hand, show acompletely different profile of a monotonically increasingjump probability. If the last node pairs to be intercon-nected are those with the lowest degrees, their additionis much more likely to case a sudden jump in the speedsof the diffusion, than if the last interconnected nodes hadhigh degree. That, however, is because in the latter casethe multiplex would already have been relatively fast.The reason lies in another curious observation shown inFigure 3e. The probability of at least a minute jump isminimized at some q < λ implies non-linear operations and onemight wonder if these results were not an artifact of thenumerical methods. To verify their relevance for dynam-ical processes occurring on the multiplexes, we run theKuramoto model and compare the variations in λ withthe synchronization timescales. Specifically, an oscillatoris placed in each node i of the multiplex with phase θ i .The oscillators have the same natural frequency and theevolution of their phases is described after linarization by the following equation [42, 43]˙ θ j = (cid:88) k (cid:54) = j W jk ( θ k − θ j ) , (2)where W jk are the elements of the adjacency matrix withthe corresponding weights for intra and inter-layer links.To ensure an easier measurement of the timescales, theinitial state of the system is set at the value of the eigen-vector associated with λ . Under these conditions, thesystem approaches synchronization as1 − R ( t ) ∝ e − λ t , (3)where R is the phase coherence of the population, R = | (cid:80) j e iθ j | , and plays the order parameter in the synchro-nization transition.In Figure 4, we show how the changes in the spectralgap λ actually translate to changes in the relaxationtimes. In Figure 4a, λ and the inverse relaxation timesare displayed as a function of q for a PC multiplex builtin descending degree sequence. Both coincide as it mustbe if the numeric estimation of λ is appropriate. Sev-eral jumps can be observed and are marked with verticallines of different colors and textures. These changes oc-cur after the introduction of single links as representedin Figure 4b. The presence of these single links bringsabout an important variation in the macroscopic systemrelaxation as can be seen in Figure 4c. The underlyingexponential decays for a select range of q values chosen tolie around two significant jumps. The observed straightlines imply that α is well-defined as an exponent even inthe vicinity of the jumps, and hence that the decrease intimescales is indeed abrupt.Our results thus show that the jumps in the algebraicconnectivity are not merely a numeric artifact, but in-stead correspond to the measured abrupt decreases inthe synchronization timescales of the dynamical system.This ratifies that microscopic changes in the topologi-cal properties of these multilayer networks lead to effectsnoticed at a global scale. We also investigated the char-acteristics of these offending node-pairs whose intercon-nection causes the jumps, but did not detect any speciallocal network feature. The lack of any correspondencebetween the degree, and the degree of the neighbors, ofthe select nodes points to the nontrivial nature of the con-nection between the network properties of the multiplexnodes and their role as crucial mediators in enhancing thesynchronization time of the multiplex, and calls for theuse of global spectral methods to determine their locationand quantify the possible effects of their introduction ordeletion.Partial financial support has been received fromthe Spanish Ministry of Economy (MINECO) andFEDER (EU) under the project ESOTECOS (FIS2015-63628-C2-2-R), from the EPSRC project GALE,EP/K020633/1 and from the EU Commission throughproject LASAGNE (318132). JJR acknowledges fundingfrom the Ram´on y Cajal program of MINECO. [1] M. Bogu˜n´a, R. Pastor-Satorras, and A. Vespignani, in S tatistical Mechanics of Complex Networks, LectureNotes in Physics, Vol. 625, edited by R. Pastor-Satorras,M. Rubi, and A. D´ıaz-Guilera (Springer, Heidelberg,2003) pp. 127–147.[2] A. Jamakovic, R. Kooij, P. van Mieghem, and E.van Dam, Proceedings of the 13th Annual Symposiumon Communications and Vehicular Technology in theBenelux , 35 (2006).[3] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes,Rev. Mod. Phys. , 1275 (2008).[4] F. M. Atay, T. Bitikoglu, and J. Jost, Physica D: Non-linear Phenomena , 35 (2006).[5] A. Arenas, A. D´ıaz-Guilera, and C. J. P´erez-Vicente,Physica D: Nonlinear Phenomena , 27 (2006).[6] J. G´omez-Garde˜nes, Y. Moreno, and A. Arenas, Phys.Rev. Lett. , 034101 (2007).[7] J. Chen, J.-A. Lu, C. Zhan, and G. Chen, in Handbook ofOptimization in Complex Networks , Springer Optimiza-tion and Its Applications, Vol. 57, edited by M. T. Thaiand P. M. Pardalos (Springer US, 2012) pp. 81–113.[8] J. A. Almendral and A. D´ıaz-Guilera, New Journal ofPhysics , 187 (2007).[9] E. Estrada, E. Vargas-Estrada, and H. Ando, Phys. Rev.E , 052809 (2015).[10] J. Mart´ın-Hern´andez, H. Wang, P. Van Mieghem, andG. D’Agostino, Physica A: Statistical Mechanics and itsApplications , 92 (2014).[11] L. Lov´asz, Tech. Report YALEU/DCS/TR-1029, De-partment of Computer Science, Yale University (1994).[12] S. Boccaletti, G. Bianconi, R. Criado, C. del Genio, J.G´omez-Garde˜nes, I. Sendi˜na-Nadal, Z. Wang, and M.Zanin, Physics Reports , 1 (2014).[13] M. Kivel¨a, A. Arenas, M. Barthelemy, J. P. Gleeson, Y.Moreno, and M. A. Porter, Journal of Complex Networks , 203 (2014).[14] K.-M. Lee, B. Min, and K.-I. Goh, The European Phys-ical Journal B , 1 (2015).[15] G. di Iasio, S. Battiston, L. Infante, and F. Pierobon,Munich Personal RePEc Archive, paper 52141, availableat https://mpra.ub.uni-muenchen.de/52141 (2013).[16] S. Poledna, J. L. Molina-Borboa, S. Mart´ınez-Jaramillo,M. van der Leij, and S. Thurner, Journal of FinancialStability , 70 (2015).[17] R. Burkholz, M. V. Leduc, A. Garas, and F. Schweitzer,available at http://arxiv.org/abs/arXiv:1506.06664(2015).[18] M. De Domenico, A. Sol´e-Ribalta, S. G´omez, and A. Are-nas, Procs. Natl. Acad. Sci. U.S.A. , 8351 (2014).[19] E. Strano, S. Shai, S. Dobson, and M. Barthelemy, Jour-nal of The Royal Society Interface , 20150651 (2015).[20] C. Granell, S. G´omez, and A. Arenas, Phys. Rev. Lett. , 128701 (2013).[21] F. Vazquez, M. A. Serrano, and M. San Miguel, availableat http://arxiv.org/abs/arXiv:1511.05606 (2015).[22] P. J. Mucha, T. Richardson, K. Macon, M. A. Porter,and J.-P. Onnela, Science (2010).[23] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, andS. Havlin, Nature , 1025 (2010).[24] P. Klimek and S. Thurner, New Journal of Physics ,063008 (2013).[25] P. Csermely, Talent Dev. and Excellence , 115 (2013).[26] Z. Wang and A. S. M. Perc, Sci. Rep. , 2470 (2013).[27] M. Diakonova, M. San Miguel, and V. M. Egu´ıluz, Phys.Rev. E , 062818 (2014).[28] P. Klimek, M. Diakonova, V. M. Egu´ıluz,M. San Miguel, and S. Thurner, available athttp://arxiv.org/abs/arXiv:1601.01576 (2016).[29] M. Diakonova, V. Nicosia, V. Latora, and M. San Miguel,New J. Phys. , 023010 (2016).[30] S. G´omez, A. D´ıaz-Guilera, J. G´omez-Garde˜nes, C. J.P´erez-Vicente, Y. Moreno, and A. Arenas, Phys. Rev.Lett. , 028701 (2013).[31] A. Sol´e-Ribalta, M. De Domenico, N. E. Kouvaris, A.D´ıaz-Guilera, S. G´omez, and A. Arenas, Phys. Rev. E , 032807 (2013).[32] F. Darabi Sahneh, C. Scoglio, and P. Van Mieghem,Phys. Rev. E , 040801 (2015).[33] F. Radicchi and A. Arenas, Nature Physics , 717 (2013).[34] I. L. Juan P. Garrahan, available athttp://arxiv.org/abs/arXiv:1406.4706 (2014).[35] N. Bastas, F. Lazaridis, P. Argyrakis, and M. Maragakis,EPL (Europhysics Letters) , 38006 (2015).[36] H. Shakeri, N. Albin, F. Darabi Sahneh, P. Poggi-Corradini, and C. Scoglio, Phys. Rev. E , 030301(2016).[37] L. V. Gambuzza, M. Frasca, and J. G´omez-Garde˜nes,EPL (Europhysics Letters) , 20010 (2015).[38] F. Radicchi, Phys. Rev. X , 021014 (2014).[39] V. Nicosia and V. Latora, Phys. Rev. E , 032805(2015).[40] M. Molloy and B. Reed, Random Structures & Algo-rithms , 161 (1995).[41] J. Aguirre, R. Sevilla-Escoboza, R. Guti´errez, D. Papo,and J. M. Buld´u, Phys. Rev. Lett. , 248701 (2014).[42] Y. Kuramoto, C hemical Oscillations, Waves, and Turbu-lence (Springer–Verlag, New York, 1984).[43] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, andC. Zhou, Physics Reports469