MultiVERSE: a multiplex and multiplex-heterogeneous network embedding approach

Network embedding approaches are gaining momentum to analyse a large variety of networks. Indeed, these approaches have demonstrated their effectiveness in tasks such as community detection, node classification, and link prediction. However, very few network embedding methods have been specifically designed to handle multiplex networks, i.e. networks composed of different layers sharing the same set of nodes but having different types of edges. Moreover, to our knowledge, existing approaches cannot embed multiple nodes from multiplex-heterogeneous networks, i.e. networks composed of several multiplex networks containing both different types of nodes and edges. In this study, we propose MultiVERSE, an extension of the VERSE framework using Random Walks with Restart on Multiplex (RWR-M) and Multiplex-Heterogeneous (RWR-MH) networks. MultiVERSE is a fast and scalable method to learn node embeddings from multiplex and multiplex-heterogeneous networks. We evaluate MultiVERSE on several biological and social networks and demonstrate its performance. MultiVERSE indeed outperforms most of the other methods in the tasks of link prediction and network reconstruction for multiplex network embedding, and is also efficient in link prediction for multiplex-heterogeneous network embedding. Finally, we apply MultiVERSE to study rare disease-gene associations using link prediction and clustering. MultiVERSE is freely available on github at https://github.com/Lpiol/MultiVERSE.


Introduction
Networks are powerful representations to describe, visualize, and analyse complex systems in many domains.Recently, machine learning techniques started to be used on networks, but these techniques have been developed for vector data and cannot be directly applied.A major challenge thus pertains to the encoding of high-dimensional and graph-based data into a feature vector.Network embedding (also known as graph representation learning) provides a solution to this challenge and allows opening the complete machine learning toolbox for network analysis.
The high efficiency of network embedding approaches has been demonstrated in a wide range of applications such as community detection, node classification, or link prediction.Moreover, network embedding can leverage massive graphs, with millions of nodes [42].Thus, with the explosion of big data, network embeddings have been used to study many different networks, such as social [58], neuronal [62] and molecular networks [72].
So far, network embedding approaches have been mainly applied to monoplex networks (i.e.single networks composed of one type of nodes and edges) [75,38,42].Current technological advances however generate a large spectrum of data, which form large heterogeneous datasets.By design, monoplex networks are not suited to represent such diversity and complexity.Multi-layer networks, including multiplex [51] and multiplex-heterogeneous [91] networks have been proposed to handle these more complex but richer heterogeneous interaction datasets.
Multiplex networks are composed of several interaction layers, each layer being a monoplex network.All the layers share the same set of nodes, but their edges belong to different categories (Figure 1A).Multiplex representation is pertinent to depict the diversity of interactions.For instance, in a molecular multiplex network, the different layers could represent physical interactions between proteins, their belonging to the same molecular complexes or the correlation of expression of the genes across different tissues.Analogously, in social multiplex networks, a person can belong to different network layers describing different types of relationships, such as friendships or common interests.
A heterogeneous network is a multi-layer network in which each layer is a monoplex network with its specific type of nodes and edges (Figure 1B).The two monoplex networks are connected by bipartite interactions, i.e. edges linking the different types of nodes belonging to the two monoplex networks.Such heterogeneous networks have been studied in different research fields.For example, in network medicine, a drug-protein target heterogeneous network has been constructed with a drug-drug similarity monoplex network, a proteinprotein interaction monoplex network and bipartite interactions between drugs

Related work in network embedding
Network embedding relies on two key components: a similarity measure between pairs of nodes in the original network and a learning algorithm.Given a network and a similarity measure, the aim of network embedding is to learn a vector representation of the nodes in a lower dimension space, while preserving as much as possible the similarity.In the next sections we will present the state-of-the-art in monoplex, multiplex and multiplex-heterogeneous network embedding.
Figure 1: In a multiplex network (A), the different layers share the same set of nodes but different types of edges.In an heterogeneous network (B), the two networks are composed of different types of both nodes and edges, connected by bipartite interactions.In a multiplex-heterogeneous network (C), all the layers of the same multiplex network share the same set of nodes.The different multiplex networks are connected by bipartite interactions between the two types of nodes.For the sake of simplicity, the figure does not represent all the possible bipartite interactions (each layer of a given multiplex is indeed linked with every layer of the other multiplex).

Monoplex network embedding
Many network embedding methods have been recently developed to study a large variety of networks, from biological to social ones.The classical method deepwalk [75] inspired a series of methods such as node2vec [38] and LINE (for Large-scale Information Network Embedding) [87].Deepwalk uses truncated random walks to compute the node similarity in the original network.Then, a combination of the skip-gram learning algorithm [66] and hierarchical softmax [68] is used to learn the graph representations.Skip-gram is a model based on natural language processing.It intends to maximize the probability of co-occurrence of nodes within a walk, focusing on a window, i.e. a section of the path around the node.Node2vec [38] upgrades deepwalk by introducing negative sampling during the learning phase [67].Moreover, node2vec allows biasing the random walks towards depth or breadth-first random walks, in order to tune the exploration of the search space.LINE [87] follows a different approach to optimize the embedding: it computes the node similarity using an adjacency-based proximity measure in association with negative sampling.Some embedding methods are based on matrix-factorization, such as GraRep [10] or HOPE [74], and others on deep neural networks such as graph convolutional networks (GCN) [50].
These embedding methods have been applied to link prediction or node labelling tasks.Their performance rely upon multiple criteria like the size of the network, its density, the embedding dimension or the evaluation metrics [37].Overall, they have been designed to handle monoplex networks.However, we now have access to a richer representation of complex systems as multiplex networks, and some recent methods have explored the embedding of such multiplex networks.

