Diffusion geometry of multiplex and interdependent systems
DDiffusion Geometry of Multiplex and Interdependent Systems
Giulia Bertagnolli ∗ CoMuNe Lab, Fondazione Bruno Kessler, Via Sommarive 18, 38123 Povo (TN), Italy andUniversity of Trento, Department of Mathematica, Via Sommarive, 14, 38123 Povo (TN), Italy
Manlio De Domenico † CoMuNe Lab, Fondazione Bruno Kessler, Via Sommarive 18, 38123 Povo (TN), Italy
Complex networks are characterized by latent geometries induced by their topology or by thedynamics on the top of them. In the latter case, different network-driven processes induce distinctgeometric features that can be captured by adequate metrics. Random walks, a proxy for a broadspectrum of processes, from simple contagion to metastable synchronization and consensus, havebeen recently used in [Phys. Rev. Lett. 118, 168301 (2017)] to define the class of diffusion geometryand pinpoint the functional mesoscale organization of complex networks from a genuine geometricperspective. Here, we firstly extend this class to families of distinct random walk dynamics –including local and non-local information – on the top of multilayer networks – a paradigm forbiological, neural, social, transportation, biological and financial systems – overcoming limitationssuch as the presence of isolated nodes and disconnected components, typical of real-world networks.Secondly, we characterize the multilayer diffusion geometry of synthetic and empirical systems,highlighting the role played by different random search dynamics in shaping the geometric featuresof the corresponding diffusion manifolds.
I. INTRODUCTION
Complex networks, an abstract representation of thestructural and functional backbone of complex systems,exhibit a wide spectrum of geometric features, from self-similarity to latent hidden metric spaces and topology [1–4], which have been successfully exploited to gain newinsights about the structural and dynamics of social [5],neural [6], transportation [7, 8] and communication sys-tems [9], to mention a few emblematic examples (seeRef. [10] for a review). More recently, it has been shownthat even the dynamics on the top of complex networkscan induce complex geometries which can not be under-stood from inspecting only structural ones. Such ge-ometries have been studied in the case of a broad classof spreading phenomena: communicability distance hasbeen introduced to study information exchange based onwalks [11]; an effective distance has been introduced forcontagion processes [12]; a diffusion distance has been in-troduced to study collective phenomena such as synchro-nization, consensus and random searches [13]; a tempo-ral distance has been defined to study how perturbationspread in biological systems [14]. The common rationaleis to model the propagation of information with networkdynamics and investigate the corresponding distances in-duced between pairs of nodes.Random walks are emblematic examples used for mod-ellng diffusion and transport dynamics from lattices todisordered media and quantum systems, across morethan a century [15–21]. They provide both an intuitive– and often analytically treatable – mathematical frame-work and a rich physical model that can be used to map ∗ Corresponding author: [email protected] † Corresponding author: [email protected] a wide spectrum of random processes and observed fea-tures, especially in the contest of complex networks wherethe interplay between the underlying topology and thedynamics on the top of it is responsible for system’s func-tion [22]. The properties of the random walks reflect par-ticular structural features of a system: for instance, theycan be used to unravel the mesoscale organization of sys-tems from a functional perspective [23–27]; localizationeffects can be observed around topological defects [28];the average number of different sites visited in a timeinterval by a random walker, the so-called coverage, canbe used to quantify the system resilience to random fail-ures [29]; the relative importance of system’s units canbe quantified in terms of their ability to attract the over-all flow, encoded by random walkers [30]; the dynamicalfeatures of the information flowing through a networkcan be understood [31], and even enhanced [32], in termsof the spectral entropy defined by statistical features ofrandom walks.Based on these dynamics, one can also define similar-ity measures that reflect the ability of the units to ex-change information [13, 33–35]. In [13] diffusion distancewas defined for single-layer networks and it has provento be useful for characterizing the functional structure ofcomplex networks, e.g. identifying functional clusters, orcentral nodes [36]. It also provides a continuous-valueddistance function which, in applications where the globalinformation about the underlying topology is missing orpartial, is more informative than the purely topologicallength of shortest-paths.However, in many real systems entities interact in mul-tiple ways, think, for example, to locations connectedthrough different transportation modes or to metabolitesin a biological network linked by various types of chem-ical reactions. The evolution of the traditional graph-theoretical tools allows us, nowadays, to represent this a r X i v : . [ phy s i c s . s o c - ph ] J un multifaceted information by means of layers defining mul-tilayer networks [37, 38]. Figure 1 shows distinct repre-sentations of the multiplexity and the interdependencyobserved in empirical complex networks. From a math-ematical perspective, taking into consideration the vari-ety of connectivity patterns and coupling among layersand dynamics [39] requires multilinear algebra and ten-sors [37] to be efficiently represented and treated analyt-ically. αβ M i α j α αβαβ D(i; α , β ) ii M i α j β M i β j β FIG. 1. Different multilayer networks and their supra-adjacency matrix representation [37, 40]: (a) an edge-coloredmultigraph, with layers corresponding to colors and no inter-layer connections; (b) a multiplex network where the replicasin the different layers are interconnected sequentially, a typeof intertwining or diagonal coupling, and (c) the most generalinterconnected case, where inter-layer connections are not re-stricted to replicas (exogenous interactions). Latin letters de-note nodes, while Greek letters are used for layers. D αβ ( i ) isthe intensity of the intertwining between state nodes i, α and i, β . C ij ( αβ ) describes the general inter-layer connectivity. In this study our contribution is twofold: on the onehand, we generalize multilayer random walk dynamicsfor the analysis of more realistic empirical systems, e.g.allowing for the existence of isolated nodes and compo-nents; on the other hand, we extend the framework ofdiffusion geometry [13] to realm of multilayer systems,highlighting the dependence of the diffusion manifold onthe interplay between structure and dynamics within andacross layers.The paper is organized as follows: in the next sectionwe recall the basic facts on random walks and diffusionprocesses on single layer networks and then generalizethem to the multilayer framework. In Sec. III we extendthe definition of diffusion distance to multilayer networks,presenting results obtained from the analysis of syntheticmodels. We conclude with applications to empirical so-cial and transportation systems in Sec. IV, followed byDiscussion and Conclusions.
II. DIFFUSION DYNAMICS IN COMPLEXNETWORKS
In a group of individuals a piece of news spreads dif-ferently if everyone knows each other or if there are sub-groups of individuals that do not communicate with theothers: the structure of this social network influences thediffusion of information. Similarly, one observes a stronginfluence of the topology on the dynamics of a broadclass of empirical systems. Because of this influence ofthe structure on the evolution of dynamical processes,the latter are often used as a proxy for probing and un-covering structural properties of complex networks. Asimple strategy to explore a network is to start from anode, to choose at random an (outgoing) edge and tofollow it toward the next node, i.e., to perform a ran-dom walk on the network with transition probabilitiesprescribed by the network connectivity. Random walks(RWs) represent a useful model for diffusion processes.Here we will focus on node-centric continuous-time ran-dom walks (CTRW)[22], continuous-time Markov chains(MC) on the vertex set of a network with transition ratesdepending on the network structure and specific naviga-tion rules.
A. Random walks and diffusion distance onsingle-layer networks
Let us consider a weighted directed network G =( V, E ) without isolated nodes and with (possiblyweighted) adjacency matrix W = { W ij } , and two nodes i, j ∈ V . In a discrete random walk the transition prob-ability from i to j in one time step is given by p ij = W ij s i .We can write it in matrix form p ( n + 1) = p ( n ) D − W ,where p is seen as a row vector and D is the diago-nal matrix of out-strengths D ii = s i = (cid:80) j W ij . Thecontinuous-time random walk corresponding to this jumpchain (also-called its embedded Markov chain) is de-scribed by the forward equation (cid:40) ˙ p ( t ) = − p ( t )˜ Lp (0) = p (1)where ˜ L = I − D − W is the random walk normalizedLaplacian and − ˜ L is the generator of the continuous-timeMarkov chain with initial distribution p ; see Appx. Afor further details. The transition probabilities are givenby the solution of (1), i.e., p ij ( t ) = ( p e − t ˜ L ) ij . Rewriting(1) as a system with P ( t ) being a matrix and with the ini-tial distribution given by the identity matrix, P (0) = I ,we obtain the unique solution P ( t ) = e − t ˜ L . Its i − th row (cid:16) e − t ˜ L (cid:17) i · = p ( t | i ) is the probability vector correspond-ing to p = e i and (cid:16) e − t ˜ L (cid:17) ij is the probability of beingin j starting in i with probability 1, after time t . The α =2 second conn. comp.LCC isolated node α =1 11 (a)
10 20 30 40 50 60102030405060 (b)
10 30 50103050 10 30 50 10 30 50 00.330.650.981.3distance c l a ss i c a l ( C R W ) d i ff u s i v e ( D R W ) p h y s i c a l w / r e l a x a t i o n ( P r R W ) (c) t = 1 t = 2 t = 5 FIG. 2. A synthetic two-layers network with unitary diagonal coupling D ( i ; αβ ) = 1, for all i and α (cid:54) = β , and its functionalcharacterization. (a) Isolated nodes and different components have been highlighted in the two layers. (b) The supra-adjacencymatrix of the multiplex. (c) the N × L diffusion supra-distance matrices for three different random walk dynamics and t = 1 , , r = 0 . diffusion distance [13] is the L − norm of the differencebetween rows of the matrix e − t ˜ L D t ( i, j ) = (cid:107) p ( t | i ) − p ( t | j ) (cid:107) . (2)Two nodes i (cid:54) = j are “near” w.r.t. D t if the probabilitythat two random walkers starting in i and j respectively,meet somewhere in the network after time t .Upon a given network, which is completely character-ized by its adjacency matrix W , we can define randomwalks with different flavors. This enables us, not onlyto model a wider range of physical spreading processes,but also to remove the rather restrictive assumption ofconnectedness, which is necessary for writing the transi-tion matrix D − W of the classical random walk. Con-sequently, also the family of diffusion distances can begeneralized to different types of random walks, as wellas to different types of networked systems. In the re-mainder of this section, we extend the diffusion distanceto multilayer networks, additionally exploring the vary-ing patterns induced by different random walks dynamics[27, 29]. Of course, there are many other types of ran-dom walks that are not mentioned here, e.g., multiplica-tive processes [41] or correlated random walks [42], and are often defined with a specific motivation. Neverthe-less, our framework is very general and, given a transi-tion matrix T , can be effortlessly extend considering thecontinuous-time Markov chain having exponential hold-ing times with rate one and T as jumping matrix, seeAppx. A.Before moving to the multilayer case, let us introducethe tensorial notation [37] for the monoplex G , whichallows us to generalize the formalism to multilayer net-works. The adjacency matrix W can be seen as a rank-2tensor W ij . In this paper, we will not use the covari-ant notation and the Einstein summation convention, sothat W ij denotes the component ( i, j ) of the tensor. Themaster equation (1) can be re-written as ˙ p j ( t ) = − N (cid:80) i =1 ˜ L ij p i ( t ) p i (0) = q i (3)with ˜ L ij = δ ij − T ij indicating the component ( i, j ) of therandom walk normalized Laplacian tensor, δ ij the Kro-necker delta and T ij the component ( i, j ) of the transitionprobability tensor. CRW PRRW DRW MERW PrRW T jβjβ M jβjβ s j ( β )+ S j ( β ) r M jβjβ s j ( β )+ S j ( β ) + (1 − r ) NL s max + M jβjβ − ( s j ( β )+ S j ( β ) ) s max M jβjβ λ max (1 − r ) M jβjβ s j ( β ) + r M jβjβ s i T jαjβ M jαjβ s j ( α )+ S j ( α ) r M jαjβ s j ( α )+ S j ( α ) + (1 − r ) NL M jαjβ s max M jαjβ λ max V jβ V jα r M jβjβ s i T iβjβ M iβjβ s i ( β )+ S i ( β ) r M iβjβ s i ( β )+ S i ( β ) + (1 − r ) NL M iβjβ s max M iβjβ λ max V jβ V iβ (1 − r ) M iβjβ s i ( β ) + r M iβjβ s i T iαjβ M iαjβ s i ( α )+ S i ( α ) r M iαjβ s i ( α )+ S i ( α ) + (1 − r ) NL M iαjβ s max M iαjβ λ max V jβ V iα r M iβjβ s i TABLE I. Transition probabilities for different random walks. (CRW) classical, (PRRW) PageRank, (DRW) diffusive, (MERW)maximal-entropy, and (PrRW) pysical with relaxation random walks. s max = max i,α { s i ( α ) + S i ( α ) } ; the jumping parameter ofthe PageRank RW is commonly indicated by α , to avoid ambiguity, we indicate it by r ; λ max is the largest eigenvalue of theadjacency tensor and V is its corresponding eigentensor, satisfying (cid:80) i,α M iαjβ V iα = λ max V jβ (see [37] for details). Note that forCRW, DRW and MERW, these transition rules generalize the one introduced in [29] for the analysis of multiplex networks.The PrRW is defined as in [27], while PRRW generalizes the walk introduced in [43]. B. Random Walks on edge-colored multigraphs
A multilayer network is defined by a set of N nodes V interacting with each other in multiple ways, simultane-ously. The different types of interaction can be encodedby colors and grouped into layers. Each node i ∈ V ex-ists at least in one layer and, if it exists in multiple layers α ∈ { , . . . , L } , we usually referred to { i, α } as a replica or state node , in contraposition to the physical node i .Interactions of the same color determine the intra-layerconnectivity, while the others are called inter-layer con-nections. The simplest multilayer structure, depicted inFig. 1(a), is the edge-colored multigraph, where there isno additional information (e.g., no order relation) on theset of layers { , . . . , L } or, equivalently, inter-layer con-nectivity is not present.There are, essentially, two ways to define random walkson these networks. One possibility is to allow the ran-dom walker to follow a sequence of edges with differentcolors: colored edges are then treated as multiple edgesand the degrees of a node counts all the edges, regardlessof their colors [44]. This choice is equivalent to aggre-gate the multilayer system to a single layer and performa classical random walk on it. This approach is not de-sirable in general, because one does not know a priori if,and to which extent, information lost while aggregatingthe structure of layers will affect the results.The second approach, used successfully in other appli-cations [32], is to run independent dynamics on each layer(color) and then integrate the transition matrices over thelayers, properly normalizing each transition probabilityas follows: (cid:104) T ij (cid:105) = L (cid:88) α =1 µ i L T iαjα . (4)Indicating by s i ( α ) the out-strength of vertex i in layer α , µ i L (cid:80) α =1 { si ( α ) (cid:54) =0 } L represents the multiplicity of node i , i.e.,the fraction of layers in which i is not isolated (a trappingnode). This approach is more desirable than the first one,since it preserves more information related to diversity of connectivity patterns across layers. Nevertheless, thechoice of adding transition matrices as in Eq. (4) can bereplaced by a more general linear combination weightedby the relative importance given to each layer.The main difference between the two approachesemerges when the edges are weighted, besides colored: a priori the scale (nominal, ordinal, ratio etc.) of theweights could be different and an inattentive summationcould lead to errors. In (4), instead, for each i ∈ V , weare taking a finite mixture of probability mass functionswith uniform weights µ i L ≥ L (cid:80) α =1 1 µ i L = 1;nothing prevents one to weight layers differently, if ad-ditional information is available. In the following, weconsider equal importance for all layers, thus avoidingthe dependence of our results on the choice of a specificset of weights to be assigned to layers.We now take advantage of the more involved notationof Sec. II A, to introduce multilayer interconnected net-works. C. Random walks on multilayer networks
In the same way a monoplex can be represented by arank-2 tensor, a multilayer network is completely char-acterized by its rank-4 adjacency tensor [37]. Let us in-dicate by M iαjβ the components of the (weighted) multi-layer adjacency tensor. In this index notation M iαjβ indi-cates the interaction between i and j in the same layer,while M iαiβ denotes the strength of the intertwining ofthe replicas of i in two distinct layers. Depending onthe inter-layer connectivity we can have different typesof multilayers, as shown in Fig. 1. When inter-layer con-nections occur only between replicas of the same physicalnode, i.e., M iαjβ = 0 for i (cid:54) = j , the network is called amultiplex. The term M iαiβ encodes the intertwining be-tween two replicas; it is a scalar depending, in general,on i, α, β and we indicate it by D ( i ; α, β ). For the easeof visualization we flatten the adjacency tensor into theso-called supra-adjacency matrix [37, 40], of dimension N L × N L . Its diagonal blocks are the adjacency matricesof the monoplexes, while its off-diagonal blocks containthe information on the inter-layer connectivity.We henceforth indicate by s i ( α ) = N (cid:80) j =1 M iαjα , the out-strength of node i in layer α and by s i = (cid:80) α s i ( α ) = L (cid:80) α =1 N (cid:80) i =1 M iαjα its multi-layer out-strength, discarding theinter-layer edges. The inter-layer strengths are obtainedas S i ( α ) = (cid:80) j (cid:80) β (cid:54) = α M iαjβ and, consequently, ( s i ( α ) + S i ( α )) N,Li =1 ,α =1 is the out-strength supra-vector with N L components obtained as the row sums of the supra-adjacency matrix.Generally, the presence of isolated nodes or compo-nents, as well as an heterogeneous or mixed distributionof inter-layer connectivity is discarded to facilitate theanalytical framework. Conversely, here we do not forceany assumption on the structure of the multilayer sys-tem: if some units are not present in all layer – whichis often the case for real data – we will add artificialreplicas, which will appear as isolates and we will ac-count for them adequately. We henceforth indicate by V the common vertex set and by N = | V | the number ofphysical nodes. Although the isolated nodes do not mod-ify the connectivity of the connected component a layer,they could be troublesome for the evaluation of transi-tion probabilities. We have to distinguish three cases (i) S i ( α ) = 0, (ii) s i ( α ) = 0, and (iii) s i = 0. The lattercorresponds to the trivial case, where node i is isolatedin every layer, so that i can be simply removed from thevertex set. We can assume, without loss of generality, s j > i ∈ V . (i) and (ii) – corresponding to nointer-layer and no intra-layer connections, respectively –constitute a problem only if S i ( α ) + s i ( α ) = 0 for some α . In this case one could see { i, α } as an absorbing statewith the probability of remaining there equal to 1; an-other option is to teleport the random walker in { i, α } to any { j, β } with uniform probability NL . As for otherapproaches [23, 27, 30], we opted for the second option,since it decreases the occupation probability of the statenode { i, α } .In the multilayer framework the probability transitionsin one time step constitute the components T iαjβ of a rank-4 tensor and we can expand the RW equation to highlightdifferent contributions of jumps and switches in the dy-namic p jβ ( t + 1) = T jβjβ p jβ ( t ) (cid:124) (cid:123)(cid:122) (cid:125) stay + L (cid:88) α =1 α (cid:54) = β T jαjβ p jα ( t ) (cid:124) (cid:123)(cid:122) (cid:125) switch ++ N (cid:88) i =1 i (cid:54) = j T iβjβ p iβ ( t ) (cid:124) (cid:123)(cid:122) (cid:125) jump + L (cid:88) α =1 α (cid:54) = β N (cid:88) i =1 i (cid:54) = j T iαjβ p iα ( t ) (cid:124) (cid:123)(cid:122) (cid:125) switch and jump α = 1 α = 2 isolated node FIG. 3. Average diffusion matrices of the two-layers syn-thetic network shown in Fig.2(a), w.r.t. three RW dynam-ics: (top) classical, (middle) diffusive, and (bottom) physicalwith relaxation r = 0 .
5. The state nodes { i, α } are coloredin the corresponding dendrogram according to the layer theybelong to, except for isolated nodes, highlighted in purple.Small squares with a black border indicate those state nodesgrouped together, e.g. { , } , { , } in the top panel. Its continuous-time version is described by the forwardequation ˙ p jβ ( t ) = − (cid:88) i,α ˜ L iαjβ p iα ( t ) (5)where ˜ L iαjβ = δ iαjβ − T iαjβ .Upon a given network structure we can define differenttypes of random walks, depending on the specific ruleswe want our random walker to explore the network. For aggregated network (b) equivalent distances (d)(c) edge-colored network (a) adjacency matrix aggregated network FIG. 4. Comparing different types of distances among physical nodes, regardless of layers. (a) The weighted adjacencymatrix of the aggregated network. The diffusion distance (w.r.t. the CRW and averaged over time) (b) of the aggregatednetwork corresponding to the synthetic multiplex of Fig. 2(a) and (c) of the edge-colored network. (d) The equivalent diffusiondistances obtained with the reduction of Eq. (7). The distances are rescaled to [0 , D t ) ≈ .
13, edge-colored: max( ¯ D t ) ≈ .
16, equivalent distances: max( ¯ D t ) ≈ . instance, in PageRank [30] a teleportation or jumping parameter gives the possibility to the walker to reachalso nodes that are not directly connected to the currentnode. The flavors of the random walks are given by theirtransition rules, which depend on the structure of thenetwork, as a function of the adjacency tensor. The RWpresented in the monoplex case is referred to as a clas-sical random walk (CRW). Additionally to the classicalrandom walk, we look here at a family of diffusion dis-tances based on four other random walk types: multilayerPageRank (PRRW) generalizing Refs. [30, 43], multilayerdiffusive (DRW) generalizing the one defined in Ref. [29],maximal-entropy (MERW) generalizing Refs. [28, 29],and physical random walk with relaxation (PrRW) [27],whose transition probabilities are shown in Tab.I.The physical random walk has been defined in [29] todescribe those dynamics where the state nodes have a“common memory”, so that the information diffuses in-stantaneously across replicas. Think, for instance, of thesystem of virtual interactions among individuals, whomay have a profile (an alter-ego) in different social net-works. A person can then exchange information in a par-ticular social network using the (intra-layer) connectionsof its alter-ego in that social system, but she/he has al-ways a complete knowledge of the information across thelayers. In this case, inter-layer connections between repli- cas of different physical nodes have no physical mean-ing and, consequently, are ignored. The physical ran-dom walk with relaxation (PrRW) [27] is a variant onthe physical random walk, where the assumption onthe complete knowledge of intertwining between layersis dropped. It can be seen from Tab.I that its transi-tion probabilities contain a trade-off between intra- andinter-links, which are followed with probability 1 − r and r respectively. If not differently stated, we consider here r = 0 . III. DIFFUSION DISTANCE IN MULTILAYERSYSTEMS
In the multilayer framework the probability of findinga random walker at a given node and layer is encoded ina time dependent tensor, whose component p jβ ( t ), corre-sponds to the probability of finding the random walker inthe state { j, β } , for a given initial distribution. Similarlyto (2), we can then define the diffusion distance betweenstate nodes { i, α } and { j, β } as D t ( { i, α } , { j, β } ) = (cid:88) k,γ ( p kγ ( t |{ i, α } ) − p kγ ( t |{ j, β } )) . (6)The diffusion distance is bounded in [0 ,
2] for all i ∈ V , α ∈ , . . . , L , and t >
0, indeed, (cid:80) k,γ ( p kγ ( t |{ i, α } )) ≤ (cid:80) k,γ p kγ ( t |{ i, α } ) = 1, and it is small if there is a largeprobability that the random walks starting in { i, α } and { j, β } meet somewhere in the multiplex by time t . Fur-thermore, as the diffusion time t increases, and assumingwalk is ergodic, p k,γ ( t ) will tend to the stationary distri-bution π k,γ and D t ( { i, α } , { i, β } ) → D t for thethree more diverse RW dynamics on a synthetic multi-layer network. In each layer, we have a network with N = 30 nodes generated from a stochastic block modelwith two blocks. We chose the probabilities in order tohave a diverse topology: dense groups, disconnected com-ponents and isolated nodes. The coupling between thelayers is D ( i ; 1 ,
2) = D ( i ; 2 ,
1) = 1, as shown in the supra-adjacency matrix of Fig. 2(b). The columns of panel (c)represent different diffusion times t . Recall that the dif-fusion time plays the role of a scale parameter [13] andthat the continuous-time Markov chain has exponentiallydistributed holding times with rate λ = 1, i.e., the ex-pected time occurring among each step of the RW is 1. D t =1 is then a function of the micro-scale structure of themultiplex and here the isolated nodes are clearly visible.Remarkably, the distances w.r.t. the diffusive randomwalk span a larger interval than the others and differen-tiate very effectively the subgroups in each layer. This isa consequence of the diffusive dynamic, where the jump-ing probabilities do not depend on the vertex, but onthe strength of hubs. It is also worth noticing that, thosenodes, which are disconnected from the largest connectedcomponent of the second layer, are generally nearer to thenodes in the first layer than to the nodes in their samelayer. To unveil the meso-scale structures of the net-work we averaged the diffusion distance matrices for upto t max = N [13] and run a hierarchical clustering on theaverage diffusion distance supra-matrices, summarizingthe results in Fig. 3. As expected, the presence of inter-links moves loosely connected nodes near to their replicain the other layer and the clustering does not separatethe two layers.Finally, we provide a grounded way to summarize thesupra-distance matrix into an N × N matrix collecting thediffusion distances among the physical nodes, regardlessof the layers. The ij − elements of the diagonal blocks ofthe supra-matrix, { D t ( { i, α } , { j, α } ) , α ∈ , . . . , L } , canbe seen, using the jargon of electrical circuits, as resis-tances in parallel between the physical nodes i, j . Theirequivalent resistance can then be found through the par- allel sum of the resistances between the replicas D t ( i, j ) = (cid:32) L (cid:88) α =1 D t ( { i, α } , { j, α } ) (cid:33) − . (7)Note that this is an effective distance obtained fromdiffusion distances across layers, and it is not related tothe concept of resistance distance [45].The resulting equivalent distance matrix is quite dif-ferent from the one obtained evaluating the diffusion dis-tance on the aggregated network, as shown in Fig. 4.Here the distance matrices have been rescaled in [0 ,
1] bydivision through their respective maximum values, re-ported in the figure caption. Looking back at Fig. 2(a),we can see that the nodes from 1 to 19 are densely con-nected in layer 1, while node 8 and nodes from 11 to30 form the largest connected component of the secondlayer. In both distance matrices there is a clear blockcorresponding to the last ten nodes. However, only theequivalent diffusion distance matrix captures the partic-ular position of the first ten nodes (removing 8): in themultiplex they are distant from the last ten, because theybelong to different communities in layer 1 and to discon-nected components in layer 2. We also compare the ag-gregated network distance matrix Fig. 4(b) to the oneof the edge-colored multigraph obtained removing theinter-layer links from the multilayer network, shown inFig. 4(c). The two matrices are very similar, with onlyminor permutations inside the clusters.
A. Multilayer Diffusion Manifolds
The use of different random walk dynamics to explore asystem has an impact on the distances between its unitsand, consequently, on how the units are distributed inthe induced diffusion spaces. Similarly, the diffusion timeshapes the pairwise distances, highlighting local featuresof complex network geometry on short time scales and itsmore persistent structures for large diffusion times. Inthe multilayer setting there is an additional level of com-plexity given by the inter-layer connections and by thelayer-layer correlations. To gain further insights, we gen-erate 3 distinct classes of synthetic multilayer networks,with system size N = 200, and analyze them through thelens of diffusion geometry.The first class consists of Barabasi-Albert scale-freenetworks [46] on each layer: we consider a linear prefer-ential attachment with 4 edges added by each new nodeduring the growth process, while setting at at 10% theedge overlapping across layers – defined in terms of thefraction of links which are present in both layers amongthe same pairs of nodes [47].The second class consists of Watts-Strogatz small-world network [48] on each layer, obtained by rewiringlattices with probability 0.2, where edge overlapping istuned similarly to the first class. classical PageRank di ff usive max-entropy physical w/ relaxation B a r a b a s i - A l b e r t W a tt s - S t r o g a t z G i r v a n - N e w m a n FIG. 5. Average diffusion distance ¯ D t on two-layers multiplexes with different topologies, for fixed values of global averageedge overlap (Barabasi-Albert and Watts-Strogatz) and partition overlap (Girvan-Newman) between layers. Dendrograms onthe right-hand side of each distance matrix represents the corresponding hierarchical clustering to highlight the meso-scaleorganization of the system, with color encoding the planted node assignment in each layer. Distance matrices have beenrescaled ¯ D t max ¯ D t . The Barabasi–Albert model is characterized by the presence of hubs, which are clearly recognizable in thesupra-distance matrix w.r.t. the diffusive random walk, despite the small overlap between layers. The Girvan-Newman two-layers multiplex has a meso-scale organized in strong communities, partly overlapping across layers. Note that at variance withedge overlapping, here partition overlapping is defined in terms of nodes belonging to the same group without requiring thosenodes to be connected by an overlapping edge. The third class consists of layers with strong meso-scale structure organized in 4 groups, like in a Girvan-Newman model [49] on each layer: the probability thattwo nodes within the same group are connected is 1,whereas cross-group connections are much sparser andpresent with probability 0.05. Group overlapping [47] –defined in terms of the fraction of nodes planted in thesame group on both layers – is again fixed at 10%.The role of layer-layer correlations and their interplaywith the distinct network topologies considered above issummarized in Fig. 5. As expected, there are relevantdifferences due to the type of random search dynamicsand to the topological features of the underlying topolo-gies. For instance, the diffusive walk for the Barabasi-Albert system leads to a high level of mixed pathwaysacross layers, resulting in nodes that do not preserve inthe diffusion space their original layer assignment. Forthe same walk, in the case of the Watts-Strogatz systemthe result is the opposite: nodes aggregate into func-tional clusters that are not well mixed up in the diffusionspace. A high amount of geometric mixing is also ob-served when the physical random walk with relaxationis used, as expected. Overall, it is not guaranteed thatmultilayer diffusion pathways layers favor the geometricmixing in the diffusion manifold: the result depends onthe type of dynamics and on layer-layer correlations. To gain additional insights, we have considered a sec-ond battery of synthetic models, where we increasinglyadd inter-layer connectivity between layers.The absence of information pathways across layers,happening for instance when two layers are not coupledtogether, leads naturally to disjoint diffusion manifolds,each one corresponding to the distinct layers. When thetwo layers are interconnected together, a trivial resultis that the strength of inter-layer connectivity facilitatethe flow of information across layers. However, the aboveprocess hides an interesting phenomenon, that is unveiledin Fig. 6. To better characterize it, we calculate theFrobenius norm, which is defined as follows, for a genericmatrix A (cid:107) A (cid:107) F = (cid:118)(cid:117)(cid:117)(cid:116) m (cid:88) i =1 n (cid:88) j =1 | A ij | = (cid:113) trace ( A T A )to quantify the overall intensity of an average diffusiondistance matrix. The Frobenius norm is computed whenthe two layers are not coupled, on the union of the twodistance matrices, i.e., (cid:112) (cid:107) A (cid:107) + (cid:107) B (cid:107) , and then on thesupra-distance matrices for increasing fraction of inter-layer connectivity: first state nodes corresponding to thesame physical node are interconnected with each otherto create an interconnected multiplex; after that this ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●● di ff usion time, t fraction of inter- and cross-links F r o b e n i u s n o r m o f D t (a) (b) (c) (d) FIG. 6. Frobenius norm of the supra-distance matrices of a synthetic two-layers network, w.r.t. a diffusive random walk.As the fraction of inter-layer links grows we move from two disconnected multiplexes to a fully inter-connected multi-layerwhit all N connections across layers. The heatmaps are four representatives of the different regimes: (a) uncoupled layers;(b) a single-inter-link between the replicas of a random nodes coupling the two layers, i.e., (cid:2) N , N (cid:3) partially interconnectedmultiplex; [1/N] fully interconnected multiplex (all state nodes corresponding to the same physical node are interconnected);(c) (cid:2) N , (cid:3) multilayer regime consisting of an interconnected multiplex with the addition of cross-links between state nodes ofdistinct physical nodes. regime is reached, the cross-links between state nodescorresponding to all other physical nodes are created,until the total of N connections is generated. Remark-ably, when one interlink is added between the layers, theFrobenius norm increases: this is due to the fact that thenew link coupling the two layers creates a bottleneck forinformation to be exchanged across layers, even for largevalues of τ . Once more inter-links are added, the Frobe-nius norm decreases, until it tends to reach a plateauwhen the fraction of inter-layer links is 100%.Our results highlight that the existence of topologicalcorrelations across layers induce changes in how infor-mation is exchanged between state nodes. Such changesalter diffusion distances and might lead to two differ-ent regimes: i) flow keeps segregated within layers andthe multilayer diffusion manifold consists of two well sep-arated sub-manifolds representing each layer separately,which are connected by weak geometric pathways; ii) flowis integrated, creating new geometric pathways for infor-mation that mix up those sub-manifolds. IV. APPLICATIONS TO EMPIRICALMULTILAYER SYSTEMS
We use the newly introduced family of metrics to studytwo real systems with multiple types of interactions: themultimodal transportation network of London [29] andthe multilayer Noordin Top terrorists network [50]. Thefirst system consists of three layers corresponding to theTube, overground, and DLR, arranged in a multiplexwith couplings D ( i ; α, β ) = 1 for α (cid:54) = β ∈ { , , } .Nodes represent stations ( N = 369 in total) and connec-tions between them are weighted and undirected. Thenetwork of interactions among 78 Indonesian terrorists( N = 79 in the data set, but actor 58 is usually re-moved since it is disconnected in all layers) is a four-layers multiplex, representing their pairwise trust (T),operational (O), communication (C) ties, and business(B) relations [51].Figure 7 shows the average diffusion manifolds, pro-jected in R through multidimensional scaling, induced0 L o n d o n t r a n s p o r t a t i o n N oo r d i n T o p t e rr o r i s t s FIG. 7. Average diffusion distance supra-matrices ¯ D t organized according to their hierarchical clustering as in previous analyses,and projection in R of the corresponding diffusion manifolds, in the case of the two real multiplex networks: the London publictransportation network (top panels) and the social relationships of Indonesian terrorists (bottom panels). In both cases, nodeshave been embedded in space, more specifically in [ − , , through multidimensional scaling. To facilitate the visualization ofthe three-dimensional embeddings, we draw the surface which better approximates the cloud of nodes in R and project it onthe plane, encoding the third dimension with colors. Nodes are shown as dots on the top of the surfaces. by different RW dynamics. As observed in [29], the bestexploration strategy, i.e. the RW to adopt to cover ef-ficiently the network, depends on the topology of themultilayer. This is reflected in the maps of Fig. 7, eventhough they are low-dimensional approximations of thetrue diffusion manifolds. As a matter of fact, for theLondon transportation network, the manifolds inducedby the classical, PageRank, and diffusive random walksappear qualitatively very similar with each other, andconsiderably different from those induced by MERW andPrRW. Instead, the supra-distance matrices and man-ifolds obtained for the terrorists network appear simi-lar in that all have a group of nodes with small pair-wise distances, and another group of nodes which aredistant from each other. To quantify more adequatelyhow diverse, or similar, the manifolds are, we comparetheir supra-distance matrices by means of Mantel’s test[52, 53], where the null hypothesis is that the pairwisedistances in one matrix are not monotonically related tothe corresponding distances in the second matrix, andshow the results of our test in Figs. 8-9. V. DISCUSSION AND CONCLUSIONS
We have considered different families of random walkdynamics, adequately extended to the realm of multilayernetworks, to introduce the multilayer diffusion geometry.Its classical counterpart, single-layer diffusion geometry,intimately relates metastable synchronization, consensusand random search dynamics, providing a novel frame-work for identifying functional clusters in complex net-works.While the framework and its validity remain the same,its natural generalization to multilayer networks – i.e.,systems consisting of multiple types of relationshipsamong their units – was missing. Here, we fill this gapand provide evidence that multilayer diffusion manifoldsencode information due to the interplay between of themultilayer structure and the dynamics on its top.From the analysis of synthetic networks with overlap-ping edges or groups across layers, we have found thatthe interplay between dynamics and topology cannot beeasily decoupled: e.g., the classical random walk reveals1 − − − − − PrRWMERWDRWPRRWCRW
FIG. 8. Estimating the similarity between diffusion mani-folds corresponding to different RW dynamics, evaluated onthe London multimodal transportation network. We calcu-late the Pearson’s correlation (encoded by size and color) be-tween the entries of pairs of supra-distance matrices (Mantel’sstatistic) which is then tested for significance by permutation(permutation test), with α = 0 . − − − − − PrRWMERWDRWPRRWCRW
FIG. 9. Same as in Fig. 8 for the Noordin Top terroristsnetwork. a manifold in which state nodes of the different layershave larger distances than the intra-layer distances, but this does not remain true in the presence of strong com-munities, which are not necessarily overlapping. Alsothe behavior of the MERW is not trivial: the top levelhierarchical structure unveiled in the Barabasi-Albert iscompatible with that of the Watts-Strogatz network, de-spite the high heterogeneity of the first, and this could besurprising since MERW is influenced by irregularities innodes degree. In another scenario, where two layers areoriginally uncoupled and do not exchange information,we add inter-layer connectivity to better understand howthe originally disjoint diffusion manifolds approach eachother because of the presence of multilayer informationpathways. Our results highlight that also in the regime ofpartial interconnected multiplex, where not all replicas ofa physical node are interacting – which in the real worldcould mean a failed connection between a bus and a trainstation – cross-layer pathways form, allowing the inter-layer information exchange. Furthermore, as we movetoward the fully-interconnected multilayer, distances be-come smaller (as shown by the Frobenius norm), but thedistinctive meso-scale structure of the Barabasi-Albertmodel, i.e., the presence of hubs, remains clearly visible.Finally, we have applied our novel framework to twoempirical multilayer systems, namely the public trans-portation of London and the social network of Noordinterrorists. The diffusion geometry corresponding to dif-ferent random walk dynamics are not necessarily distinct,and we have developed a quantitative method to assessthe correlation between the underlying multilayer diffu-sion manifolds. In the case of the transportation system,we find that the MERW, which in the synthetic networkswas able to separate the layers, highlights two groups ofnear nodes, that are not captured by other dynamics.This may suggest that (i) the structure of this system, atdifferent scales, has features that are characteristic of dif-ferent models (e.g., heterogeneity and communities); (ii)different dynamics induce different manifolds and conse-quently, the analysis of networked systems embedded intospace cannot exclude the analysis of the dynamics itself.Conversely, in the case of the social system, we find thatthe metrics have generally higher correlations, so thattheir latent diffusion spaces may be likewise similar. Thehierarchical structure unveiled by the supra-distance ma-trices seems to suggest a cross-layer core-periphery func-tional organization, which we leave to future work.Our work provides a novel tool for the analysisof multilayer systems from network geometry perspec-tive [10]. Since the latent diffusion geometry is inducedby network-driven processes, our framework provides alsoa complementary view to structural analysis, such as theone provided by hyperbolic network geometry [5, 54, 55],recently used to analyze multilayer networks [56], andhigher-order analysis [57]. [1] C. Song, S. Havlin, and H. A. Makse, Nature , 392(2005). [2] C. Song, S. Havlin, and H. A. Makse, Nature physics ,
275 (2006).[3] M. A. Serrano, D. Krioukov, and M. Bogun´a, Physicalreview letters , 078701 (2008).[4] D. Taylor, F. Klimm, H. A. Harrington, M. Kram´ar,K. Mischaikow, M. A. Porter, and P. J. Mucha, Naturecommunications , 1 (2015).[5] F. Papadopoulos, M. Kitsak, M. ´A. Serrano, M. Bogun´a,and D. Krioukov, Nature , 537 (2012).[6] A. Allard and M. ´A. Serrano, PLoS computational biol-ogy , e1007584 (2020).[7] M. Boguna, D. Krioukov, and K. C. Claffy, NaturePhysics , 74 (2009).[8] M. Xu, Q. Pan, A. Muscoloni, H. Xia, and C. V. Can-nistraci, Nature Communications , 1 (2020).[9] M. Bogun´a, F. Papadopoulos, and D. Krioukov, Naturecommunications , 1 (2010).[10] M. Boguna, I. Bonamassa, M. De Domenico, S. Havlin,D. Krioukov, and M. A. Serrano, arXiv:2001.03241 (2020), arXiv:2001.03241.[11] E. Estrada, Linear Algebra and its Applications ,4317 (2012).[12] D. Brockmann and D. Helbing, science , 1337 (2013).[13] M. De Domenico, Physical review letters , 168301(2017).[14] C. Hens, U. Harush, S. Haber, R. Cohen, and B. Barzel,Nature Physics , 403 (2019).[15] M. Kac, The American Mathematical Monthly , 369(1947).[16] J.-P. Bouchaud and A. Georges, Physics reports , 127(1990).[17] S. Havlin and D. Ben-Avraham, Advances in physics ,187 (2002).[18] J. D. Whitfield, C. A. Rodr´ıguez-Rosario, andA. Aspuru-Guzik, Physical Review A , 022323 (2010).[19] S. E. Venegas-Andraca, Quantum Information Process-ing , 1015 (2012).[20] V. Zaburdaev, S. Denisov, and J. Klafter, Reviews ofModern Physics , 483 (2015).[21] L. Giuggioli, Physical Review X , 021045 (2020).[22] N. Masuda, M. A. Porter, and R. Lambiotte, PhysicsReports , 1 (2017).[23] M. Rosvall and C. T. Bergstrom, Proceedings of the Na-tional Academy of Sciences , 1118 (2008).[24] J.-C. Delvenne, S. N. Yaliraki, and M. Barahona, Pro-ceedings of the national academy of sciences , 12755(2010).[25] M. Faccin, P. Migda(cid:32)l, T. H. Johnson, V. Bergholm, andJ. D. Biamonte, Physical Review X , 041012 (2014).[26] R. Lambiotte, J.-C. Delvenne, and M. Barahona, IEEETransactions on Network Science and Engineering , 76(2014).[27] M. De Domenico, A. Lancichinetti, A. Arenas, andM. Rosvall, Physical Review X , 011027 (2015).[28] Z. Burda, J. Duda, J.-M. Luck, and B. Waclaw, Physicalreview letters , 160602 (2009).[29] M. De Domenico, A. Sol´e-Ribalta, S. G´omez, andA. Arenas, Proceedings of the National Academy of Sci-ences , 8351 (2014).[30] S. Brin and L. Page, Computer Networks and ISDN Sys-tems , 107 (1998).[31] M. De Domenico and J. Biamonte, Physical Review X ,041062 (2016).[32] A. Ghavasieh and M. D. Domenico, Physical Review Re- search (2020), 10.1103/physrevresearch.2.013155.[33] R. R. Coifman and S. Lafon, Applied and computationalharmonic analysis , 5 (2006).[34] M. Cao, H. Zhang, J. Park, N. M. Daniels, M. E. Crov-ella, L. J. Cowen, and B. Hescott, PloS one (2013).[35] D. K. Hammond, Y. Gur, and C. R. Johnson, in (IEEE, 2013) pp. 419–422.[36] G. Bertagnolli, C. Agostinelli, and M. De Domenico,Journal of Complex Networks (2019), 10.1093/com-net/cnz041.[37] M. De Domenico, A. Sol´e-Ribalta, E. Cozzo, M. Kivel¨a,Y. Moreno, M. A. Porter, S. G´omez, and A. Arenas,Physical Review X , 041022 (2013).[38] M. Kivel¨a, A. Arenas, M. Barth´elemy, J. P. Gleeson,Y. Moreno, M. A. Porter, and E. Estrada, Journal ofComplex Networks , 203 (2014).[39] M. De Domenico, C. Granell, M. A. Porter, and A. Are-nas, Nature Physics , 901 (2016).[40] S. Gomez, A. Diaz-Guilera, J. Gomez-Gardenes, C. J.Perez-Vicente, Y. Moreno, and A. Arenas, Physical re-view letters , 028701 (2013).[41] S. Havlin, R. B. Selinger, M. Schwartz, H. Stanley, andA. Bunde, Physical Review Letters , 1438 (1988).[42] J. Gillis, in Mathematical Proceedings of the CambridgePhilosophical Society , Vol. 51 (Cambridge UniversityPress, 1955) pp. 639–651.[43] M. De Domenico, A. Sol´e-Ribalta, E. Omodei, S. G´omez,and A. Arenas, Nature communications , 1 (2015).[44] F. Battiston, V. Nicosia, and V. Latora, New Journal ofPhysics (2016), 10.1088/1367-2630/18/4/043035.[45] D. J. Klein and M. Randi´c, Journal of mathematicalchemistry , 81 (1993).[46] A.-L. Barab´asi and R. Albert, science , 509 (1999).[47] M. De Domenico, V. Nicosia, A. Arenas, and V. Latora,Nature communications , 1 (2015).[48] D. J. Watts and S. H. Strogatz, nature , 440 (1998).[49] M. Girvan and M. E. Newman, Proceedings of the na-tional academy of sciences , 7821 (2002).[50] N. Roberts and S. Everton, “Terrorist data: Noordin topterrorist network (subset).[machine-readable data file],”(2011).[51] F. Battiston, V. Nicosia, and V. Latora, Phys. Rev. E , 032804 (2014).[52] N. Mantel, Cancer research , 209 (1967).[53] P. Legendre and L. F. Legendre, Numerical ecology (El-sevier, 2012).[54] D. Krioukov, F. Papadopoulos, M. Kitsak, A. Vahdat,and M. Bogun´a, Physical Review E , 036106 (2010).[55] G. Bianconi and C. Rahmede, Scientific reports , 41974(2017).[56] K.-K. Kleineberg, M. Bogun´a, M. ´A. Serrano, and F. Pa-padopoulos, Nature Physics , 1076 (2016).[57] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lu-cas, A. Patania, J.-G. Young, and G. Petri, Physics Re-ports (2020).[58] J. R. Norris, Markov chains , 2 (Cambridge universitypress, 1998). Appendix A: Random walks and Markov chains
A Poisson process is a right-continuous process ( X t ) t ≥ with values { , , , . . . } with holding times S , S , . . . ( S i = J i − J i − is the time occurring between the ran-dom jump times J i − and J i ) that are independent ex-ponential random variables of rate 0 < λ < ∞ . In a gen-eralized Poisson process (or birth process) the parameter λ is allowed to depend on the current state of the pro-cess. Given its birth rates ≤ q i < ∞ for i = 0 , , , . . . ,( S i ) i ∈ N + are independent exponential random variableswith rates q i . Finally, a continuous-time Markov chain(MC) ( X t ) t ≥ on a finite set I with generator Q and ini-tial distribution p can be described in terms of a Poissonprocess. Each state i ∈ I of the process is a chamber anddoors close the passage to the other states. From time totime a single door opens (events cannot be simultaneous)allowing the process to change state and the doors openat the jump times of a Poisson process of rate q ij [58].The generator of the MC is indicated by Q , because it isa particular matrix, called Q − matrix in [58], satisfyingthree conditions(i) 0 ≤ − q ii < ∞ (ii) q ij ≥ ∀ i (cid:54) = j (iii) (cid:80) j q ij = 0 ∀ i .To recap, a continuous-time MC can be (equivalently)defined in terms of its jump chain and holding times, orof its transition probabilities given by the solution of theforward equation [58, Thm.2.8.2]. Appendix B: Maximal-entropy RW and thePerron–Frobenius theorem
An essential condition for the definition of themaximal-entropy random walk is that every componentof the eingenvector ψ , corresponding to the leading eigen-value λ max , be strictly positive. This is guaranteed by thePerron–Frobenius theorem for irreducible non-negativematrices, where a matrix A is said to be irreducibleif ∀ i, j = 1 , . . . , N there exists an integer m such that A mij >
0, which is exactly the irreducibility of a randomwalk on { , . . . , N }}