Universal multilayer network exploration by random walk with restart

ARTICLE Universal multilayer network exploration by random walk with restart Anthony Baptista 1,2✉, Aitor Gonzalez 2& Anaïs Baudot 1,3✉ The amount and variety of data have been increasing drastically for several years. These data are often represented as networks and explored with approaches arising from network theory. Recent years have witnessed the extension of network exploration approaches to capitalize on more complex and richer network frameworks. Random walks, for instance, have been extended to explore multilayer networks. However, current random walk approaches are limited in the combination and heterogeneity of networks they can handle. New analytical and numerical random walk methods are needed to cope with the increasing diversity and complexity of multilayer networks. We propose here MultiXrank, a method and associated Python package that enables Random Walk with Restart on any kind of multilayer network. We evaluate MultiXrank with leave-one-out cross-validation and link prediction, and measure the impact of the addition or removal of network data on prediction performances. Finally, we measure the sensitivity of MultiXrank to input parameters by in-depth exploration of the parameter space. https://doi.org/10.1038/s42005-022-00937-9 OPEN 1Aix-Marseille Univ, INSERM, MMG, Turing Center for Living Systems, CNRS, Marseille, France. 2Aix-Marseille Univ, INSERM, TAGC, Turing Center for Living Systems, Marseille, France. 3Barcelona Supercomputing Center, Barcelona, Spain. ✉email: [email protected];[email protected] COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys 1 1234567890():,; Data amount and variety have soared as never seen before, offering a unique opportunity to better understand complex systems. Among the different modes of representation of data, networks appear as particularly successful. Networks are indeed interesting to refine raw data and extract relevant features, patterns, and classes. They are exploited for years to study complex systems, and a wide and powerful range of tools from graph theory are available for their exploration. However, the integrated exploration of large multidimensional datasets remains a major challenge in many scientificfields. For instance, a comprehensive understanding of biological systems would require the integrated analysis of dozens of different datasets produced at different molecular, cellular or tissular scales. Recently, multilayer networks emerged as essential players in the analysis of such complex systems. Multilayer networks allow integrating more than one network in a unified formalism, in which the different networks are considered as layers1.For instance, Duran-Frigola et al.2combined 25 different networks of chemical compounds and their relationships, gathering relationships from chemical structures to clinical outcomes. This multilayer framework allows an integrated study of chemical compounds and their biological activities. Another example is given by the Hetionet project. The authors collected dozen of heterogeneous networks, i.e networks with various types of nodes such as genes, drugs or diseases, to prioritize drugs for repurposing3. Several definitions of multilayer networks have been proposed, based on the (in)homogeneity of the layers and the properties of the connections between layers4–6. For instance, multiplex networks are multilayer networks composed of different layers containing the same nodes (called replica nodes) but different types of edges, and thereby different topologies. Heterogeneous networks link networks composed of different types of nodes thanks to bipartite interactions. Temporal networks follow the dynamic of a network over time: all the layers have the same nodes, but each layer represents the interaction state at a given time7. We will here consider universal multilayer networks, which can be defined as multilayer networks composed of any number of multiplex (or monoplex) networks (with edges that can be directed and/or weighted), linked by bipartite networks (with edges that can be directed and/or weighted) (Fig. 1). A wide range of methods have been developed in the recent years to analyze multilayer networks. For instance, different network metrics have been adapted to multilayer networks8, as well as various network clustering algorithms for community detection9–11 or random walk for network exploration12–15. Random walks are iterative stochastic processes widely used to explore network topologies. They can be described as simulated particles that walk iteratively from one node to one of its neighbors with some probability16. The PageRank algorithm, for instance, is based on a random walk simulating the behavior of an internet user walking from one page to another thanks to hyperlinks. The user can also restart the walk on any arbitrary page17. In this particular random walk strategy, the restart prevents the random walker from being trapped in dead-ends18. An interesting alternative strategy restricts the restart to specific node(s), called the seed(s)19. In this strategy, named Random Walk with Restart (RWR) or Personalized PageRank, the random walk represents a measure of proximity from all the nodes in the network to the seed(s). RWR can also be described as a diffusion process, in which the objective is to determine the steady-state of an initial probability distribution20. RWR are widely used to exploit large-scale networks. In computational biology, for instance, RWR strategies have been shown to significantly outperform methods based on local distance measures for the prioritization of gene-disease associations21. Importantly, different upgrades of the RWR approach have been implemented during the last decade, including its extension to (i) heterogeneous networks12, (ii) multiplex networks13 and (iii) multiplex-heterogeneous networks15. In RWR, the degrees of freedom are summarized in the Transition rate matrix, and correspond to the available transitions between the different nodes of the graph. The extensions of RWR are challenging because the Transition rate matrices need to be normalized. To the best of our knowledge, this normalization is currently only solved for multilayer networks composed of two heterogeneous multiplex networks15,22 and the more universal case of Nmultiplex networks remains unsolved. We propose here MultiXrank, a framework composed of a method and a Python package to execute RWR on universal multilayer networks. We first introduce the mathematical bases of this RWR for universal multilayer networks, which correspond to a generalization of the approach from12. We evaluate MultiXrank with leave-one-out cross-validation and link prediction protocols. These evaluations reveal that more network data is not always better and highlight the critical influence of the bipartite networks. We finally present an in-depth exploration of the parameter space to measure the stability of the RWR output scores under variations of the input parameters. The MultiXrank Python package is freely available at https://github.com/anthbapt/ multixrank, with an optimized implementation allowing its application to large multilayer networks. Results Random walk with restart (RWR). Let us consider an irreducible and aperiodic Markov chain, for instance a network composed of a giant component with undirected edges, G=(V,E), where Vis the set of vertices and E⊆(V×V) is the set of edges. In the case of irreducible and aperiodic Markov chains, a stationary probability Fig. 1 A universal multilayer network. A universal multilayer network composed of three multiplex networks (green, blue and red multiplex networks). Each multiplex network contains different types of nodes (denoted 1 to 4, αto ϵ, and a to d, respectively). Their corresponding Supraadjacency matrices are denoted by Ai. The three multiplex networks are linked by six bipartite networks (represented here as bipartite interactions for the sake of visualization). The corresponding Bipartite network matrices are denoted by B i,j . It is to note that a connection between a node iin a first multiplex network and αand a node jin a second multiplex network β imposes the creation of edges between all replicas of node ipresent in the different layers of the multiplex network αand all replicas of node jpresent in the different layers of multiplex network β. All the edges of the universal multilayer networks can be weighted and/or directed. ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 2COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys p*exists and satisfies the following properties: pðiÞ>0;8i2V ∑i2VpðiÞ¼1 ð1Þ We next introduce the probability defining the walk from one node to another. Let us define x, a particle that explores the network, x t its position at time tand x t+1 its position at time t+1. Considering two nodes iand j: Pðxtþ1¼jjxt¼iÞ¼ 1 diif ði;jÞ2E 0 Otherwise (ð2Þ with d i being the degree of the node i. All the normalized possible transitions can be included in the Transition rate matrix. This Transition rate matrix, noted M, can be seen as the matrix of the degrees of freedom of the particle in the system. It is useful to note that the Transition rate matrix is equal to the columnnormalized Adjacency matrix. The distribution denoted by pt¼ ðptðiÞÞi2Vdescribes the probability of being in the node iat time t, and the stationary distribution p*is obtained thanks to the homogeneous linear difference equation [3]18,23: pT tþ1¼MpT tð3Þ with pT tdenoting the transpose of the vector p t . Moreover, we can introduce a non-homogeneous linear difference equation [4]23 to take into account the restart on the seed(s). When the Transition rate matrix is a Stochastic matrix, the stationary distribution is reached18 (Supplementary Note 1.A.1 for elements of proof of convergence) and this distribution can be seen as a measure of proximity of all the network nodes with respect to the seed(s). pT tþ1¼ð1rÞMpT tþrpT 0ð4Þ The distribution p 0 corresponds to the initial probability distribution, where only the seed(s) have non-zero values; r represents the restart probability. RWR on multiplex networks. The RWR method has been extended to multiplex networks, i.e., multilayer networks with a one-to-one mapping between the (replica) nodes of the different layers (Fig. 1)1,13,14. Multiplex networks can be represented by Supra-adjacency matrices, which correspond to a generalization of the standard Adjacency matrix. In the following, we will use several multiplex networks, indexed by k. We denoted by Akthe Supra-adjacency matrix of the multiplex network indexed by k. The Adjacency matrix of the layer lof the multiplex network kis denoted by A½l k. The element of this adjacency matrix from node i to node jis defined as ðA½l kÞi;j≥0. The dimension of the Supraadjacency matrix Akof the multiplex network kis equal to (L k *n k )*(L k *n k ), with n k the number of nodes in each layer of the multiplex network kand L k the number of layers in the multiplex network k. The Supra-adjacency matrix Akis defined as follows: ðAkÞil;jm¼A½l k  i;jif l¼m δi;jif l≠m 8 < :ð5Þ where δdefines the Kronecker delta (i.e., 1 if iequal jand 0 otherwise), and land mrepresent the layers of the multiplex network k. We can also define a multiplex network as a set of nodes, VAkand a set of edges, EAk: GAk¼ðVAk;EAkÞ VAk¼fvl i;i¼1;¼;nk;l¼1;¼;Lkg EAk¼fell i;j;i;j¼1;¼;nk;l¼1;¼;Lk;ðA½l kÞi;j≠0g ∪felm i;i;i¼1;¼;nk;l≠mg 8 > > > > > < > > > > > : ð6Þ Importantly, we need to column-normalize the Supra-adjacency matrix defined in the equations [5–6] in order to converge to the steady-state, as defined in15. This normalization requires including the parameters δ k related to the jumps from one layer to another inside the matrix representation, as described in13 (Fig. 2). In the next section, we need to index by kall the parameters that are dedicated to the multiplex network k. The Supra-adjacency matrix representing the multiplex network kcan be written as described in equation [7]. The matrix I k represents the Identity matrix of size n k . Ak¼ ð1δkÞA½1 k δk ðLk1ÞIk¼ δk ðLk1ÞIk δk ðLk1ÞIkð1δkÞA½2 k¼ δk ðLk1ÞIk . . .. . ... .. . . δk ðLk1ÞIk δk ðLk1ÞIk¼ð1δkÞA½Lk k 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 ð7Þ Fig. 2 MultiXrank Random Walk with Restart parameters. Parameters of the Random Walk with Restart allowing to explore universal multilayer networks composed of Nmultiplex networks (each composed of several layers containing the same set of (replica) nodes but different edges). The parameters δare associated with the probability to jump from one layer to another in a given multiplex network, λwith the probability to jump from one multiplex network to another multiplex network, τwith the probability to restart in a given layer of a given multiplex network, and ηwith the probability to restart in a given multiplex network. COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 ARTICLE COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys 3 RWR on universal multilayer networks.Weheredefine a RWR method that can be applied to universal multilayer networks. Universal multilayer networks are composed of any combination of multiplex networks, linked by any combination of bipartite networks (Fig. 1). All network edges can also be weighted and/or directed. The formalism for the application of RWR on multiplex networks is described in the previous section. We will now detail the Bipartite network matrices, and how to combine intra- and intermultiplex networksinformationtoobtaintheSupra-heterogeneous adjacency matrix. The Supra-heterogeneous adjacency matrix will embed all the possible transitions in a universal multilayer network. Bipartite networks connect heterogeneous nodes. The Bipartite network matrices contain the transitions between different types of nodes present in different networks. If the network αhas n α nodes, and the network βhas n β nodes, the Bipartite network matrix denoted b α,β has a size equal to n α *n β . Now, let us define Aαand Aβ, two Supra-adjacency matrices representing the multiplex networks αand β. The Bipartite network matrix B α,β represents the transitions from the nodes of the multiplex network αto the nodes of the multiplex network β. The size of the Bipartite network matrix B α,β is equal to (L α *n α )*(L β *n β ). The Bipartite network matrices are composed of (L α *L β ) times the Bipartite network matrix b α,β (equation [8]). The matrix b α,β is composed of all the transitions from one layer of the multiplex network αto one layer of the multiplex network β. We extended the formalism used in15 in order to consider more than two different multiplex networks. Bα;β¼ bα;βbα;β¼bα;β bα;βbα;β¼bα;β . . .. . ... .. . . bα;βbα;β¼bα;β 2 6 6 6 6 6 4 3 7 7 7 7 7 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Lβtimes 9 > > > > > > > > = > > > > > > > > ; Lαtimes ð8Þ The representation of the bipartite networks as a set of nodes VBand a set of edges EBcan be written as: GB¼ðVB;EBÞ VB¼fvα k;k¼1;¼;nαg∪fvβ l;l¼1;¼;nβg EB¼feαβ k;lk¼1;¼;nα;l¼1;¼;nβ;ðbα;βÞk;l≠0g 8 > > < > > : ð9Þ It is to note that if the bipartite networks are undirected, bT β;α¼ bα;βand BT β;α¼Bα;β. Universal multilayer networks unify the representation of heterogeneous multiplex networks. We previously defined the Supraadjacency matrices of each multiplex network and the Bipartite network matrices connecting the different multiplex networks. We now introduce the Supra-heterogeneous adjacency matrix, denoted by S.Thismatrix,defined in equation [10], collects the NSupra-adjacency matrices representing each multiplex network, A1;A2;¼;AN,andtheN*(N−1) Bipartite network matrices connecting each multiplex network, B 1,2 , B 1,3 ,…,B 1,N ,B 2,1 ,…,B N,N−1 . S¼ A1B1;2¼B1;N B2;1A2¼B2;N . . .. . ... .. . . BN;1BN;2¼AN 2 6 6 6 6 6 4 3 7 7 7 7 7 5 ð10Þ We can also define the Supra-heterogeneous adjacency matrix as a set of nodes and edges: GS¼VS;ES  VS¼S N k¼1 fvαk k;i;i¼1;¼;nk;αk¼1;¼;Lkg ES¼S N k¼1 feαk;αk i;j;i;j¼1;¼;nk;A½αk k  i;j ≠0g  ∪feαk;βk i;i;i¼1;¼;nk;αk≠βk;αk;βk¼1;¼;Lkg ∪S N k;l¼1;k≠l feαk;αl i;j;i¼1;¼;nk;j¼1;¼;nl;Bk;l  i;j ≠0g 8 > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > :ð11Þ The normalization of the Supra-heterogeneous adjacency matrix ensures the convergence of the RWR to the steady-state.Themost complex issue is the normalization of the Supra-heterogeneous adjacency matrix into a Transition rate matrix that can be used in equation [4]. The normalization allows obtaining a Stochastic matrix that guarantees the convergence of the RWR to the steady-state18 (see elements of proof in Supplementary Note 1.A.1). It is important to note that we have chosen a column normalization. The resulting normalized matrix, denoted by b Sis definedinequation[12]. We generalized the formalism of Li and Patra12 established for two heterogeneous monoplex networks (Supplementary Note 1.D). This generalization to universal multilayer networks is done thanks to the intra- and intermultiplex network normalizations defined in equations [13–14], with α∈[[1, N]], β∈[[1, N]]. In addition, ciαis the number of bipartite networks in which the node i α appears as source of the multiplex network αdenoted by M α . b S¼bS11 bS12 ¼bS1N bS21 bS22 ¼bS2N . . .. . ... .. . . bSN2bSN2¼bSNN 2 6 6 6 6 6 4 3 7 7 7 7 7 5 ð12Þ In equation [13], bSαα defines the transition probabilities inside a given multiplex network. In the case of a multiplex network, if a node has no bipartite interactions with nodes from another multiplex networks, we can use the standard normalization. If bipartite interactions exist, then the normalization takes into account the probability that the walker can stay in the multiplex network ð1∑ciα β¼1λαβÞ.Inequation[14], bSαβ defines the transition probability between two different multiplex networks. There are here three possibilities. If the node has no bipartite interactions, the transition probability is equal to zero. If the node has bipartite interactions, the transition probability is equal to the standard normalization weighted by the jump probability (λ αβ ). Finally, if the node exists only in the bipartite network, the normalization corresponds to the standard normalization weighted by a modified jump probability. This normalization takes into account all the bipartite interactions of the considered node. bSααðiα;jαÞ¼ Aαðiα;jαÞ ∑ nα kα¼1 Aαðiα;kαÞ if 8β:∑ nβ kβ¼1 Bα;βðiα;kβÞ¼0 1∑ ciα β¼1 λαβ  Aαðiα;jαÞ ∑ nα kα¼1 Aαðiα;kαÞ Otherwise 8 > > > > > > > < > > > > > > > : ð13Þ ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 4COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys bSαβðiα;jβÞ¼ λαβBα;βðiα;jβÞ ∑ nβ kβ¼1 Bα;βðiα;kβÞ if ∑ nβ kβ¼1 Bα;βðiα;kβÞ≠0 λαβ ∑ c β¼1 λαβ ∑ c iα¼1Bα;βðiα;jβÞ ∑ c iα¼1 ∑ nβ kβ¼1 Bα;βðiα;kβÞ if iαnot in Mα 0 Otherwise 8 > > > > > > > > > > < > > > > > > > > > > : ð14Þ The normalization allows including the parameters λ αβ to jump between the multiplex networks (Fig. 2). In other words, these parameters weight the jumps from one multiplex network αto another multiplex network β, if the bipartite interaction exists. Moreover, the standard probability condition of normalization imposes that ∑N α¼1λαβ ¼1;8β, where Nrepresents the number of multiplex networks. Finally, the RWR equation on universal multilayer networks is defined as: pT tþ1¼ð1rÞbSpT tþrpT 0:ð15Þ RWR initial probability distribution in universal multilayer networks. The initial probability distribution p 0 from equation [15], which contains the probabilities to restart on the seed(s), can be written in its general form as follows: pT 0¼ η1 v1 0 η2 v2 0 ¼ ηN vN 0 2 6 6 6 43 7 7 7 5ð16Þ where η k is the probability to restart in one of the layers of the multiplex network k,and vk 0is the initial probability distribution of the multiplex network k.Thesizeof vk 0is equal to (L k *n k ), where L k isthenumberoflayersinthemultiplexnetworkkand n k is the number of nodes in the multiplex network k. We constraint the parameter ηwith the standard condition of normalization of the probability that imposes ∑N k¼1ηk¼1. We defined another parameter, τ, to take into account the probability of restarting in the different layers of a given multiplex network. This parameter includes τ kj ,wherekcorresponds to the index of the multiplex network, and jto the index of the layer of the multiplex network k (Fig. 2). In other words, τ kj corresponds to the probability to restart in the jth layer of the multiplex network k.Finally, vk 0is defined as follows:  vk 0¼½τk1vk 0;τk2vk 0;¼;τkLkvk 0T,withvk 0being a vector with 1/ω k in the position(s) of seed(s) and zeros elsewhere, and ω k being the number of seeds in the multiplex network k. The standard condition of normalization of the probability gives the constraint: ∑Lk j¼1τkj ¼1, ∀k. Numerical implementation: multiXrank. Our RWR on universal multilayer networks is implemented as a Python package called MultiXrank (Supplementary Note 2). MultiXrank has an optimized implementation. Default parameters allow exploring homogeneously the multilayer network (Supplementary Note 1.B). The running time of the package depends on the number of edges of the multilayer network (complexity analyses in Supplementary Note 2.A). The package is available on GitHub https://github.com/anthbapt/multixrank, and can be installed with standard pip installation command: https://pypi.org/project/ MultiXrank. Evaluations. We evaluated the performances of MultiXrank using two different multilayer networks. The first one is a large biological multilayer network composed of two multiplex networks and one monoplex network. It contains a gene multiplex network gathering gene physical and functional relationships, a drug multiplex network containing drug clinical and chemical relationships, and a disease monoplex network representing disease phenotypic similarities. Each monoplex/multiplex network is connected to the others thanks to bipartite networks containing gene-disease, drug-gene, and drug-disease interactions (Supplementary Note 3.B). The second multilayer network is composed of three multiplex networks. It contains a French airports multiplex network, a British airports multiplex network, and a German airports multiplex network. In each multiplex network, the nodes represent the airports of each country and the edges represent the national flight connections between these airports for three different airline companies. The three multiplex networks are linked with bipartite networks corresponding to transnational flight connections (Supplementary Note 3.A). We designed a Leave-One-Out Cross-Validation (LOOCV) protocol inspired by F.Mordelet and J.P.Vert24 and A.Valdeolivas et al.15. In this protocol, we systematically leave-out some known associations and assess the reconstruction of this left-out data using the data remaining in the network (Supplementary Note 4.A and Fig. S9). In the case of the biological multilayer network,we systematically left-out known gene-disease associations. More specifically, for each disease associated with at least two genes, each gene is remove one-by-one and considered as the left-out gene. The remaining gene(s) associated with the same disease are used as seed(s). When the disease network is considered in the evaluation, the disease node is used as seed together with the gene node(s). The RWR algorithm is then applied, and all the network nodes are scored according to their proximity to the seed(s). The rank of the gene node that was left-out in the ongoing run is recorded. The perfect ranking for the left-out gene is 1; the closer the rank is to 1, the better the prediction. The gene left-out process is repeated iteratively for all the genes. Finally, the Cumulative Distribution Function (CDF) of the ranks of the leftout genes is plotted (Fig. 3). The CDF displays the ratio of left-out genes that are ranked by the RWR within the top-Kranked gene nodes. The CDFs are used to evaluate and compare the performance of the RWR applied to different combinations of biological networks: the protein-protein interactions (PPI) network alone, the gene multiplex network, the multilayer network composed of the gene multiplex and the disease monoplex networks, and the multilayer network composed of the gene and drug multiplex networks and the disease monoplex network (Fig. 3a). We observed that considering multiple sources of network data is always better than considering the PPI alone. In addition, considering multilayer information is better than considering only the gene multiplex network. However, the increased performances in the LOOCV seem to arise only from combining the gene multiplex network with the disease monoplex network (and associated gene-disease bipartite network). Indeed, the addition of the drug multiplex network (and associated drug-gene and drug-disease bipartite networks) to the multilayer system does not increase the performances (Fig. 3a). We repeated the same LOOCV protocol for the airports multilayer network, in which the left-out nodes are French airport nodes associated with a given British airport node. Here, the behavior is different, as adding the third multiplex network containing German airports connections (and associated French- German and British-German bipartite networks) increases the performances of the RWR to predict the associations between French and British airports (Fig. 3b). To better understand these different behaviors, we examined in detail the amount of common nodes (called overlaps) existing between the nodes of the different bipartite networks. We COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 ARTICLE COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys 5 observed that only 23% of the genes from the gene-disease bipartite network are present in the drug-gene bipartite network. Similarly, only 5% of the diseases from the gene-disease bipartite network are present in the disease-drug bipartite network (Fig. S10). Given these low overlaps, the drug multiplex network might not contribute significantly to connecting gene and disease nodes during the random walks. This might explain why adding the drug multiplex network does not improve the performances of the LOOCV. Contrarily, the bipartite networks of the airport multilayer network displays high overlaps (Fig. S10). These high overlaps might explain why the addition of the third multiplex network in this case increases the predictive power (Fig. 3b). To validate the proposed central role of bipartite networks in the RWR performances, we artificially increased the connectivity of the gene-drug and disease-drug bipartite networks before applying the same LOOCV protocol. To this goal, we added artificial transit drug nodes linking existing gene-disease associations (strategy described in Supplementary Note 4C and Fig. S12). We observed that these artificially added transit nodes increased drastically the performances of the LOOCV (Fig. 3c). The same phenomenon is observed for the airports multilayer network (Fig. 3d). In addition, we checked if random perturbations in these artificially enhanced bipartite networks would decrease the performances of the LOOCV. To do so, we progressively randomized the edges in the bipartite networks with artificially increased connectivity, until obtaining completely random bipartite networks. We observed that the progressive randomization of the bipartite networks continuously decreases the Fig. 3 Evaluation and comparison of multiXrank performances on different combinations of multilayer networks. a,bCumulative Distribution Functions (CDFs) representing the ranks of the left-out nodes in the Leave-One-Out Cross-Validation (LOOCV) protocol. a: focus on different combinations of biological networks: protein-protein interactions network alone (PPI), gene multiplex network (multi-1), multilayer network composed of the gene multiplex network and the disease monoplex network (multi-2), and multilayer network composed of the gene and drug multiplex networks and the disease monoplex network, for two different sets of parameters (multi-3, multi-3 bis). The multilayer networks are connected by the bipartite networks described in the Evaluations section. b: focus on different combinations of airports networks: French multiplex network (multi-1), multilayer network composed of the French and British airports multiplex networks (multi-2), and multilayer network composed of the French, British, and German airports multiplex networks, for two different sets of parameters (multi-3, multi-3 bis). These multilayer networks are connected by the bipartite networks described in the Evaluations section. c,dCDFs representing the ranks of the left-out nodes in the LOOCV protocol for the multi-3 multilayer networks described previously with artificially increased connectivity in the gene-drug and disease-drug bipartite networks. cThe connectivity is artificially increased thanks to the addition of 1 (multi3+1), 2 (multi3+2) or 5 (multi3+5) transit drug nodes for each gene-disease association. d: In the airport multilayer network, the connectivity is artificially increased in the French-German and British-German bipartite networks thanks to the addition of 1 (multi3+1), 2 (multi3+2) or 5 (multi3+5) transit German nodes for each French-British airports association. The parameters of the Random Walk with Restart (RWR) are detailed in Supplementary Tables S5–S6. ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 6COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys predictive power of the RWR up to obtaining the same performances as with only two multiplex networks (Fig. S13.A for the airport multilayer networks and S13.B for the biological multilayer networks). Finally, we repeated all these evaluations using a standard Link Prediction (LP) protocol (Supplementary Note 4.B). LP has already been used to measure the predictive power of RWR methods25.In the LP protocol, we systematically removed gene-disease edges from the gene-disease bipartite network, and predicted the rank of the removed gene using the disease as seed in the RWR. The LP protocol is applied on the airport multilayer network by removing a French-British edge from the French-British bipartite network, and predicting the rank of the French airport using the British airport node as a seed in the RWR. We overall observed similar behaviors as in the LOOCV (Fig. S11 and S14). Importantly, the LOOCV and LP protocols can be used to evaluate the pertinence of adding new multiplex networks in a multilayer network or new network layers in a multiplex network. Both evaluation protocols are available within the MultiXrank package. Parameter space exploration. We next evaluated the stability of MultiXrank output scores upon variations of the input parameters. We illustrate this exploration of the parameter space with the biological multilayer network composed of the gene multiplex network and the disease monoplex network. We first compared the top-5 and top-100 gene and disease nodes prioritized by MultiXrank using 125 different sets of parameters (see Supplementary Note 5 for the definition of the sets of parameters). We observed that the top-ranked gene nodes vary more depending on the input parameters than the top-ranked disease nodes (Fig. 4a). To better understand the stability of the output scores upon variations of the input parameters, we proposed a protocol based on 5 successive steps: (i) definition of the sets of parameters, (ii) construction of a matrix containing the similarities of the RWR output scores obtained with each set of input parameters, using a the similarity measure defined in equation [17]. The similarities are computed for each type of node independently (i.e., for gene and disease nodes independently). Θk γσ ¼∑ nk j¼1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 rk γ  jrk γσ  j hi 0 @1 A 2 þ1 rk σ ðÞ jrk σγ  j hi 0 @1 A 2 v u u u trk γ  jþrk σ ðÞ j 2 ! 2 ð17Þ where γand σdefine two sets of parameters, n k is the number of nodes associated with the multiplex network k. In addition, rk γ (resp. rk σ) is the rank output scores distribution that associates with each node its rank given by the RWR with the set of parameters γ(resp. σ) for the multiplex network k. Finally, rk γσ (resp. rk σγ) gives to each node of the output scores distribution obtained by the set of parameters γ(resp. σ) (in the multiplex network k) their rank in the distribution σ(resp. γ). We next computed a consensus Similarity matrix with a normalized euclidean norm of each individual Similarity matrix (equation [18]). Θγσ ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∑ N k¼1 ðΘk γσÞ2 nk v u u tð18Þ where Nis the number of multiplex networks. The next step is (iii) projection of the consensus Similarity matrix into a Principal Component Analysis (PCA) space (Fig. 4b). In this PCA space, each dot represents the output scores resulting from a set of parameters. Then, (iv) clustering (using k-means on the two first principal components) to identify sub-regions containing similar RWR output scores. Finally, (v) comparing the top-ranked nodes obtained with the set of parameters belonging to each cluster (Fig 4c, Supplementary Note 5). We applied this protocol to evaluate the output scores obtained by MultiXrank on the previously defined biological multilayer network composed of the gene multiplex network and the disease monoplex network, using 125 different combinations of parameters (Fig. 4, supplementary Fig. S16). We projected the consensus Similarity matrix into a PCA space and identified 8 clusters (Fig. 4b). To illustrate the behavior inside clusters, we concentrated our analyses on the two clusters defined in the bottom left subspace (clusters number 4 and 6, zoom-in Fig. 4b). The top-100 ranked gene and disease nodes inside each of the two clusters are overall similar (Fig. 4c). This means that, even if the node prioritization can be sensitive to input parameters, we can identify regions of stability in the parameter space. Moreover, the protocol allows identifying the monoplex/multiplex networks that generate most variability in the output scores upon changes in the input parameters. We applied the parameter space exploration protocol to other multilayer networks and observed diverse behaviors, from highly variable top-rankings and scattered projections in the PCA space for the airport multilayer network (Supplementary Fig. S15) to robust top-rankings with well-clustered projections in the PCA space for the biological multilayer network composed of 3 types of nodes (genes, diseases and drugs, Supplementary Fig. S16). Overall, our parameter space study reveals different sensitivities to input parameters depending on the multilayer network explored. The protocol is available within the MultiXrank package and can be used to characterize in-depth the sensitivity to input parameters of any multilayer network. Discussion Multilayer networks are nowadays very popular, in particular because they allow capturing a larger part of real and engineered systems. In biology, multilayer networks integrating multiscale sources of heterogeneous interactions provide a more comprehensive picture of biological system functionalities. However, data representation as multilayer networks must be accompanied by the development of tools allowing their exploration. Many efforts are thereby dedicated to extend classical network theory algorithms to multilayer systems5,26. These algorithms include for instance clustering algorithms27, Graph Convolutional Networks28,29 or meta-path based methods3,30. Other important network exploration algorithms, such as diffusion kernels or methods based on random walk, are based on the principle of network propagation26.Themethods based on random walk, such as PageRank, biased random walk or Random Walk with Restart (RWR), are widely used in network science. They are indeed versatile: the random walk output scores can be used directly for node prioritization and subnetwork extraction, but can also be used as input for downstream analyses, for instance for supervised classification or node embedding22. Different random walk methods have been adapted to consider multilayer networks. However, a large variety of multilayer networks exist, from multiplex to temporal networks, for instance. To the best of our knowledge, network exploration algorithms that have been adapted to handle multilayer networks can usually be applied only to specificcategoriesof COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 ARTICLE COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys 7 multilayer networks, such as multiplex networks composed of thesamesetofnodes. We present here MultiXrank, a tool that proposes an optimized and general formalism for RWR on universal multilayer networks. MultiXrank can be applied to explore multilayer networks composed of any combination of multiplex, monoplex or bipartite networks, and all the network edges can be directed and/or weighted. To the best of our knowledge, any type of multilayer networks could be represented with our formalism, even if it might sometimes require some adaptations. We illustrated the use of MultiXrank with RWR on biological and airport multilayer networks and thereby provide guidelines for users. Even if one’s initial intuition in data analysis could be that “more data is better”, the addition of interaction network layers also brings additional degrees of freedom5. To evaluate the pertinence of the addition of multiplex networks or the addition of layers in a multilayer system, MultiXrank includes a systematic evaluation protocol based on Leave-One-Out-Cross-Validation and Link Prediction. Overall, our results show that adding networks data does not always increase the predictive power of the RWR, as already suggested by previous studies11. Our evaluation protocol can be used, for the first time to our knowledge, to evaluate indepth the signal-to-noise of multilayer system combinations. Finally, we complemented MultiXrank with a parameter space exploration protocol to measure the influence of varying the input parameters on the global stability of the output scores. It is to note that this parameter space exploration protocol is universal and can be used to study any complex system exploration approach providing scores as outputs. The output scores of MultiXrank can be used in a wide variety of downstream analyses. For instance, shallow embedding methods need similarity measures for the optimization of the loss function22,31. MultiXrank can produce such a similarity measure respecting the global topology of the multilayer network. An interesting application could be to use MultiXrank output scores for embedding and evaluate the predictive power of the gene-disease association prediction task. Indeed, the embedding is expected to be more robust to the noise than the direct network space32. The MultiXrank package can be applied to any kind of multilayer network such as social, economic, or ecological multilayer networks. MultiXrank is optimized and can handle multilayer networks containing up to millions edges. To consider billion-scale network problems, several strategies could be considered, such as the Block Elimination Approach for RWR (BEAR) that can be exact or approximate33 or the Best of Preprocessing and Iterative approaches (BEPI) that is an approximate approach34. Data availability All the data and the code used in the article are available on an OSF repository: https:// osf.io/zsmua (DOI 10.17605/OSF.IO/ZSMUA). This repository includes all the results obtained in the article. Code availability The package is available on GitHub https://github.com/anthbapt/multixrank, can be installed with standard pip installation command: https://pypi.org/project/MultiXrank, and is associated with complete documentation: https://multixrank-doc.readthedocs.io/ en/latest. Fig. 4 Exploration of multiXrank parameter space. a Comparison of the top-5 and top-100 nodes ranked by MultiXrank using a biological multilayer networks composed of the gene multiplex network and the disease monoplex network for 125 different sets of parameters. The top-5 or top-100 ranked nodes for each set of parameters are merged, and the number of occurrences of each node are counted. The nodes are represented in bars colored in red when the node is found in all top-5 or top-100 scores, and in blue otherwise. bClustering in the Principal Component Analysis (PCA) space of the output scores obtained with MultiXrank on the biological multilayer network composed of the gene multiplex network and the disease monoplex network using 125 different sets of parameters. The zoom-in emphasizes the clusters number 4 and 6. cComparison of the top-100 nodes retrieved for the sets of parameters belonging to clusters 4 and 6 defined in (b). The bar is colored in red when a node is found in all top-100 scores, and in blue otherwise. The parameters of the Random Walk with Restart (RWR) are detailed in Supplementary Table S7. ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 8COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys Received: 13 September 2021; Accepted: 8 June 2022; References 1. Bianconi, G. Multilayer Networks: Structure and Function. (Oxford University Press, Oxford, 2018). 2. Duran-Frigola, M. et al. Extending the small-molecule similarity principle to all levels of biology with the chemical checker. Nat. Biotechnol. 38, 1087–1096 (2020). 3. Himmelstein, D. S. et al. Systematic integration of biomedical knowledge prioritizes drugs for repurposing. eLife 6, e26726 (2017). 4. De Domenico, M. et al. Mathematical formulation of multilayer networks. Phys. Rev. X 3, 041022 (2013). 5. Kivelä, M. et al. Multilayer networks. J. Complex Netw. 2, 203–271 (2014). 6. Lee, B., Zhang, S., Poleksic, A. & Xie, L. Heterogeneous multi-layered network model for omics data integration and analysis. Front. Genet. 10, 1381 (2020). 7. Holme, P. & Saramäki, J. Temporal networks. Phys. Rep. 519,97–125 (2012). 8. Battiston, F., Nicosia, V. & Latora, V. Structural measures for multiplex networks. Phys. Rev. E 89, 032804 (2014). 9. Mucha, P. J., Richardson, T., Macon, K., Porter, M. A. & Onnela, J.-P. Community structure in time-dependent, multiscale, and multiplex networks. Science 328, 876–878 (2010). 10. Didier, G., Brun, C., Baudot, A. & Gomez, S. Identifying communities from multiplex biological networks. PeerJ 3, e1525 (2015). 11. Choobdar, S. et al. Assessment of network module identification across complex diseases. Nat. Methods 16, 843–852 (2019). 12. Li, Y. & Patra, J. C. Genome-wide inferring gene-phenotype relationship by walking on the heterogeneous network. Bioinformatics 26, 1219–1224 (2010). 13. De Domenico, M., Solé-Ribalta, A., Gómez, S. & Arenas, A. Navigability of interconnected networks under random failures. Proc. Natl Acad. Sci. 111, 8351–8356 (2014). 14. Cho, H., Berger, B. & Peng, J. Compact integration of multi-network topology for functional analysis of genes. Cell Syst. 3, 540–548.e5 (2016). 15. Valdeolivas, A. et al. Random walk with restart on multiplex and heterogeneous biological networks. Bioinformatics 35, 497–505 (2018). 16. Lovász, L. Random walks on graphs: a survey. Combinatorics, Paul. Erdos is. Eighty 2,1–46 (1993). 17. Brin, S. & Page, L. The anatomy of a large-scale hypertextual web search engine. Computer Netw. ISDN Syst. 30, 107–117 (1998). Proceedings of the Seventh International World Wide Web Conference. 18. Langville, A. N. & Meyer, C. D. Google’s PageRank and Beyond: The Science of Search Engine Rankings. (Princeton University Press, USA, 2006). 19. Pan, J.-Y., Yang, H.-J., Faloutsos, C. & Duygulu, P. Automatic multimedia cross-modal correlation discovery. In Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’04, 653–658 (Association for Computing Machinery, New York, NY, USA, 2004). https://doi.org/10.1145/1014052.1014135. 20. Gómez, S. et al. Diffusion dynamics on multiplex networks. Phys. Rev. Lett. 110, 028701 (2013). 21. Köhler, S., Bauer, S., Horn, D. & Robinson, P. N. Walking the interactome for prioritization of candidate disease genes. Am. J. Hum. Genet. 82, 949–958 (2008). 22. Pio-Lopez, L., Valdeolivas, A., Tichit, L., Remy, E. & Baudot, A. Multiverse: a multiplex and multiplex-heterogeneous network embedding approach. Sci. Rep. 11, 8794 (2021). 23. Meyer, C. D. Matrix Analysis and Applied Linear Algebra. (Society for Industrial and Applied Mathematics, USA, 2000). 24. Mordelet, F. & Vert, J.-P. Prodige: prioritization of disease genes with multitask machine learning from positive and unlabeled examples. BMC Bioinforma. 12, 389 (2011). 25. Zhou, M., Zheng, C. & Xu, R. Combining phenome-driven drug-target interaction prediction with patients’electronic health records-based clinical corroboration toward drug discovery. Bioinformatics 36, i436–i444 (2020). 26. Boccaletti, S. et al. The structure and dynamics of multilayer networks. Phys. Rep. 544,1–122 (2014). 27. Huang, X., Chen, D., Ren, T. & Wang, D. A survey of community detection methods in multilayer networks. Data Min. Knowl. Discov. 35,1–45 (2021). 28. Ghorbani, M., Baghshah, M. S. & Rabiee, H. R. Mgcn: Semi-supervised classification in multi-layer graphs with graph convolutional networks. In Proceedings of the 2019 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining, ASONAM ’19, 208-211 (Association for Computing Machinery, New York, NY, USA, 2019). https://doi.org/10. 1145/3341161.3342942. 29. Shanthamallu, U. S., Thiagarajan, J. J., Song, H. & Spanias, A. Gramme: Semisupervised learning using multilayered graph attention models. IEEE Trans. Neural Netw. Learn. Syst. 31, 3977–3988 (2020). 30. Zhang, X., Zou, Q., Rodríguez-Patón, A. & ZENG, X. Meta-path methods for prioritizing candidate disease mirnas. IEEE/ACM Trans. Computational Biol. Bioinforma. 16, 283–291 (2019). 31. Hamilton, L., Ying, W., R. & Leskovec, J. Representation learning on graphs: Methods and applications (v3). https://arxiv.org/abs/1709.05584 (2018). 32. Nelson, W. et al. To embed or not: Network embedding as a paradigm in computational biology. Front. Genet. 10, 381–381 (2019). 33. Shin, K., Jung, J., Lee, S. & Kang, U. Bear: Block elimination approach for random walk with restart on large graphs. In Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data, SIGMOD ’15, 1571-1585 (Association for Computing Machinery, New York, NY, USA, 2015). https://doi.org/10.1145/2723372.2723716. 34. Jung, J., Park, N., Lee, S. & Kang, U. Bepi: Fast and memory-efficient method for billion-scale random walk with restart. In Proceedings of the 2017 ACM International Conference on Management of Data, SIGMOD ’17, 789-804 (Association for Computing Machinery, New York, NY, USA, 2017). https://doi.org/10.1145/3035918.3035950. Acknowledgements The project leading to this preprint has received funding from the ≪Investissements d’Avenir ≫French Government program managed by the French National Research Agency (ANR-16-CONV-0001), from Excellence Initiative of Aix-Marseille University - A*MIDEX and from the Inserm Cross-Cutting Project GOLD. Author contributions A.Bap. and A.Bau. designed research; A.Bap. performed research; A.Bap. analyzed data; A.Bap. and A.G. contributed to packaged code; A.Bap. and A.Bau. wrote the paper. Competing interests The authors declare no competing interests. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s42005-022-00937-9. Correspondence and requests for materials should be addressed to Anthony Baptista or Anaïs. Baudot. Peer review information Communications Physics thanks Albert Solé-Ribalta, Joao Gama Oliveira and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available. Reprints and permission information is available at http://www.nature.com/reprints Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2022 COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00937-9 ARTICLE COMMUNICATIONS PHYSICS | (2022) 5:170 | https://doi.org/10.1038/s42005-022-00937-9 | www.nature.com/commsphys 9