Multiplex network embedding
The most straightforward approach to deal with multiplex networks is to merge the different layers into a monoplex network [7].However, this approach creates a new network with its own topology, and loses the topological features of the individual layers.This new topology is logically biased towards the initial topology of the denser layers [18].Different network embedding methods have been introduced in order to avoid merging multiplex network layers and take advantage of the multiplex structure [98,102,2,97].Overall, these approaches are based on truncated random walks to compute the similarity in the multiplex network.Ohmnet [102] relies on node2vec [38] and requires the definition of a hierarchy of layers to model dependencies between the layers.This additional information is an advantage when the layer hierarchy is known or easy to establish, but this is not the case for many multiplex networks, such as social or molecular networks.In 2018, another multiplex network embedding method, named Scalable Multiplex Network Embedding (MNE), have been developed [98].This approach is also based on node2vec [38].For each network node, it extracts one high-dimensional common embedding shared across all the layers of the multiplex network.In addition, MNE computes a lower-dimensional embedding for every node in each layer of the multiplex network.Another method based on node2vec is multi-node2vec [97].In multi-node2vec, the multiplex embedding is constructed by the random walks jumping from one layer to another.Multi-Net [2] also proposes a random walks procedure in the multiplex network, inspired from [39].Similarly to multi-node2vec, the random walks can jump from one layer to another.Multi-Net learns the embeddings using stochastic gradient descent.The performances of Ohmnet [102], Multi-net [2] and MNE [98] have been compared in the context of network reconstruction [2].In this task, the aim is to reconstruct one layer of the multiplex network from the embeddings of the other layers.The results show better performances for Multi-net on a set of social and biological multiplex networks [2].

Multiplex-heterogeneous network embedding
Some methods can perform the embedding of heterogeneous networks [25,83].One famous approach is metapath2vec [25].It extends skip-gram to learn node embeddings for heterogeneous networks using meta-paths, which are predefined composite relations between different types of nodes.For instance, in the context of a drug-protein target heterogeneous network, the meta-path drug-protein target-drug in the network could bias the random walks to extract the information related to drug combinations.
Nevertheless, to our knowledge, no approach is specifically dedicated to the embedding of different types of nodes from multiplex-heterogeneous networks.In the next section, we present formally MultiVERSE, a new method for multiplex and multiplex-heterogeneous network embedding relying on VERSE [89] and coupled with Random Walks with Restart extended to Multiplex (RWR-M) and Multiplex-Heterogeneous graphs (RWR-MH) [91].

MultiVERSE
In this section, we define the key components of MultiVERSE: the VERSE general framework in the context of multiplex graphs, the learning objective, and our particular implementation with Random Walk with Restart for Multiplex networks (RWR-M) and Random Walk with Restart for Multiplex-Heterogeneous networks (RWR-MH) (Figure 2).We finally describe the MultiVERSE algorithm.

VERSE: a general framework for network embedding
The aim of VERSE network embedding is to learn a low-dimensional nonlinear representation w i of the nodes v i to a d-dimensional continuous vector, where d < n, using Kullback-Leibler optimization [89].We denote d the dimension of the embedding space, and n the dimension of the adjacency matrix.VERSE was originally developed for the embedding of monoplex networks [89].The VERSE framework is nevertheless general and versatile enough to be expanded to multiplex and multiplex-heterogeneous networks, as shown in the following sections.

Similarity distributions
Consider an undirected graph G = (V, E) with V = {v i , i = 1, . . ., n} the set of nodes (|V | = n), and E ⊆ V × V the set of edges, and sim G : Hence, the similarity for any node v is expressed as a probability distribution sim G (v, .).We note w i the vector representation of node i in the embedding space (W is a (n × d)-matrix).
The (non-normalized) distance or similarity between two nodes embeddings w u and w v is defined as the dot product w u • w T v .Using the softmax function, we obtain the normalized similarity distribution in the embedding or vector space: Finally, the output of any network embedding method is a matrix of embeddings W such as, ∀v ∈ V , sim Emb (v, .):sim Emb (v, .)≈ sim G (v, .).This requires a learning phase, which is described in the next section.

Learning objective
This step updates the embeddings at each iteration in order to project sim G into the embedding space leading to the preservation of the topological structure of the graph.In the framework of VERSE, as sim E and sim G are both probability distributions, this optimization phase aims to minimize the Kullback-Leibler divergence (KL-divergence) between these two similarities: We can keep only the parts related to sim Emb as it is the target to optimize and sim G is constant.This leads to the following objective function: sim Emb is defined as a softmax function and needs to be normalized over all the nodes of the graph at each iteration, which is computationally heavy.Therefore, following the VERSE algorithm [89], we used Noise Contrastive Estimation (NCE) to compute this objective function [40,69].NCE trains a binary classifier to distinguish node samples coming from the distribution of similarity in the graph sim G and those generated by a noise distribution Q.We define D as the random variable representing the classes, D = 0 for a node if it has been drawn from the noise distribution Q or D = 1 if it has been drawn from the empirical distribution and E is the expected value.With u a node drawn from P and v drawn from sim G (u, .), with NCE we draw s < n negative samples v neg from Q(u).
In this framework, the objective function becomes the negative log-likelihood that we want to minimize via logistic regression: where P W is computed as the sigmoid (σ(x) = (1 + e −x ) −1 ) of the dot product of the embeddings w u and w v , and sim E (u, .) is computed without normalization.It has been proven that the derivative of NCE converges to gradient of cross-entropy when s increases, but in practice small values work well [69].Therefore, we are minimizing the KL-divergence from sim G .
Overall, VERSE is a general framework for network embedding with the only constraint that sim G must be defined as a probability distribution.In this work, we computed sim G using Random Walks with Restart on Multiplex (RWR-M) and Random Walks with Restart on Multiplex-Heterogeneous (RWR-MH) networks [91].We describe this particular implementation in the next section.Let us consider a finite graph, G = (V, E), with adjacency matrix A. In a classical RW, an imaginary particle starts from a given initial node, v 0 .Then, the particle moves to a randomly selected neighbour of v 0 with a probability defined by its degree.We can define p t (v) as the probability for the random walk to be at node v at time t.Therefore, the evolution of the probability distribution, p t = (p t (v)) v∈V , can be described as follows:

Random
where M denotes a transition matrix that is the column normalization of A.
The stationary distribution of Equation ( 6) represents the probability for the particle to be located at a specific node for an infinite amount of time [60].Random Walk with Restart (RWR) additionally allows the particle to jump back to the initial node(s), known as seed(s), with a probability r ∈ (0, 1) at each step.In this case, the stationary distribution can be interpreted as a measure of the proximity between the seed(s) and all the other nodes in the graph.We can formally define RWR by including the restart probability in Equation ( 6): The vector p 0 is the initial probability distribution.Therefore, in p 0 , only the seed(s) have values different from zero.Equation ( 7) can be solved in a iterative way [91].
In our previous work, we expanded the Random Walk with Restart algorithm to Multiplex (RWR-M) and Multiplex-Heterogeneous networks (RWR-MH) [91].Below, we show how the output of RWR-M and RWR-MH can easily be adapted to produce sim G , the required input for the VERSE framework.

Random Walk with Restart on Multiplex networks (RWR-M)
We define a multiplex graph as a set of L undirected graphs, termed layers, which share the same set of n nodes [51,22].The different layers, α = 1, . . ., L, are defined by their respective n × n adjacency matrices, A [α] = (A [α] (i, j)) i,j=1,...,n .A [α] (i, j) = 1 if node i and node j are connected on layer α, and 0 otherwise [3].We do not take into account potential self-interactions and therefore set In addition, we consider that v α i represents the node i in layer α.
Thus, we can represent a multiplex graph by its adjacency matrix: and define it as G M = (V M , E M ), where: RWR-M should ideally explore in parallel all the layers of a multiplex graph to capture as much topological information as possible.Therefore, a particle located in a given node, v α i , may be able to either walk to any of its neighbours within the layer α or to jump to its counterpart node in another layer, v β i with β = α [21].Additionally, the particle can restart in the seed node(s) on any layer of the multiplex graph.In order to match these requirements, we previously defined a multiplex transition matrix and expanded the restart probability vector, allowing us to apply Equation ( 6) on multiplex graphs [91].
In this study, we independently run the RWR-M algorithm n times, using each time a different node as seed.This allows measuring the distance from every single node to all the other nodes in the multiplex graph.This node-to-node distance matrix is actually a probability distribution describing the particle position in the steady state, where therefore fulfilling the requirements of the VERSE input.We set the RWR-M parameters to the same values used in our original study (r = 0.7, τ = (1/L, 1/L, . . ., 1/L), δ = 0.5) [91].

Random Walk with Restart on Multiplex-Heterogeneous networks (RWR-MH)
A heterogeneous graph is composed of two graphs with different types of nodes and edges.In addition, it also contains a bipartite graph in order to link the nodes of different type (bipartite edges) [54].In our previous study [91], we described how to extend the RWR to a graph which is both multiplex and heterogeneous.However, this study considered only one multiplex graph in the multiplex-heterogeneous graph.For the present work, we additionally expanded RWR-MH to a complete multiplex-heterogeneous graph, i.e. both components of the heterogeneous graph can be multiplex (Figure 1, C), based on the work of [27].Let us consider a L-layers multiplex graph, We also define a second L-layers multiplex graph, with m×L nodes, U M = u α j , j = 1, . . ., m, α = 1, . . ., L .We additionally need a bipartite graph The edges of the bipartite graph only connect pairs of nodes from the different sets of nodes, V M and U M .It is to note that the bipartite edges should link nodes with every layer of the multiplex graphs.We therefore need L identical bipartite graphs, G B ) to define the multiplex-heterogeneous graph.We can then describe a multiplex-heterogeneous graph, G M H = (V M H , E M H ), as: In the RWR-MH algorithm, the particle should be allowed to move in any of the multiplex graphs as described in the RWR-M section.In addition, it may be able to jump from a node in one multiplex graph to the other multiplex graph following a bipartite edge.We also have to bear in mind that the particle could now restart in different types of node(s), i.e. we can have seed(s) of different category (see Figure 1, C).We accordingly defined a multiplex-heterogeneous transition matrix and expanded the restart probability vector.This gave us the opportunity to extent and apply Equation ( 6) on multiplex-heterogeneous graphs [91,27].
In the context of MultiVERSE, we independently run the RWR-MH algorithm n + m times.In each execution, we select a different seed node until all the nodes from both multiplex graphs have been used as individual seeds.As a result, we can define a node-to-node distance matrix matching VERSE input criteria, i.

MultiVERSE algorithm
Algorithm 1 presents the pseudo-code of MultiVERSE based on the random walks with restart on multiplex and multiplex-heterogeneous networks [91] and Our implementation of VERSE with NCE is slightly different from the original one.We perform sequentially first the RWR-M or RWR-MH for all the nodes of the graphs in order to obtain the similarity distribution sim G M .The output of this step is the probability matrix p, where p u is the probability vector representing the similarities between u and all the other nodes.The matrix of the embedded representation of the nodes, W , is randomly initialized.For each iteration, from one node u sampled randomly from a uniform distribution U, we filter the probability vector p u .We keep the N max highest probabilities because the shape of the distribution of probabilities falls very fast to very low probabilities.Doing so, we can speed up the calculation by filtering out the lowest probabilities.We normalize this resulting probability vector ṗu , and sample one node v according to its probability in p u .We set empirically the parameter N max = 300 for large graphs (number of nodes exceeding 5000).For smaller graphs, we set this parameter to 10% − 20% of the number of nodes of the graphs, depending on the shape of the distribution.These two steps (lines 6 and 7) were not in the original VERSE.We parallelized the repeat loop (line 4) and added a parallelized for loop after line 5 in order to run the code from line 6 to 12 in parallel P times.In our simulations, we set P = 100.
Then, we update W u and W v according to algorithm 2 by reducing their distances in the embedding space.We added the bias for NCE: bias pos = log(N ) and bias neg = log(N/s).
Then, s negative nodes are sampled from Q(u) and we update the corresponding embeddings by increasing their distances in the embedding space.The update can also be seen as the training part with lr as the learning rate of the for i=1, ... , s do 10: end for 13: until Maximum step reached binary classifier of the NCE estimation as described in equation 5.The whole process is repeated until the maximum steps are reached.MultiVERSE is freely available on github at https://github.com/Lpiol/MultiVERSE.

Experimental setup
We propose a benchmark to compare the performance of MultiVERSE and other embedding methods for multiplex and multiplex-heterogeneous networks.The performances are evaluated with link prediction for both multiplex and multiplex-heterogeneous networks, and with network reconstruction for multiplex networks.

Evaluation of multiplex network embedding
In the next sections, we describe the datasets, the evaluation tasks and the methods used for comparison.

Multiplex network datasets
We used 7 multiplex networks (4 biological, 3 social) to evaluate the different approaches of multiplex network embedding.The networks CKM, LAZEGA, C.ELE, ARXIV, and HOMO have been extracted from the CoMuNe lab database https://comunelab.fbk.eu/data.php.We constructed the other two networks, DIS and MOL.
• CKM physician innovation (CKM): This multiplex network describes how physicians in four towns in Illinois used the new drug tetracycline [14].
It is composed of 3 layers corresponding to three questions asked to the physicians: i) to whom do you usually turn when you need information or advice about questions of therapy?ii) who are the three or four physicians with whom you most often find yourself discussing cases or therapy in the course of an ordinary week -last week for instance?iii) would you tell me the first names of your three friends whom you see most often socially?
• Lazega network (LAZEGA): This multiplex social network is composed of 3 layers based on co-working, friendship and advice between partners and associates of a corporate law partnership [29].
• ArXiv network (ARXIV): This multiplex network is composed of 8 layers corresponding to different ArXiv categories.The dataset has been restricted to papers with 'networks' in the title or abstract, up to may 2014 [19].The original data from the CoMuNe Lab database is divided in 13 layers.We extracted the 8 layers (1-2-3-5-6-8-11-12) having more than 1000 edges.
• Disease multiplex network (DIS): We constructed a disease multiplex network composed of 3 layers: i) A disease-disease network based on a projection of a disease-drug network from the Comparative Toxicogenomics Database (CTD) [16] extracted from BioSNAP [103].In this network, an edge between two diseases is created if the Jaccard Index between the neighborhoods of the two nodes in the original bipartite network is superior to 0.4.Two diseases are thereby linked if they share a similar set of drugs.This projection has been done using NetworkX [41].ii) A disease-disease network where the edges are based on shared symptoms.The network has been constructed from the bipartite disease-symptoms network from [101].
Similarly to [101], we use the cosine distance to compute the symptom-based diseases similarity for this network.We kept for the disease-disease network all interactions with a cosine distance superior to 0.5 iii) A comorbidity network from epidemiological data extracted from [45].
• Homo sapiens network (HOMO): We collected 4 layers from the original network on CoMuNe Lab [20], keeping physical association, direct interaction, association and co-localization layers.The data are initially extracted from BioGRID [85] • Human molecular multiplex network (MOL): We also constructed another molecular network composed of 3 layers: i) A protein-protein interaction (PPI) layer corresponding to the fusion of 3 datasets: APID (apid.[77] and corresponding to the human Reactome data [15].iii) A molecular complexes layer constructed from the fusion of Hu.map [26] and Corum [34], using OmniPathR [90].
The number of nodes and edges of the different layers composing these multiplex networks are detailed in table 1.

Methods for comparison
We compare MultiVERSE with 6 methods designed for monoplex network embedding (deepwalk, node2vec, LINE) and multiplex network embedding (Ohmnet, MNE, Multi-node2vec), and 4 link prediction heuristic scores (in the link prediction task).

Network embedding methods
• deepwalk [75]: This method is based on non-biased random walks, and apply the skip-gram algorithm [66] to learn the embeddings.We set the context window to 10, and the number of random walks to start at each node to 10.
• node2vec [38]: This method is an extension of deepwalk with a pair of parameters p and q that will bias the random walks for Breadth-first Sampling or Depth-first Sampling.We set p = 2 and q = 1 to promote moderate explorations of the random walks from a node, as stated in [38], and the other parameters as for deepwalk.
• LINE [87]: LINE is not based on random walks, but computes the similarities using an adjacency-based proximity measure in association with negative sampling.It approximates the first and second order proximities in the network from one node.First order proximity refers to the local pairwise proximity between the vertices in the network (only neighbours), and second order proximity look for nodes sharing many connections.We set the negative ratio to 5.
• OhmNet [102]: This approach takes into account the multi-layer structure of multiplex networks.It is a random walk-based method that uses node2vec to learn the embeddings layer by layer.We applied the same parameters as in node2vec.The user has to define a hierarchy between layers.We created a 2-level hierarchy for all multiplex networks, in the same way as [2].
• MNE [98]: This method is also designed for multiplex networks and uses node2vec to learn the embeddings layer by layer.For each node, MNE computes a high-dimensional common embedding and a lower-dimensional additional embedding for each type of relation of the multiplex network.The final embedding is computed using a weighted sum of these two high-dimensional and low-dimensional embeddings.We used the default parameters.
• Multi-node2vec [97]: This multiplex network embedding method is also based on node2vec.The random walks can jump to different layers and explore in this way the multiplex neighborhood.The length of the random walks is set to 100.
We used OpenNE (https://github.com/thunlp/OpenNE) to implement deepwalk, node2vec and LINE.The other methods have been implemented from the source code associated to the different publications.

Link prediction heuristics
In order to evaluate the relevance of the aforementioned network embedding methods, we also compared them with four classical and straightforward link prediction heuristic scores for node pairs [38].Table 2 provides formal definitions of these heuristic scores.

Score Definition
Jaccard Coefficient (JC) Table 2: Link prediction heuristic scores defined in terms of neighborhood for two nodes u and v.We denote N (u) the set of neighbour nodes of node u ∈ V in the graph G(V, E)

Evaluation tasks
We evaluate the different methods on multiplex networks by measuring their performances in two different tasks: link prediction and network reconstruction.For all the evaluations, we set the embedding dimension to d = 128 as in [38,102,75], and used the package EvalNE v0.3.1 [63].

From node embeddings to edge predictions
MultiVERSE and the other embedding methods allow learning vector representations of nodes from networks.We aim here to test their performance on link prediction and network reconstruction.We hence need to predict whether an edge exists between every pairs of node embeddings.To do so, given two nodes u and v, we define an operator • over the corresponding embeddings f (u) and f (v).This gives a representation g : V × V → R d , with d the dimension of the embeddings, V the set of nodes and g(u, v) = f (u) • f (v).Our test set contains both true and false edges (present and absent edges, respectively).We apply five different operators •: Hadamard, Average, Weighted-L1, Weighted-L2 and Cosine (Table 3)).The outputs of the embedding operators are used to feed a binary classifier for the evaluation tasks.This classifier aims to predict if there is an edge or not between two nodes embeddings.Similarly, we use the output of the four link prediction heuristic scores described in Table 2 with a binary classifier to predict edges in a multiplex network.

Link prediction
We first evaluate the performance of the different methods to predict correctly edges removed from the original multiplex networks (Figure 3).We remove 30% of the links in each layer of the original networks.We applied the Andrei Broder algorithm [8] in order to randomly select the links to be removed while keeping a connected graph in each layer.This step provides the multiplex training network, to which we apply the 3 categories of methods (see Figure 3): • The methods specifically designed for monoplex network embedding (node2vec, deepwalk and LINE) are applied individually on each layer of the multiplex Table 3: Embedding operators used to predict edges in the tasks of link prediction and network reconstruction.The definitions describe the ith components of g(u, v).
networks.We thereby obtain one embedding per layer and average them (arithmetic mean) in order to obtain a single embedding for each node.
We then apply the embedding operators.We refer to these approaches in the results section as node2vec-av, deepwalk-av and LINE-av.gned for multiplex network embedding (Ohmnet, MNE, Multi-node2vec) are applied directly on the training multiplex network.We then apply the embedding operators.
• The link prediction heuristic scores JC, CN, AA and PA are applied individually on each layer of the multiplex networks.We then average the scores, as JC-av, CN-av, AA-av, and PA-av.
From the outputs of the embedding operators and heuristic scores, we feed and train a binary classifier and then test it on the 30% of test edges that have been removed previously.The binary classifier is a logistic regressor.
The evaluation metrics for link prediction is ROC-AUC as it is commonly used for embedding evaluation on link prediction and to validate network embedding [98,38].The ROC-AUC is computed as the area under the ROC curve, which plots the true positive rate (TPR) against the false positive rate (FPR) at various threshold settings.An AUC value of 1 represent a model that classifies perfectly the samples.

Network reconstruction
Network reconstruction is another approach to evaluate network embedding methods [95,2,36].In this case, the goal is to quantify the amount of topological information captured by the embedding methods.This is equivalent to predict if we can go back from the embedding to the original adjacency matrix of each layer of the multiplex graph.
Theoretically, to reconstruct the networks, one would need to apply link prediction to every possible edge in the graphs.This is however in practice not scalable to large graphs.Indeed, it would correspond to n(n − 1)/2 potential edges to classify (for undirected networks of n nodes without self-loops).In addition, the networks in our study are sparse, with much more false (absent) Figure 3: General approach for link prediction on multiplex networks than true (present) edges, leading to large class imbalance.In this context, ROC-AUC can be misleading, as large changes in the ROC Curve or ROC-AUC score can be caused by a small number of correct or incorrect predictions [30].In order to account for class imbalance, we used the precision@K [95].This evaluation metric is based on the sorting in descending order of all predictions and consider the first K best predictions to evaluate how many true edges (the minority class) are predicted correctly by the binary classifier.From the outputs of the embedding operators, we perform network reconstruction by training a binary classifier on a subset of the original networks (Figure 4).We choose a subset of 95% of the edge pairs from the original adjacency matrix of each layer for the smaller multiplex networks (CKM, LAZEGA and C.ELE) to construct the training graph.As the class imbalance increases with the number of nodes and sparsity of the networks, we choose smaller subsets for the largest networks, respectively 5% of edges for the ARXIV network and 2.5% for the other networks, as in previous publications [36,95].For each layer, K is defined as the maximum of true edges in this subset of edge pairs.We use a Random Forest algorithm as a binary classifier for network reconstruction, as it is known to be less sensitive to class imbalance [71].In network reconstruction, the results correspond to the training phase of the classifier, there is no test phase.

Multiplex-heterogeneous network datasets
• Gene-disease multiplex-heterogeneous network: We use the two multiplex networks presented in the previous sections: the disease (DIS) and molecular multiplex networks (MOL) (Table 1).In addition, we extracted the curated gene-disease bipartite network from the DisGeNET database [76] in order to connect the two multiplex networks.This bipartite interaction network contains 75445 interactions between 5188 diseases and 9179 genes.We obtain a multiplex-heterogeneous network, as represented in Figure 1C.
• Drug-target multiplex-heterogeneous network: We use the same molecular multiplex network (MOL) from the previous multiplex-heterogeneous network.We constructed the following 3-layers drug multiplex network: (i) the first layer (2795 edges, 877 nodes) has been extracted from Bionetdata (https://rdrr.io/cran/bionetdata/man/DD.chem.data.html)and the edges correspond to Tanimoto chemical similarities between drugs if superior to 0.6, (ii) the second layer (678 edges, 362 nodes) comes from [12] and the edges are based on drug combinations as reported in clinical data, (iii) the third layer (13397 edges, 658 nodes) is the adverse drug-drug interactions network available in [12].The drug-target bipartite network has been extracted from the same publication [12], and contains 15030 bipartite interactions between 4412 drugs and 2255 protein targets.

Evaluation task
We validate the multiplex-heterogeneous network embedding using link prediction.We remove randomly 30% of the edges but only from the bipartite interactions to obtain a training graph.We then train a Random Forest on the training graph, and test on the 30% removed edges.Based on the multiplexheterogeneous networks described previously, the idea behind this evaluation is to test if we can predict gene-disease and drug-gene links.However, comparisons with other approaches are not possible as, to our knowledge, no existing multiplex-heterogeneous network embedding method are currently available in the literature.

Link prediction
Our aim for this case-study is to predict new gene-disease links.We thereby applied MultiVERSE on the full gene-disease multiplex-heterogeneous network without removing any edges, and trained a binary classifier (Random Forest) using edges from the bipartite interactions.Then, we test all possible gene-disease edges that are not in the original bipartite interactions and involve Progeria and Xeroderma pigmentosum VII diseases.Finally, we select the top 5 new gene-disease associations for each disease.

Clustering
We also applied MultiVERSE to the gene-disease multiplex-heterogeneous genedisease network, followed by spherical K-means [9] to cluster the vector representations of nodes.K-means clustering is an approach based on cosine similarity and well-adapted to high-dimensional clustering [99].We define the number of clusters for spherical K-means to 350, in order to obtain cluster sizes that can be analysed from a biological point of view.

Link prediction
We evaluate the performance of the different methods (link prediction heuristics and network embedding) on the task of link prediction applied to the set of multiplex networks.First, we can observe that the heuristics are not efficient for link prediction, with ROC-AUC only slightly better than random classification (Table 4).The methods based on embedding always perform better than these baselines.In addition, the ROC-AUC is in most of the cases higher when the models take into account the multiplex network structure rather than the monoplex-average, as observed in [98].For instance, using the Hadamard operator, the ROC-AUC average over all the networks of the three monoplex-average approaches (node2vec-av, deepwalk-av, LINE-av) is 0.8025, whereas the average of the three multiplex-based approaches (Ohmnet, MNE, Multi-node2vec) is 0.8381.The ROC-AUC score average of MultiVERSE in this context is 0.9011.Nevertheless, node2vec-av and deepwalk-av perform very well and even outperform multiplexbased approaches on various scenarios, for instance on the C.ELE and ARXIV networks.
MultiVERSE combined with the Hadamard operator outperforms the other methods for all the tested networks but CKM.In addition, MultiVERSE is the best approach when combined with three out of five operators (Hadamard,  4: ROC-AUC scores for link prediction on the 7 reference multiplex networks, for link prediction heuristics (CN-av, AA-av, JC-av, PA-av) and network embedding methods combined with different operators (Hadamard, Weighted-L1, Weighted-L2, Average, and Cosine).For each network, the best score is in bold; for each operator, the best scores are underlined.
Average, Cosine).These results suggest that RWR-M is able to better capture the topological features of the networks under study.

Network reconstruction
We next evaluate the performances of the different embedding methods on the task of network reconstruction applied to multiplex networks.As described in section 4.1.3,we now rely on the evaluation metric, precision@K.The experimental results are shown in Table 5.
On one hand, for the small networks (i.e., CKM, LAZEGA and C.ELE), the best precision is achieved with LINE-av in combination with any of the operators but Cosine.In particular, LINE-av obtains a perfect score for the CKM network using the Weighted-L2 or Hadamard operators.MNE is in second position with more than 99% of precision using the Weighted-L1 or Weighted-L2 operators.LINE-av also presents good performances for the C.ELE network with a precision of 93.67% using the Weighted-L2 operator, almost 20% higher than the second best method on this network (Multi-node2vec with a score of 0.7568 using the Weighted-L2 operator).
On the other hand, we can group together the results obtained for the large networks (DIS, ARXIV, HOMO and MOL).In this case, MultiVERSE achieves the best performance in combination with different operators.Large networks are sparse, leading to high class imbalance (4.1.3).Still, MultiVERSE achieves a good score for the HOMO and DIS networks, with precision@K of 0.8729 and 0.6784, respectively.The precision obtained on the molecular network (MOL) is the lowest, with a precision@K of 0.4143.The complexity of the task is possibly higher as the number of nodes and class imbalance increase.
Overall, the lowest scores are obtained by MNE and, in general, the Cosine operator performs poorly for all methods.The network reconstruction process is a complex task, and the performance depends on the size and density of the different layers composing the multiplex network.Nevertheless, MultiVERSE obtains good results for most of the networks without any processing of the imbalanced data.

Evaluation results for multiplex-heterogeneous network embedding
The task of link prediction on multiplex-heterogeneous networks is applied to MultiVERSE only, as to our knowledge no other methods exist for the embedding of multiple nodes from multiplex-heterogeneous networks.MultiVERSE has a score of ROC-AUC superior to 0.9 with the Hadamard and Average operators (Table 6), meaning that the method can predict with high precision the genedisease and drug-disease links from the corresponding multiplex-heterogeneous networks.
5.3 Case study results: discovery of new gene-disease associations

Discovery of new gene-disease associations with link prediction
The results of the evaluations on multiplex-heterogeneous network link prediction show that MultiVERSE combined with the Hadamard and Average operators reach ROC-AUC scores superior to 0.9 (Table 6).We here investigate in detail the top 5 new gene-disease associations predicted by MultiVERSE combined with these operators for Hutchinson-Gilford Progeria Syndrome (HGPS) and Xeroderma pigmentosum VII (Table 7) .MultiVERSE top predictions reveal interesting candidate genes, including NOS2, TNF, KRAS, CDC42BPB, and HSPB1.NOS2 is predicted by both operators.It encodes a nitric oxide synthase expressed in liver.It has been associated with longevity [70].TNF is a member of the tumor necrosis factor superfamily, and produces a multifunctional proinflammatory cytokine.TNF is also known to be involved in aging [17] and has been previously linked to Progeria [73].KRAS is an oncogene of the RAS family, and, as LMNA, it is farnesylated [6].CDC42BPB produces a serine/threonine protein kinase, and is a downstream effector of CDC42.CDC42BPB is linked to Myotonic Dystrophy [56], which presents some progeroid features [65].This gene has also been recently related to a neurodevelopmental disorder associated to hypolasia [13], one of the major phenotype of HGPS.Finally, HSPB1 encodes a small heat shock protein involved in stress resistance and actin organization [49].HSPB1 is also linked to Charcot-Marie Tooth disease [43], as LMNA [4].

Xeroderma pigmentosum VII
Xeroderma Pigmentosum (XP) is characterized by extreme sensitivity to sunlight, resulting in sunburns, pigment changes in the skin and a greatly elevated incidence of skin cancers.It is a genetically heterogeneous autosomal recessive disorder.Several types exist, and the MeSH term Xeroderma pigmentosum VII corresponds to the group G, caused by mutations in the ERCC5 gene, and with symptoms that overlap Cockayne syndrome [52,93].
MultiVERSE identified various interesting candidates for this disease (Table 7).PTEN is a tumor suppressor mutated in a large number of cancers.It has been linked to melanoma and XP [64].TP53 also a tumor suppressor and is related to many cancers, in particular skin cancers [33].It has also been associated to XP [81].IFNG encodes a cytokine that is a part of the type II interferon class.Altered immune reactivity has been shown in XP patients, and the production of interferon-gamma in peripheral blood cells is impaired in the absence of UV exposure for XP patients [31].The gene SOD1 encodes the superoxide dismutase, which has the function to destroy free superoxide radicals.SOD1 has been associated to spontaneous skin damage and delays wound healing in SOD1-deficient mices [44].ERCC2 produces the XPD protein, involved in XP group D [88].FOXC1 is a member of Forkhead box family transcription factor.It is a DNA binding transcriptional factor that is involved in a wide range of cellular and developmental processes including skin development [5].In addition, FOXC1 is involved in melanoma [96], one of the major phenotypes of XP.POMC encodes a protein called Pro-opiomelanocortin, also related to melanoma [82].

Discovery of new gene-disease associations with clustering
Another illustration of the multiplex-heterogeneous network embedding is clustering.To this goal, we identify clusters with K-means (see section 4.3), and focus more particularly on the clusters containing HGPS and Xeroderma pigmentosum VII disease nodes.Clustering is particularly interesting as it can be applied directly on the embeddings without supervised training.This is also in line with the results from [72] where clustering from embeddings outperforms the other methods for the detection of biological communities.

Cluster containing the HGPS disease node
The cluster containing the HGPS disease node (see Figure 5) contains the LMNA gene node.LMNA mutations have been observed in many different diseases that also belong to the identified cluster, for example the Heart-hand syndrome (Slovenian type) [79], the Charcot-Marie-Tooth disease type I [24], LMNA-related muscular distrophy, and different cardiac diseases related to LMNA mutations [92].
We also analysed the list of genes in the cluster with g:Profiler ([78], default parameters).We found significant enrichments in several annotations related to cardiac functions, such as regulation of heart rate and contraction.The most significant enrichment is cardiac muscle contraction and involves the following genes: GAA, NOS1AP, PLN, RYR2, SCN1B, SCN4B, SCN5A, TNNT2, and CASQ2.In HGPS, death occurs typically as a result of cardiac or cerebrovascular complications [35].

Cluster containing the Xeroderma pigmentosum VII disease node
We also analysed the cluster containing the Xeroderma pigmentosum VII disease node (Figure 6).The cluster contains different diseases, including sunburns, Xeroderma pigmentosum type IV, trichothiodystrophies and Pena-Shokeir syndrome.All these diseases have been previously related to Xeroderma phenotypes [52,53,48].
Several genes known for their implications in XP are present in the cluster, such as ERCC1 to ERCC6 or XPC [86].Using the complete list of genes in the cluster as an input for g:Profiler ( [78], default parameters), we identified several significantly enriched annotations.Among them, we can cite nucleotide-excision DNA repair or response to ultraviolet radiation.XP patients show important impairments in these biological processes [52].
CD151 is another gene of interest.It has no edges with the genes known to be mutated in XP, but has been shown to be involved in skin repair [32], squamous cell carcinoma [57] and is essential for the correct assembly of human basement membranes in kidney and skin [47].In addition, CD151 could be involved in inner ear development [47].Xeroderma is associated to inner ear malfunction [94].

Discussion and conclusion
We present in this study MultiVERSE, a new approach for multiplex and multiplex-heterogeneous network embedding.MultiVERSE is fully parallelized and scalable, even if the current implementation requires the generation of dense matrices, which can raise memory issues when dealing with very large networks.
For multiplex network embedding, we compared MultiVERSE with state-ofthe-art methods using link prediction and network reconstruction.We show that MultiVERSE outperforms various methods specifically developed for multiplex network embedding.As RWR-M applies a random walk in pseudo-infinite time, it might allow MultiVERSE to effectively capture node properties and a better representation of the topological structure of the multiplex network.
A natural extension of this work would be to consider multiplex networks composed of both directed and undirected layers.In a biological context, this would allow considering metabolic and signalling pathways networks into a multiplex structure without losing the information about the information flow.In addition, for the optimization phase, we set a neighborhood parameter N max that depends on the size of the network.A potential improvement could be to develop an adaptive version of the parameter N max that would depend on node topological properties.
For multiplex-heterogeneous network embedding, MultiVERSE allows the embedding different types of nodes.We demonstrate its efficiency for link prediction and illustrates its usefulness for the study of gene-disease associations.We here limited the multiplex-heterogeneous network to two multiplex and one bipartite network.Another natural extension of our work would be to generalize RWR for multiplex-heterogeneous for n multiplex networks and n(n − 1)/2 bipartite linking them (n ∈ N).Doing so, one could easily integrate many different types of nodes.The previous discussion about directed networks is in addition also valid for multiplex-heterogeneous network embedding.
By integrating different types of edges for multiplex network embedding or by integrating different types of both edges and nodes for multiplex-heterogeneous network embedding, MultiVERSE could have a wide variety of applications in diverse domains such as network biology and medicine, social science, computer science, neuroscience or physics.Our illustration of MultiVERSE embedding to study gene-disease associations could easily be applied to drug repositioning and drug discovery, for instance with a multiplex drug-drug network, a drug-target bipartite and a molecular multiplex.In this way, genes, diseases and drugs could be projected in the same vector space for further studies.In neuroscience, multiplex-heterogeneous network embedding could be applied to study the links between genes and neurons [1].In social science, multiplex networks are gaining interest to understand human behaviour [84].Multiplex-heterogeneous network embedding could give insights on epidemic spread [46,59], socio-economic systems [80] or socio-ecological systems [55].comprehensive resource of mammalian protein complexes-2019.Nucleic acids research 47, D1 (2019), D559-D563.

Figure 5 :
Figure 5: Cluster containing the HGPS disease node.Disease-Disease edges from the disease multiplex network are represented in green (shared symptoms) and blue (CTD projection).Gene-Gene edges from the molecular multiplex network are represented in pink (Reactome pathways), red (protein-protein interactions) and orange (molecular complexes).Gene-Disease bipartite interactions are represented with black dashed lines.

Figure 6 :
Figure 6: Cluster containing the Xeroderma pigmentosum VII disease node.Disease-Disease edges from the disease multiplex network are represented in blue (CTD projection).Gene-Gene edges from the molecular multiplex network are represented in pink (Reactome pathways), red (protein-protein interactions) and orange (molecular complexes).Gene-Disease bipartite interactions are represented with black dashed lines.

Table 5 :
precision@K scores for network reconstruction on the 7 reference multiplex networks, for the network embedding methods combined with different embeddings operators (Hadamard, Weighted-L1, Weighted-L2, Average, and Cosine).For each network, the best score is in bold; for each operator, the best score is underlined.The percentage of edges used for the reconstruction is indicated under the name of the network.

Table 6 :
ROC-AUC scores for link prediction using MultiVERSE on 2 multiplexheterogeneous reference networks.Link predictions are computed for the bipartite interactions of the multiplex-heterogeneous networks.The scores higher than 0.9 are highlighted in bold.

Table 7 :
Top 5 predictions of new gene-disease associations for HGPS and Xeroderma pigmentosum VII by MultiVERSE combined with Average and Hadamard operators.