Capucine Jourdain De Muizon

M Srid'evi

Ramanoudjame

Lucie Esteoulle

Claude Ling

Germain Brou

Nicolas Anton

Thierry Vandamme

Marc -Andr'e Delsuc

Dominique Bonnet
email: dbonnet@unistra.fr

Bruno Kieffer
email: kieffer@igbmc.fr

] C Ling

M A Delsuc

] S M Ramanoudjame

T Vadamme

Self-organization properties of a GPCR binding peptide incorporating fluorinated tail studied by fluorine NMR

Keywords: Fluorine NMR, GPCR ligand, Micellization, apelin-17, peptide

The conjugation of the bioactive apelin-17 peptide with a fluorocarbon chain results in higher potency and self-organisational property of the peptide into micelles. Fluorine NMR studies show that fluoropeptide's micelles are mono-disperse while proton NMR indicates that the peptide moiety remains largely disordered despite micellization. A very fast exchange rate is measured between the free and micellar states of the peptide enabling to estimate the number of molecules present in the micelle as two hundred, in agreement with values found by Dynamic Light Scattering measurements.

 The substitution of natural amino-acids by their fluorinated variants or the addition of trifluoromethyl groups at the protein surface have indeed recently expanded the panel of protein NMR

technics to provide insights on G Protein Coupled Receptor (GPCR) function, [3] investigate protein-ligand interactions, [4] measure an enzymatic activity [5] or observe protein folding inside living Escherichia coli cells. [6] While these methods see fluorine as a label and assume a minimal effect of the fluorine substitution on the protein state, other applications take advantage of the fluorine properties to perturb the conformation and dynamics of the molecular system in order to gain insights on the related molecular mechanisms. One compelling example is prolines with one of the two geminal Hγ hydrogens substituted by a fluorine atom, a modification that induces a bias of the ring conformation and the cistrans equilibrium of the X-Pro peptide bond. [7] This property has been used to investigate the role of hydroxyproline in the self-association mechanism of collagen [8] and the human amyloidogenic β2-microglobulin. [9] Another way to modulate biological activities of peptides and proteins is the conjugation with a fluorocarbon chain (F-chain) to one extremity. Indeed, due to their very specific properties and their ability to self-assemble in aqueous solution, the addition of F-chains to organic molecules changes their physical and chemical properties. Thereby, the conjugation of F-chain to insulin has been shown to have an impact on the native hexamer structure and could represent a new way to control self-assembly of proteins and their in vivo stabilities. [10] The addition of F-chain to peptides was applied in the development of vaccines as well. The authors recently demonstrated that the presence of F-chain in peptides derived from HIV, influenza and hepatitis C virus increased their cellular immunogenicity by controlling their ability to form micelles in solution. [11] They also showed that the conjugation of therapeutic peptides targeting GPCR with F-chain provided a promising approach to enhance their plasma stability and in vivo activity. [12] This approach has been successfully applied to apelin-17. This peptide together with its cognate GPCR receptor dubbed APJ [13] has emerged as novel targets for the treatment of cardiovascular diseases. [14] However, the low stability of apelin-17 in plasma precluded its therapeutic development. To circumvent this limitation, a perfluoroundecanoic acyl group (F8-chain) was incorporated to the N-terminal part of apelin-17 (LIT01-196) resulting in a dramatic enhancement of its half-life in plasma, as well as its efficacy on the cardiac contractility and arterial blood pressure on rat model. [12] For all those applications, it is crucial to characterize and to monitor the self-organization of fluoropeptides in order to strictly control their biological activities.

Thanks to its ability to describe both dynamical and structural features of micelles without requesting any chemical modifications, fluorine NMR is best suited to study the properties of fluorosurfactants. [15] Hence, NMR was successfully used to study micellization of gemini surfactants or perfluorooctylsulfonates. [16] In this paper, using Fluorine NMR as a main investigation tool, we report for the first time on the influence of F8-chain addition to apelin-17 on both the structural and self-assembling properties of the resulting fluoroapelin-17 (F8-Ap). To assess the potential structural changes induced by the addition of the F8-chain at the N-terminal position of apelin-17, preliminary proton NMR spectra of both peptides dissolved in PBS buffer were first recorded at 100 µM and 420 µM of peptide concentrations. At 420 µM, the resonances of F8-Ap display significant line broadening while those of apelin-17 remains sharp (data not shown). In order to compare structural propensities of the two peptides, homonuclear ROESY, TOCSY and 1 H-13 C HSQC spectra of apelin-17 and F8-Ap were recorded at 2 mM in water and 298 K, where no line broadening is observed. Frequencies of proton and carbon nuclei were obtained using sequential assignment (Figure S1 and Table S1). The chemical shifts of H α , H β and the corresponding carbon nuclei that offer sensitive probes of the local peptide conformation were used to characterize the potential conformational differences between the two peptides (Figure 1). The condensation of the F8-chain to apelin-17 induces significant chemical shifts changes for the first four N-terminal residues (KFFR, Figure 1A), while the C-terminal display similar chemical shifts.

The differences observed for the first four residues are due to the substitution of the free NH3 group by a peptide bond and probably a reorientation of the aromatic ring of the phenylalanine upon fluorine-pi interaction with the F8chain. [17] The chemical shifts of other amino acids measured for both peptides are very similar and close to theoretical values obtained for a random coil peptide as computed by the Potenci software [18] indicating a lack of secondary structure. Those results are in accordance with previous NMR studies on the structure of apelin-17 by NMR highlighting a predominantly random coil character of the peptide. [19] Thereby, we show here that the incorporation of the F8-chain onto peptide does not impact the structure of apelin-17. The self-association property of the F8-Ap was then studied by recording fluorine spectra for peptide concentrations ranging from 30 µM to 1.3 mM in phosphate buffer. Across this range of concentrations, the fluorine peaks experienced a constant broadening and an additional set of resonances arose for concentrations above 200 µM, indicating the presence of a second population of molecules in the solution (Figure 2). Owing to their significantly larger line-widths, these resonances were assigned to a pool of peptides organized into micelles. Integrating the trifluoromethyl resonances at -82.8 and -80.66 ppm was used to measure the relative fractions of the self-organized and the monomer forms of the peptides, respectively. Assuming that the fraction of the monomers f mono is described by the following fraction

f mono = CMC C0 (1)
where C 0 is the total concentration of peptide and CMC, the Critical Micellar Concentration. [16b] The latter is usually determined from the dependence of the monomeric fraction upon the inverse of the total peptide concentration (Figure S2). This analysis led to a CMC value of 80 µM, in full agreement with the CMC of 84 µM determined by surface tension using a drop tensiometer (Figure S3). This value is significantly lower than the CMC reported for perfluorooctylsulfonates surfactants (in the mM range), indicating a contribution of the peptide moiety to the interactions within the micelle. Indeed, we observed a large dependence of the fluorine line broadening upon phosphate concentration as expected from the sensitivity of self-association properties of surfactant molecules on ionic strength. [20] In the particular case of F8-Ap, screening electrostatic interactions between the positively charged residues neighbouring the F8 chain may be critical to modulate the peptide self-association properties.

The observation of cross-peaks between the two sets of resonances in a fluorine NOESY spectrum recorded at a concentration of 200 µM of F8-Ap indicated that exchange between the monomeric and the micellar states occured in the millisecond time scale (Data not shown). The F8 fluoromethyl group experi-enced a chemical shift variation of 2.1 ppm (1200 Hz) between the two states. This remarkable shift enabled the selective manipulation of the corresponding populations of spins using frequency selective soft pulses (Figure S3). A selective inversion of the fluoromethyl broad peak at -82.8 ppm resulted into a fast decay of the resonance at -80.66 ppm corresponding to the monomeric form of the peptide (Figure 3). After 10 ms, the signal became negative before a slow recovery, indicating that the magnetization transfer between the micelle and the monomer resonances through chemical exchange is significantly faster that the longitudinal relaxation rate (R 1 ) of the fluoromethyl resonance of the monomeric peptide.

Figure 3.

Selective inversion of the trifluoromethyl signal at -82.8 ppm (blue dots) leads to fast magnetization transfer to the signal at -80.5 ppm (red dots) corresponding to the tri-fluoromethyl of the free peptide. Fitting of the peak intensities using a three exponentials model allows the determination of the selective R 1 of the fluoromethyl in both states and the exchange rate between these states Fitting the time evolution of signal intensities with a three exponentials model accounting for two specific longitudinal relaxation and the chemical exchange rates provided values of 173 ± 5 s -1 for the exchange rate between the monomeric and the micellar forms, and 7.5 ± 1.8 s -1 and 1.9 ± 0.6 s -1 for fluoromethyl's R 1 in the monomeric and micellar environment, respectively. As expected, the rate of chemical exchange is significantly slower than frequency difference between the two states (1200 Hz), explaining why no chemical shift difference is observed upon variation of the total peptide concentration. The lower R 1 value observed for the fluoromethyl resonance in the micellar state is consistent with a longer correlation time expected in micelles and the expected evolution of the dipolar relaxation upon increasing correlation time. The micelle dissociation rate obtained from the fit corresponded to a mean lifetime of 15.6 ms, a value significantly longer than the 4.2 µs observed for perfluorooctylsulfonates. [16b] In order to get an insight into the size of the F8-Ap micelles, translational diffusion coefficients were measured by 19 F DOSY NMR. Values of magnetization exchange rates between the micellar and free forms of the peptide suggested a complete averaging of the measured diffusion coefficients during the 150 ms diffusion time of the DOSY while integration of the two trifluoromethyl resonances provided an accurate estimate of the relative populations in both monomeric and micellar states. The 19 F DOSY implementing bipolar gradient pulses were measured for several concentrations of peptide ranging from 100 to 200 µM where the relative fractions of micelles and free peptide undergo a large variation (Figure 4). Consistent with an exchange between two states with distinct diffusion properties, the diffusion coefficient of the free form displayed a linear dependence on the free peptide fractions, while the diffusion of added trifluoroacetic acid TFA remained constant with a value of 1230.0 ± 96 µm 2 .s -1 . As expected from an averaging between two states, the translational diffusion coefficient follows a linear dependence on the relative fraction of free fluoroapelin. The extrapolated values for solutions containing only monomers or micelles provided us with translational diffusion coefficients of 49.8 ± 4 and 289 ± 23 µm 2 .s -1 . Using the Stokes-Einstein expression for the self-diffusion coefficient D of a hard spherical particle of radius R h in a medium of viscosity η

D = k B T 6πηR h (2)
where k B is the Boltzmann constant and T the temperature, a hydrodynamic radius of 49 Å was obtained for the F8-Ap micelles. This value is consistent with the with Dynamic Light Scattering data measured on F8-Ap reporting a micelle's diameter of 10 nm (Figure S5). The ratio of diffusion coefficients obtained for the monomeric and the micellar states of F8-Ap gave a ratio of the hydrodynamic radii of 5.8 ± 1.0. The corresponding ratio of hydrodynamic volumes was 198, providing an estimate of the number of molecules in the micelle. It should be noted that the low peptide concentrations used in this study due to the low CMC value of F8-Ap led to a negligible value of the volume fraction (Φ) of the micelles that usually limits the validity of equation 2 for large solute concentrations. [15] Nevertheless, the use of the third power of the ratio of the diffusion constants is based on serious assumptions that may alter the accuracy of the aggregation factor from NMR diffusion data. First, the shape factor, assumed to be identical for the F8-Ap in the monomeric and micellar forms is probably different for a disordered peptide and an organized micelle. Second, the third power used in our calculation probably overestimates the micelle's size. Using the Flory theory for polymers, it was shown that translational diffusion coefficients measured for globular proteins scales with their size as a power of 2.56 while a power value of 1.7 was found for disordered proteins. [START_REF] Aug'e | Delsuc in NMR Measure of Translational Diffusion and Fractal Dimension[END_REF] Although slightly overestimated, the number of molecules F8-Ap is consistent with the size of the micelles. Indeed, assuming a spherical shape and a partial specific volume of 0.7 cm 3 .g -1 typical of proteins led to an estimate of 172 molecules per micelle. Furthermore, such large aggregation numbers (125) were reported for single chain phospholipid Lyso PG, characterized by a low CMC (18 µM) using translational diffusion measurements by NMR. [22] For F8-Ap micelles, both translational diffusion measurements in the vicinity of CMC and size estimates from DLS suggest significantly larger aggregation numbers than those reported for the sole perfluorooctylsulfonate surfactant (36). [16b] In this work, we described for the first time a structural and dynamical characterization of a bioactive peptide grafted with a fluorocarbon chain at its N-terminal extremity using 19 F NMR methods. On apelin-17, this addition led to self-organisational properties of the peptide into micelles of defined size with little perturbation of the peptide's structure. The self-association of F8-chain is a highly dynamic process that depends on the composition of the solution. Thereby, we suggest that fine-tuning these self-association properties based on 19 F NMR characterization could provide a novel interesting route for designing potent bioactive fluoropeptides.

Experimental Section F8-Ap (LIT01-196) was synthetized as described in [12] . The structural characterization was performed using samples of 3 mg and 2.4 mg of apelin-17 and F8-Ap, respectively, dissolved either in 450 µL of miliQ water and 50 µL D 2 O for NMR field lock. Dissolved samples were subsequently subjected to 5 minutes of vortex and 5 minutes of sonication. Resonance assignments were based on 2D spectra recorded at 298K. Through-bonds connectivities were obtained from TOCSY (80 ms mixing time) and 1 H-13 C HSQC spectra while through-space sequential connectivities were obtained from a ROESY spectrum (300 ms mixing time) on a 700 MHz BRUKER® Avance 3 spectrometer equipped with a TCI cryogenic probe.

Self-association properties of F8-Ap were investigated by dissolving 2.6 mg of dry peptide in 300 µL of PBS buffer at pH 7.4 to get a stock solution at 3 mM. Cascade dilutions were performed to obtain lower concentrations, and each diluted sample was submitted to two cycles of 30 s sonication and 30 s vortex. Fluorine NMR Spectra were recorded using 3 mm tubes on an Avance BRUKER® spectrometer equipped with a cryogenic QCI fluorine probe, operating at 600 MHz proton frequency. The DOSY spectra were recorded using the bipolar-pulse pair stimulated echo pulse sequence [23] with 24 z-gradient amplitudes up to 50 gauss.cm -1 . The diffusion time and gradient pulse length were 150 ms and 2 ms, respectively. The DOSY were processed using PALMA [24] while quantifications of 1D fluorine spectra were performed using SPIKE. [25] Modelling the fluorine signal time evolution after a selective trifluoromethyl signal inversion was performed using a non-linear least square of peak intensities using the Levenberg-Marquart algorithm implemented in the Numpy/Scipy leastsq function.

The F8-Ap, the bioactive apelin-17 peptide grafted with a fluorocarbon chain at its N-terminal extremity presents self-organisational properties. 19 F NMR exchange spectroscopy and translational diffusion measurements is used to investigate this phenomenon. F8-Ap forms micelles of defined size with a fast exchange of peptides between free and micellar states.

Figure 1 .

 1 Figure 1. Chemical structure of fluorocarbon conjugated apelin-17 (F8-Ap) (A). Chemical shift differences between the perfluorinated and the natural Apelin-17 peptide for the Hα (grey) and Cα (black) nuclei measured in water and 298 K (B).

Figure 2 .

 2 Figure 2. 1D 19 F spectrum of F8-Ap in PBS buffer as a function of total peptide concentration C 0 indicated in the right panel. The trifluoromethyl resonance used for quantification is shown in the left panel.

Figure 4 .

 4 Figure 4. Diffusion coefficients of free

Acknowledgements

This work was supported by the Agence Nationale de la Recherche (FluoroPEP, ANR-16-CE18-0030), the French Infrastructure for Integrated Structural Biology (FRISBI, ANR-10-INSB-05-01), Instruct-ERIC, the LabEx MEDALIS (ANR-10-LABX-0034), the Centre National de la Recherche Scientifique, University of Strasbourg. L.E. and S. M. R. was supported by a fellowship from the Minist'ere de l'Education Nationale and ANR.


Léo Pio-Lopez
email: leo.pio.lopez@gmail.com

Alberto Valdeolivas

Laurent Tichit

Élisabeth Remy

Anaïs Baudot

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/ Multi VERSE.

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 highdimensional 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 approaches can exploit very large graphs, with millions of nodes 1 . Thus, with the explosion of big data, network embeddings have been used to study many different networks, such as social 2 , neuronal 3 and molecular networks 4 .

 . Current technological advances however generate a large spectrum of data, which form large heterogeneous datasets. Single monoplex networks are not suited to represent such diversity. Therefore, multi-layer networks, including multiplex 7 and multiplex-heterogeneous 8 networks have been proposed to handle these richer sets of relationships.

). Multiplex representation is pertinent to depict the diversity of interactions between the same nodes. 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 layers describing different types of relationships, such as friendships or common interests.

). 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 protein-protein interaction monoplex network and bipartite interactions between drugs and their target proteins [START_REF] Luo | A network integration approach for drug-target interaction prediction and computational drug repositioning from heterogeneous information[END_REF] . In social science, citation networks are constructed with author-author and document-document monoplex networks connected by author-documents bipartite interactions, as in [START_REF] Zhou | Co-ranking authors and documents in a heterogeneous network[END_REF] .

A multiplex-heterogeneous network is a combination of heterogeneous and multiplex networks by connecting several multiplex networks through bipartite interactions (Fig. 1C). The multiplex-heterogeneous structure is expected to provide a richer view on biological [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] , social [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF] or other real-world systems describing complex relations among different components.

Recently, different studies proposed embedding approaches for multiplex networks [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF][START_REF] Zhang | Scalable multiplex network embedding[END_REF][START_REF] Zitnik | Predicting multicellular function through multi-layer tissue networks[END_REF][START_REF] Wilson | Fast embedding of multilayer networks: An algorithm and application to group fmri[END_REF] and heterogeneous networks [START_REF] Dong | metapath2vec: Scalable representation learning for heterogeneous networks[END_REF][START_REF] Shi | Heterogeneous information network embedding for recommendation[END_REF] . A recent method uses multiplex-heterogeneous information to embed one category of nodes [START_REF] Dursun | Gene Embeddings of Complex network (GECo) and hypertension disease gene classification[END_REF] . However, to our knowledge, no embedding methods are specifically dedicated to the embedding of nodes of different types from multiplex-heterogeneous networks. In this paper, we present MultiVERSE, a fast, scalable and versatile embedding approach to learn node embeddings on multiplex and multiplex-heterogeneous networks. MultiVERSE is based on the VERSE framework [START_REF] Tsitsulin | Versatile graph embeddings from similarity measures[END_REF] , and coupled with Random Walks with Restart on Multiplex (RWR-M) and on Multiplex-heterogeneous (RWR-MH) networks [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] . Our contributions are the following: (C) A multiplex-heterogeneous network composed of two multiplex networks. The multiplex networks are connected by bipartite interactions (dashed lines). For the sake of simplicity, the figure does not represent all the possible bipartite interactions (each layer of a given multiplex is in reality linked with every layer of the other multiplex).

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 vector representations 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 of monoplex, multiplex and multiplexheterogeneous network embedding.

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 [START_REF] Perozzi | Online learning of social representations[END_REF] inspired a series of methods such as node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] and LINE (for Large-scale Information Network Embedding) [START_REF] Tang | Line: Large-scale information network embedding[END_REF] . Deepwalk uses truncated random walks to compute the node similarity in the network. Then, a combination of the skip-gram learning algorithm [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF] and hierarchical softmax [START_REF] Mnih | A scalable hierarchical distributed language model[END_REF] 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 6 upgrades deepwalk by introducing negative sampling during the learning phase [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF] . 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 [START_REF] Tang | Line: Large-scale information network embedding[END_REF] follows a different approach to optimize the embedding: it computes the node similarity using an adjacency-based proximity measure in association with negative sampling. Other embedding methods are based on matrix-factorization, such as GraRep [START_REF] Cao | Grarep: Learning graph representations with global structural information[END_REF] or HOPE [START_REF] Ou | Asymmetric transitivity preserving graph embedding[END_REF] . It has been shown that random-walk based methods for network embedding can be expressed in terms of matrix-factorization [START_REF] Liu | A general view for network embedding as matrix factorization[END_REF] . Another class of methods are based on neural networks such as GraphSAGE [START_REF] Hamilton | Inductive representation learning on large graphs[END_REF] , graph convolutional networks (GCN) [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF] or graph auto-encoders (GAE/VGAE) [START_REF] Kipf | Variational graph auto-encoders[END_REF] .

These embedding methods have been applied to link prediction or node labelling tasks. Their performance rely upon multiple criteria such as the size of the network, its density, the embedding dimension or the evaluation metrics [START_REF] Goyal | Benchmarks for graph embedding evaluation[END_REF] . 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 [START_REF] Boccaletti | The structure and dynamics of multilayer networks[END_REF] . However, this merging 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 [START_REF] De Domenico | The physics of spreading processes in multilayer networks[END_REF] . Different network embedding methods have been introduced in order to avoid merging multiplex network layers and take advantage of the multiplex structure [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF][START_REF] Zhang | Scalable multiplex network embedding[END_REF][START_REF] Zitnik | Predicting multicellular function through multi-layer tissue networks[END_REF][START_REF] Wilson | Fast embedding of multilayer networks: An algorithm and application to group fmri[END_REF] . Overall, these approaches are based on truncated random walks to compute the similarity in the multiplex network. Ohmnet [START_REF] Zitnik | Predicting multicellular function through multi-layer tissue networks[END_REF] relies on node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] and requires the definition of a hierarchy of layers to model dependencies between them. But usually, this layer hierarchy information is not known or easy to establish, particularly for multiplex networks such as social or molecular networks. The Scalable Multiplex Network Embedding (MNE) method [START_REF] Zhang | Scalable multiplex network embedding[END_REF] is also based on node2vec [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] . 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. Multi-node2vec [START_REF] Wilson | Fast embedding of multilayer networks: An algorithm and application to group fmri[END_REF] is another method based on node2vec that constructs the multiplex embedding with the random walks jumping from one layer to another. Multi-Net 11 also proposes a random walks procedure in the multiplex network, inspired from [START_REF] Guo | Levy random walks on multiplex networks[END_REF] . 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 13 , Multi-net [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF] and MNE [START_REF] Zhang | Scalable multiplex network embedding[END_REF] have been compared in the context of network reconstruction [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF] . 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 [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF] .

Multiplex-heterogeneous network embedding. Some methods can perform the embedding of heterogeneous networks [START_REF] Dong | metapath2vec: Scalable representation learning for heterogeneous networks[END_REF][START_REF] Shi | Heterogeneous information network embedding for recommendation[END_REF] . A famous approach is metapath2vec [START_REF] Dong | metapath2vec: Scalable representation learning for heterogeneous networks[END_REF] . 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

MultiVERSE

In this section, we present the key components of MultiVERSE: the VERSE general framework, 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) (Fig. 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 [START_REF] Tsitsulin | Versatile graph embeddings from similarity measures[END_REF] . 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 [START_REF] Tsitsulin | Versatile graph embeddings from similarity measures[END_REF] . The VERSE framework is nevertheless general and versatile enough to be expanded to multiplex and multiplex-heterogeneous networks.

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 : V × V → R a given similarity measure on G such that 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) 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 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 embed- ding space leading to the preservation of the topological structure of the graph. In the framework of VERSE,

(1) ∀v ∈ V , u∈V sim G (v, u) = 1 . (2) sim Emb (v, .) = exp(w v • w T ) n i=1 exp(w v • w i )
. as sim Emb 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 itera- tion, which is computationally heavy. Therefore, following the VERSE algorithm [START_REF] Tsitsulin | Versatile graph embeddings from similarity measures[END_REF] , we used Noise Contrastive Estimation (NCE) to compute this objective function [START_REF] Gutmann | Noise-contrastive estimation: A new estimation principle for unnormalized statistical models[END_REF][START_REF] Mnih | A fast and simple algorithm for training neural probabilistic language models[END_REF] . 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 Emb (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 [START_REF] Mnih | A fast and simple algorithm for training neural probabilistic language models[END_REF] . 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 [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] . We describe this particular implementation in the next section.

Random walk with restart on multiplex and multiplex-heterogeneous networks. Random walk (RW) and random walk with restart (RWR).

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 prob- ability 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:

where M denotes a transition matrix that is the column normalization of A. The stationary distribution of Eq. ( 6) represents the probability for the particle to be located at a specific node when times tends to infinity [START_REF] Lovász | Random walks on graphs: A survey[END_REF] .

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 Eq. ( 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 [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] .

In our previous work, we expanded the Random Walk with Restart algorithm to Multiplex (RWR-M) and Multiplex-Heterogeneous networks (RWR-MH) [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] . 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 [START_REF] Kivelä | Multilayer networks[END_REF][START_REF] De Domenico | Navigability of interconnected networks under random failures[END_REF] . 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 [START_REF] Battiston | Structural measures for multiplex networks[END_REF] . We do not take into account potential self-interactions and there- fore set A [α] (i, i) = 0 ∀ i = 1, . . . , n . In addition, we consider that v α i represents the node i in layer α. Thus, we can represent a multiplex graph by its adjacency matrix:

(3) [START_REF] Hamilton | Methods and applications[END_REF] , . . . , A [L] 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 β = α 38 . 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 Eq. ( 6) on multiplex graphs [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] .

v∈V M KL(sim G (v, .) � sim Emb (v, .)) (4) L = - v∈V M sim G (v, .) log(sim Emb (v, .)) (5)
L NCE = u ∼ P v ∼ sim G (u, .) logP W (D = 1 | sim Emb (u, v)) + s.E v neg ∼Q(u) logP W (D = 0 | sim Emb (u, v)) (6) p T t+1 = Mp T t (7) p T t+1 = (1 -r)Mp T t + rp T 0 (8) A = A
In this study, we independently run the RWR-M algorithm n times, using each time a different node as seed. As a result, we obtain a n × n matrix in which each column describes the probability of finding the particle in every network node when the steady state is reached. We use this probability distribution as a measure of similarity between a given node and all the other nodes of the multiplex graph. Hence, we have

u∈V M sim G (v, u) = 1 ∀ v ∈ V M ,
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) [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] .

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) [START_REF] Lee | Ranking nodes on a heterogeneous graph for flexible hybrid recommender systems[END_REF] . In our previous study [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] , 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 (Fig. 1C), based on the work of [START_REF] Dursun | A tool for gene prioritization using complete multiplex heterogeneous networks[END_REF] . Let us consider a L-layers multiplex graph,

G M = (V M , E M ) , with n × L nodes, V M = v α i , i = 1, . . . , n, α = 1, . . . , L .
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 G B = (V M ∪ U M , E B ) with E B ⊆ V M × U M .
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 = (V M ∪ U M , E [α] B )
to define the multiplex-heterogeneous graph. We can then describe a multiplex-heterogeneous graph, G MH = (V MH , E MH ) , 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 Fig. 1C). We accordingly defined a multiplex-heterogeneous transition matrix and expanded the restart probability vector. This gave us the opportunity to extent and apply Eq. ( 6) on multiplex-heterogeneous graphs [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF][START_REF] Dursun | A tool for gene prioritization using complete multiplex heterogeneous networks[END_REF] .

In the context of MultiVERSE, we independently run the RWR-MH algorithm n + m times. In each execu- tion, 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 similarity matrix matching VERSE input criteria, i.e

u∈V MH sim G (v, u) = 1 ∀ v ∈ V EM .
We set the RWR-MH parameters to the same values used in the original study ( r = 0.7

, τ = (1/L, 1/L, . . . , 1/L) , δ = 0.5 , = 0.5 , η = 0.5) 8 .

MultiVERSE algorithm. Algorithm 1 presents the pseudo-code of MultiVERSE based on RWR on multi-

plex and multiplex-heterogeneous networks [START_REF] Valdeolivas | Random walk with restart on multiplex and heterogeneous biological networks[END_REF] and Kullback-Leibler optimization from the VERSE algorithm [START_REF] Tsitsulin | Versatile graph embeddings from similarity measures[END_REF] .

Our implementation of VERSE with NCE is slightly different from the original. We perform first the RWR-M or RWR-MH for all the nodes of the network 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 truncate 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 and reduce memory constraints by filtering out the lowest probabilities and by reducing the size of the similarity matrix. 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 networks with more than 5000 nodes. For smaller networks, we set this parameter to 10% -20% of the number of nodes of the network, depending on the shape of the distribution. These choices of N max have been done for memory and quality of embeddings reasons. Indeed, these values are sufficient to obtain high quality embeddings, and avoid storing and manipulating the whole output of RWR-M(H), which is a n × n matrix. We store with this truncated sampling strategy a n × N m ax matrix. These two steps (lines 6 and 7) were not in the

V M = v α i , i = 1, . . . , n, α = 1, . . . , L , E M = (v α i , v α j ), i, j = 1, . . . , n, α = 1, . . . , L, A [α] (i, j) � = 0 (v α i , v β i ), i = 1, . . . , n, α � = β . V MH ={V M ∪ U M } E MH = ∪ α=1,...,L E [α] B ∪ E V M ∪ E U M
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 parameter s has been set to s=10 for networks with a number of nodes superior to 5000 and to s = 3 as in the VERSE original algorithm for smaller networks. The precision of the NCE depends on the parameter s, and small values work well in practice [START_REF] Mnih | A fast and simple algorithm for training neural probabilistic language models[END_REF] . The update can also be seen as the training part with lr as the learning rate of the binary classifier of the NCE estimation as described in Eq. ( 5). The whole process is repeated until the maximum steps are reached.

Regarding computational time, it depends mainly on the available number of cores and number of nodes ((as RWR has a time complexity of O(n 2 ))=. On a i7-6820HQ CPU @2.70GHz with 8 cores and 48 Gb of RAM, the whole computation of MultiVERSE for the molecular multiplex network (see next section) with d = 128 and s = 10 takes 45 minutes.

MultiVERSE is freely available on github at https:// github. com/ Lpiol/ Multi VERSE.

Evaluation protocol

We propose a benchmark to compare the performance of MultiVERSE and other embedding methods for multiplex and multiplex-heterogeneous networks. The performances are evaluated through 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 evaluations.

Multiplex network datasets. We used 7 multiplex networks (2 molecular, 1 disease, 1 neuronal, 1 co-authorship and 2 social networks) 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:// comun elab. fbk. eu/ data. php. We constructed the other two networks, DIS and MOL. A description of each of these multiplex networks follows. The number of nodes and edges of the different layers are detailed in Table 1.

• CKM physician innovation (CKM) is a multiplex network describing how physicians in four towns in Illinois used the new drug tetracycline [START_REF] Coleman | The diffusion of an innovation among physicians[END_REF] . 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) is a multiplex social network composed of 3 layers based on co-working, friend- ship and advice between partners and associates of a corporate law partnership [START_REF] Emmanuel | The collegial phenomenon. the social mechanisms of cooperation among peers in a corporate law partnership[END_REF] .

• Caenorabidis Elegans connectome (C.ELE) is a neuronal multiplex network composed of 3 layers correspond- ing to different synaptic junctions [START_REF] Chen | Wiring optimization can relate neuronal structure and function[END_REF][START_REF] De Domenico | A tool for multilayer analysis and visualization of networks[END_REF] : electrical, chemical poladic and chemical monadic. • ArXiv network (ARXIV) 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 [START_REF] De Domenico | Identifying modular flows on multilayer networks reveals highly overlapping organization in interconnected systems[END_REF] . 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) containing more than 1000 edges. • Homo sapiens network (HOMO) is composed of 4 layers extracted from the original network on CoMuNe Lab [START_REF] De Domenico | A tool for multilayer analysis and visualization of networks[END_REF] , keeping physical association, direct interaction, association and co-localization layers. The data are initially extracted from BioGRID 46 • Disease multiplex network (DIS) has been constructed, composed of 3 layers: i) A disease-disease network based on a projection of a disease-drug network from the Comparative Toxicogenomics Database (CTD) [START_REF] Davis | The comparative toxicogenomics database: Update 2019[END_REF] extracted from BioSNAP [START_REF] Zitnik | BioSNAP Datasets: Stanford biomedical network dataset collection[END_REF] . 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 [START_REF] Hagberg | Exploring network structure, dynamics, and function using networkx[END_REF] . 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 [START_REF] Zhou | Human symptoms-disease network[END_REF] . Similarly to [START_REF] Zhou | Human symptoms-disease network[END_REF] , 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 [START_REF] Jensen | Temporal disease trajectories condensed from population-wide registry data covering 6.2 million patients[END_REF] . • Human molecular multiplex network (MOL) is a molecular network, consisting of 3 layers: (i) A protein- protein interaction (PPI) layer corresponding to the fusion of 3 datasets: APID (apid.dep.usal.es) (Level 2, human only), Hi-Union and Lit-BM (http:// www. inter actome-atlas. org/ downl oad). (ii) A pathways layer extracted from NDEx [START_REF] Pratt | the network data exchange[END_REF] and corresponding to the human Reactome data [START_REF] Croft | The reactome pathway knowledgebase[END_REF] . iii) A molecular complexes layer constructed from the fusion of Hu.map [START_REF] Drew | Integration of over 9,000 mass spectrometry experiments builds a global map of human protein complexes[END_REF] and Corum [START_REF] Giurgiu | The comprehensive resource of mammalian protein complexes-2019[END_REF] , using OmniPathR [START_REF] Turei | Integrated intra-and intercellular signaling knowledge for multicelullar omics analysis[END_REF] . Methods implemented for comparisons. 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 (only in the link prediction task).

Monoplex network embedding methods.

• deepwalk 5 : This method is based on non-biased random walks, and apply the skip-gram algorithm [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF] 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 6 : This method is an extension of deepwalk with a pair of parameters p and q that biases 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 [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] . We set the other parameters as for deepwalk. • LINE [START_REF] Tang | Line: Large-scale information network embedding[END_REF] : LINE is not based on random walks, but computes the similarities using an adjacency-based proxim- ity 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 nodes in the network (only neighbours), and second order proximity look for nodes sharing many connections.

We set the negative ratio to 5.

Multiplex network embedding methods.

• OhmNet 13 : 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 with first layer as the higher in the hierarchy and the other layers are defined at the second level of the hierarchy, in the same way as [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF] . • MNE 12 : 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 (https:// github. com/ HKUST-KnowC omp/ MNE). • Multi-node2vec 14 : 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 [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] . Table 2 provides formal definitions of these heuristic scores.

Evaluation tasks. On multiplex networks, we evaluate the different methods 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 [START_REF] Perozzi | Online learning of social representations[END_REF][START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF][START_REF] Zitnik | Predicting multicellular function through multi-layer tissue networks[END_REF] for fair comparisons, and used the package EvalNE v0.3.1 [START_REF] Mara | A framework for evaluating network embeddings on link prediction[END_REF] . EvalNE is a package dedicated to the evaluation of network embedding.

From node embeddings to edges. 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 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 network 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 links removed from the original multiplex networks (Fig. 3). We remove 30% of the links in each layer of the original networks. We applied the Andrei Broder algorithm 58 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 Fig. 3):

• The methods specifically designed for monoplex network embedding (node2vec, deepwalk and LINE) are applied individually on each layer of the multiplex 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. • Methods specifically designed 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 initially. 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 [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF][START_REF] Zhang | Scalable multiplex network embedding[END_REF] . 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. Table 2. Definition of the heuristic scores of a link (u, v) in the graph G(V, E). N (u) denotes the set of neighbour nodes of node u ∈ V in G(V, E).

Score Definition

Jaccard Coefficient (JC)

|N (u)∩N (v)| |N (u)∪N (v)| Common neighbours (CN) |N (u) ∩ N (v)| Adamic Adar (AA) t∈|N (u)∩N (v)| 1 log |N (t)| Preferential attachment (PA) |N (u)|.|N (v)|
Table 3. Embedding operators used to predict edges in the tasks of link prediction and network reconstruction. The definitions describe the i th components of g(u, v).

Operators Symbol Definition

Hadamard

[f (u) f (v)] i = fi (u) * fi (v) 2 Average ⊞ [f (u) ⊞ f (v)] i = f i (u) + f i (v) Weighted-L1 . 1 � f i (u).f i (v) � 1 i =| f i (u) -f i (v) | Weighted-L2 . 2 � f i (u).f i (v) � 2 i =| f i (u) -f i (v) | 2 Cosine cos cos[f (u), f (v)] i = fi (u) * fi (v) �fi (u)��fi (v)� Figure 3.
General approach for link prediction on multiplex networks: (top) for the link prediction heuristics, we apply them to each layer and average them across all layers; (center) for monoplex-based methods, we embed each layer with the given method, then average it; (bottom) for multiplex-based methods, we apply the specific embedding method to the network. The embedding operators are then applied to monoplex-and multiplexbased method embeddings. The three types of methods are finally evaluated for link prediction using a binary classifier and a ROC-AUC is computed. Network reconstruction. Network reconstruction is another approach to evaluate network embedding methods [START_REF] Bagavathi | Multi-net: A scalable multiplex network embedding framework[END_REF][START_REF] Wang | Structural deep network embedding[END_REF][START_REF] Goyal | Graph embedding techniques, applications, and performance: A survey[END_REF] . 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 network.

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) 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 [START_REF] Fernández | Learning From Imbalanced Data Sets[END_REF] . In order to account for class imbalance, we used the precision@K 59 . 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 (Fig. 4). We 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 [START_REF] Wang | Structural deep network embedding[END_REF][START_REF] Goyal | Graph embedding techniques, applications, and performance: A survey[END_REF] . 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 [START_REF] More | Review of random forest classification techniques to resolve data imbalance[END_REF] . In network reconstruction, the results correspond to the training phase of the classifier, there is no test phase.

Evaluation of multiplex-heterogeneous network embedding. 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 [START_REF] Piñero | The disgenet knowledge platform for disease genomics: 2019 update[END_REF] 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 Fig. 1C. • Drug-target multiplex-heterogeneous network We use the molecular multiplex network (MOL) from the pre- vious 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/ bione tdata/ 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 [START_REF] Cheng | Network-based prediction of drug combinations[END_REF] and the edges are based on drug combinations as reported in clinical data, (iii) the third layer (13397 edges, 658 nodes) is the adverse drugdrug interactions network available in [START_REF] Cheng | Network-based prediction of drug combinations[END_REF] . The drug-target bipartite network has been extracted from the same publication [START_REF] Cheng | Network-based prediction of drug combinations[END_REF] , 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 multiplex-heterogeneous networks described previously, the idea behind this evaluation is to test if we can predict gene-disease and drug-gene links. Comparisons with other approaches are not possible as, to our knowledge, no existing multiplex-heterogeneous network embedding method are currently available in the literature. www.nature.com/scientificreports/ Case study: discovery of new gene-disease associations. Link prediction. Our aim for this casestudy is to predict new gene-disease links. We thereby applied MultiVERSE on the full gene-disease multiplexheterogeneous 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 disease nodes. Finally, we select the top 5 new gene-disease associations for each disease.

Clustering. We also applied MultiVERSE to the gene-disease multiplex-heterogeneous gene-disease network, followed by spherical K-means [START_REF] Buchta | Spherical k-means clustering[END_REF] to cluster the vector representations of nodes. Spherical K-means clustering is well-adapted to high-dimensional clustering [START_REF] Zhong | Efficient online spherical k-means clustering[END_REF] . We define the number of clusters for spherical K-means to 500, in order to obtain cluster sizes that can be analysed from a biological point of view.

Results

Evaluation results for multiplex network embeddings. 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 the heuristic 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 [START_REF] Zhang | Scalable multiplex network embedding[END_REF] . 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 deepwalkav perform very well and even outperform multiplex-based approaches on various scenarios, for instance link prediction 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, 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 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. 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 gene-disease and drug-target links from the corresponding multiplex-heterogeneous networks.

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 Mul-tiVERSE combined with these operators for Hutchinson-Gilford Progeria Syndrome (HGPS) and Xeroderma pigmentosum VII (Table 7) .

Hutchinson-Gilford Progeria Syndrome. Hutchinson-Gilford Progeria Syndrome (HGPS) is a rare premature aging genetic disease characterized by postnatal growth retardation, midface hypoplasia, micrognathia, premature atherosclerosis, coronary artery disease, lipodystrophy, alopecia and generalized osteodysplasia [START_REF] Sandre-Giovannoli | Lamin a truncation in Hutchinson-Gilford progeria[END_REF] . HGPS is caused by mutations in the LMNA genes that cause the production of a toxic form of the Lamin A protein called Progerin.

MultiVERSE top predictions reveal interesting candidate genes (Table 7). In particular, NOS2 encodes a nitric oxide synthase expressed in liver. It has been associated with longevity [START_REF] Montesanto | Common polymorphisms in nitric oxide synthase (nos) genes influence quality of aging and longevity in humans[END_REF] . 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 [START_REF] Davizon-Castillo | Tnf-α-driven inflammation and mitochondrial dysfunction define the platelet hyperreactivity of aging[END_REF] and has been previously linked to Progeria [START_REF] Osorio | Nuclear lamina defects cause atm-dependent nf-κ b activation and link accelerated aging to a systemic inflammatory response[END_REF] . TERF1 and TERF2 both encode telomerebinding proteins and TERF2IP encodes a protein that is part of a complex involved in telomere length regulation. HGPS patients show increased activation of DNA damage signalling at telomeres associated to reduced telomere length [START_REF] Decker | Telomere length in Hutchinson-Gilford progeria syndrome[END_REF] . In addition, it has been reported DNA damage accumulation and TRF2 degradation in atypical Werner syndrome (adult Progeria) fibroblasts with LMNA mutations [START_REF] Saha | Dna damage accumulation and trf2 degradation in atypical werner syndrome fibroblasts with lmna mutations[END_REF] . POT1 also produces a telomeric protein that has been linked to the Werner syndrome [START_REF] Sowd | Mechanism and substrate specificity of telomeric protein pot1 stimulation of the werner syndrome helicase[END_REF] , the maintenance of haematopoeitic stem cell activity during aging [START_REF] Hosokawa | The telomere binding protein pot1 maintains haematopoietic stem cell activity with age[END_REF] and cellular senescence [START_REF] Li | Seryl trna synthetase cooperates with pot1 to regulate telomere length and cellular senescence[END_REF] . IL6 encodes a cytokine involved in inflammation, which have also been linked to aging [START_REF] Maggio | Interleukin-6 in aging and chronic disease: A magnificent pathway[END_REF] . 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 multiplex 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 in parenthesis under the name of the network. In the case of large networks (DIS, ARXIV, HOMO and MOL) MultiVERSE achieves the best performance in combination with different operators. Finally,SOD1 and SOD2 are members of the superoxide dismutase multigene family that destroy free superoxide radicals. They both have been associated to aging and cellular senescence [START_REF] Zhang | A new role for oxidative stress in aging: The accelerated aging phenotype in sod1-/-mice is correlated to increased cellular senescence[END_REF][START_REF] Velarde | Mitochondrial oxidative stress caused by sod2 deficiency promotes cellular senescence and aging phenotypes in the skin[END_REF] .

Xeroderma pigmentosum VII. Xeroderma Pigmentosum (XP) is characterized by extreme sensitivity to sunlight, resulting in sunburns, pigment changes in the skin and a highly elevated incidence of melanoma. It is a genetically heterogeneous autosomal recessive disorder. Several XP 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 [START_REF] Kraemer | Xeroderma pigmentosum: Cutaneous, ocular, and neurologic abnormalities in 830 published cases[END_REF][START_REF] Vermeulen | Xeroderma pigmentosum complementation group g associated with cockayne syndrome[END_REF] . MultiVERSE identified various candidates for this disease (see Table 7), we detail here the most interesting ones. TNF has been related to XP [START_REF] Capulas | Ultraviolet-b-induced apoptosis and cytokine release in xeroderma pigmentosum keratinocytes[END_REF] and skin tumour development [START_REF] Arnott | Expression of both tnf-α receptor subtypes is essential for optimal skin tumour development[END_REF] . IL6 is involved in melanoma, one of the major phenotypes of XP [START_REF] Lu | Interleukin 6: a fibroblast-derived growth inhibitor of human melanoma cells from early but not advanced stages of tumor progression[END_REF] . SOD2, also predicted as candidate for HGPS, have been recently associated to melanoma [START_REF] Yuan | Braf mutant melanoma adjusts to braf/mek inhibitors via dependence on increased antioxidant sod2 and increased reactive oxygen species levels[END_REF] . TP53 is a tumor suppressor implicated in many cancers, in particular melanomas [START_REF] Giglia-Mari | Tp53 mutations in human skin cancers[END_REF] . It has also been associated to XP [START_REF] Sarasin | Familial predisposition to tp53/complex karyotype mds and leukemia in dna repair-deficient xeroderma pigmentosum[END_REF] . VCAM1 encodes the Vascular Cell Adhesion Molecule, associated to melanoma [START_REF] Klemke | High affinity interaction of integrin α4β 1 (vla-4) and vascular cell adhesion molecule 1 (vcam-1) enhances migration of human melanoma cells across activated endothelial cell layers[END_REF] , and NUP62 encodes the Nuclear pore glycoprotein p62, involved in cell carcinoma proliferation [START_REF] Hazawa | Rock-dependent phosphorylation of nup 62 regulates p63 nuclear transport and squamous cell carcinoma proliferation[END_REF] . ERCC2 produces the XPD protein, mutated in XP group D [START_REF] Taylor | Xeroderma pigmentosum and trichothiodystrophy are associated with different mutations in the xpd (ercc2) repair/transcription gene[END_REF] .

Discovery of new gene-disease associations with clustering. Another illustration of the advantages of multiplexheterogeneous network embedding is clustering. We identify clusters with spectral K-means, 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. In addition, it has been shown that clustering from embeddings outperforms the other methods for the detection of biological communities [START_REF] Nelson | To embed or not: Network embedding as a paradigm in computational biology[END_REF] .

Cluster containing the HGPS disease node. The cluster containing the HGPS disease node (see Fig. 5) contains the LMNA node. LMNA mutations have been observed in many diseases that also belong to the identified cluster, including the Heart-hand syndrome (Slovenian type) [START_REF] Renou | Heart-hand syndrome of Slovenian type: A new kind of laminopathy[END_REF] , lipodystrophy associated with mandibuloacral dysplasia [START_REF] Agarwal | Severe mandibuloacral dysplasia-associated lipodystrophy and progeria in a young girl with a novel homozygous arg527cys lmna mutation[END_REF] , the Charcot-Marie-Tooth disease, type 2B1 [START_REF] Sinha | Progeria: A rare genetic premature ageing disorder[END_REF] , LMNA-related muscular distrophy, and different cardiac diseases caused by LMNA mutations [START_REF] Van Tintelen | High yield of LMNA mutations in patients with dilated cardiomyopathy and/or conduction disease referred to cardiogenetics outpatient clinics[END_REF] . We also analysed the cluster's genes annotations with g:Profiler [START_REF] Raudvere | Profiler: A web server for functional enrichment analysis and conversions of gene lists (2019 update)[END_REF] (default parameters). We found significant enrichments in several annotations related to cell nuclear organization. One of the most significant enrichments is nuclear envelope, involving the following genes: EMD, LEMD2, LMNA, KPNA1, MLIP, TMEM43, and ZMP-STE24. HGPS is a disorder of the nuclear envelope [START_REF] Worman | Diseases of the nuclear envelope[END_REF] .

Cluster containing the Xeroderma pigmentosum VII disease node. We also analysed the cluster containing the Xeroderma pigmentosum VII disease node (Fig. 6). The cluster contains different diseases, including Xeroderma pigmentosum with normal DNA repair rates, and Cerebro-oculo-facio-skeletal syndrome 4, which is also a nuclear-excision repair disorder [START_REF] Graham | Cerebro-oculo-facio-skeletal syndrome with a nucleotide excision-repair defect and a mutated xpd gene, with prenatal diagnosis in a triplet pregnancy[END_REF] .

Several genes known for their implication in XP are present in the cluster, such as ERCC1, ERCC4, ERCC5, XPA and XPC [START_REF] Sugasawa | Xeroderma pigmentosum genes: functions inside and outside DNA repair[END_REF] . Using the complete list of genes in the cluster as an input for g:Profiler [START_REF] Raudvere | Profiler: A web server for functional enrichment analysis and conversions of gene lists (2019 update)[END_REF] (default parameters), we identified several significantly enriched annotations. Among them, we can cite nucleotide-excision DNA repair, defective DNA repair after ultraviolet radiation damage or response to ultraviolet radiation. XP patients show important impairments in these biological processes [START_REF] Kraemer | Xeroderma pigmentosum: Cutaneous, ocular, and neurologic abnormalities in 830 published cases[END_REF] . SPRTN is another gene of interest. It encodes a metalloprotease that repairs DNA-protein crosslinks. SPRTN does not share interactions with genes known to be mutated in XP, but has been shown to be involved in UV sensibilization and cancer [START_REF] Hiom | Sprtn is a new player in an old story[END_REF] .

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-of-the-art methods using link prediction and network reconstruction. We show that MultiVERSE outperforms various methods specifically developed for multiplex network embedding. Our results suggest that there are several advantages to use RWR for the computation of the similarity between nodes: we experimentally demonstrate that methods based on the SkipGram approach with the truncated random walks (such as node2vec-av, deepwalk-av, ohmnet, multi-node2vec and MNE) are less effective to learn node embedding from multiplex networks than MultiVERSE in several contexts. In particular, multi-node2vec, which is an extension of node2vec (and therefore SkipGram) to multiplex networks, shows reduced performance for our evaluation tasks than MultiVERSE. In addition, node2vec-av is also less effective to learn node embeddings from multiplex networks. SkipGram uses a truncated random walk whereas MultiVERSE uses a (non-truncated) global random walk with RWR to compute the similarity. 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.

The methods we found in the literature with available code for multiplex network embedding (i.e., ohmnet, multi-node2vec, MNE) are all based on node2vec. Recently, several graph neural networks approaches have been developed for network embedding, including GCN and GAE/VGAE [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF][START_REF] Kipf | Variational graph auto-encoders[END_REF] . These methods allow learning high quality embeddings for monoplex networks [START_REF] Khosla | A comparative study for unsupervised network representation learning[END_REF] . However, to the best of our knowledge, they have not yet been extended to multiplex or multiplex-heterogeneous network embedding. It will be interesting to compare MultiVERSE with this class of graph neural networks approaches once they will be developed for multiplex and multiplex-heterogeneous network embedding.

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 of different types of nodes. We demonstrate its effectiveness for link prediction and illustrate its usefulness for the study of genedisease 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 multiplexheterogeneous 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 100 . In social science, multiplex networks are gaining interest to understand human behaviour 101 . Multiplex-heterogeneous network embedding could give insights on epidemic spread 102,103 , socio-economic systems 104 or socio-ecological systems 105 .

Figure 1 .

 1 Figure 1. Illustrations of the different types of networks. (A) A multiplex network. The different layers share the same set of nodes but different types of edges. (B) A heterogeneous network. The two networks are composed of different types of nodes and edges, connected by bipartite interactions (black dashed lines). (C) A multiplex-heterogeneous network composed of two multiplex networks. The multiplex networks are connected by bipartite interactions (dashed lines). For the sake of simplicity, the figure does not represent all the possible bipartite interactions (each layer of a given multiplex is in reality linked with every layer of the other multiplex).

Figure 2 .

 2 Figure 2. Overview of the MultiVERSE pipeline. Starting from a multiplex-heterogeneous network, we represent its structure through an adjacency matrix (size |V | × |V | ); we then compute a similarity matrix using Random Walk with Restart algorithm, and apply an optimized version of the VERSE algorithm to compute the embeddings. The resulting matrix of embeddings will be used for the applications.

Figure 4 .

 4 Figure 4.General approach for network reconstruction on multiplex networks: (top) for monoplex-based methods, embed each layer with the given method, then average it; (bottom) for multiplex-based methods, apply the specific embedding method to the network. Embedding operators are then applied to monoplex-and multiplex-based method embeddings. The three types of methods are finally evalutated for network embedding using a binary classifier and a precision@K score is computed.

Table 7 .Figure 5 .

 75 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 .

 6 Figure 6. Cluster containing the Xeroderma pigmentosum VII disease node. 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 1 .

 1 Description of the 7 multiplex networks used for the evaluation protocol.

	Dataset	Layers Nodes	Edges
		1	215	449
	CKM	2	231	498
		3	228	423
		1	71	717
	LAZEGA	2	69	399
		3	71	726
		1	253	514
	C.ELE	2	260	888
		3	278	1703
		1	1558	3013
		2	5058	14,387
		3	2826	6074
	ARXIV	4 5	1572 3328	4423 7308
		6	1866	4420
		7	1246	1947
		8	4614	11,517
		1	12,,345 48,528
	HOMO	2 3	14,770 1626	83,414 1953
		4	5680	18,381
		1	3891	117,527
	DIS	2	4155	101,104
		3	434	3137
		1	14,704	122,211
	MOL	2	7926	194,500
		3	8537	63,561

Table 4 .

 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 multiplex network, the best score is in bold; for each operator, the best scores are underlined. Overall, the MultiVERSE algorithm combined with the Hadamard shows the best scores.

	Operators	Method	CKM	LAZEGA C.ELE ARXIV DIS	HOMO MOL
		CN-av	0.4944 0.6122	0.5548 0.5089	0.5097 0.5113	0.5408
	Link prediction heuristics	AA-av JC-av	0.4972 0.6105 0.4911 0.523	0.549 0.5424 0.5113 0.5081	0.5428 0.5112 0.5425 0.5113	0.5404 0.5433
		PA-av	0.5474 0.6794	0.5634 0.5139	0.496	0.5185	0.5278
		node2vec-av	0.7908 0.6372	0.8552 0.9775	0.9093 0.8638	0.8753
		deepwalk-av	0.7467 0.6301	0.8574 0.9776	0.9107 0.8638	0.8763
		LINE-av	0.5073 0.4986	0.5447 0.8525	0.9013 0.8852	0.8918
	Hadamard	Ohmnet	0.7465 0.7981	0.833	0.9605	0.9333 0.9055	0.8613
		MNE	0.5756 0.6356	0.794	0.9439	0.9099 0.8313	0.8736
		Multi-node2vec 0.8182 0.7884	0.8375 0.9581	0.8528 0.8592	0.8835
		MultiVERSE	0.8177 0.8269	0.8866 0.9937	0.9401 0.917	0.9259
		node2vec-av	0.7532 0.737	0.8673 0.9738	0.885	0.6984	0.7976
		deepwalk-av	0.7226 0.7094	0.8635 0.9751	0.8888 0.7142	0.8089
		LINE-av	0.6091 0.5776	0.6192 0.7539	0.8586 0.7439	0.7792
	Weighted-L1	Ohmnet	0.7421 0.7849	0.8128 0.8488	0.8503 0.7007	0.6983
		MNE	0.6289 0.6523	0.8019 0.7805	0.8313 0.7619	0.8182
		Multi-node2vec 0.8611 0.8089	0.8261 0.9659	0.8628 0.8472	0.8997
		MultiVERSE	0.7043 0.7789	0.7516 0.8647	0.7754 0.683	0.7273
		node2vec-av	0.7556 0.6851	0.8691 0.9743	0.8867 0.7048	0.8028
		deepwalk-av	0.7221 0.6904	0.864	0.9771	0.8891 0.7145	0.813
		LINE-av	0.5851 0.5756	0.6275 0.7609	0.8621 0.7429	0.7835
	Weighted-L2	Ohmnet	0.7505 0.7788	0.8166 0.8439	0.8599 0.7041	0.6992
		MNE	0.601	0.5397	0.7999 0.7815	0.8333 0.7483	0.8122
		Multi-node2vec 0.8637 0.8091	0.8282 0.968	0.8675 0.8525	0.9004
		MultiVERSE	0.7125 0.7801	0.7441 0.8661	0.7808 0.6918	0.7475
		node2vec-av	0.59	0.6596	0.6842 0.6615	0.8256 0.8308	0.777
		deepwalk-av	0.5954 0.657	0.6784 0.6582	0.8267 0.8307	0.7737
		LINE-av	0.5465 0.6581	0.6699 0.6465	0.8477 0.8653	0.8276
	Average	Ohmnet	0.5764 0.656	0.7334 0.6772	0.8533 0.8825	0.7962
		MNE	0.5882 0.6615	0.7028 0.6723	0.8242 0.8024	0.783
		Multi-node2vec 0.5571 0.6584	0.7365 0.6657	0.8222 0.8216	0.7589
		MultiVERSE	0.5963 0.6728	0.7438 0.6752	0.8586 0.8643	0.812
		node2vec-av	0.7805 0.7335	0.8515 0.9711	0.8643 0.7368	0.8105
		deepwalk-av	0.7465 0.7066	0.8416 0.9724	0.8667 0.7512	0.8079
		LINE-av	0.545	0.5126	0.5477 0.8198	0.7409 0.6745	0.816
	Cosine	Ohmnet	0.7898 0.7352	0.8094 0.9642	0.859	0.7829	0.7909
		MNE	0.6203 0.6506	0.7877 0.8951	0.8347 0.6474	0.8102
		Multi-node2vec 0.8532 0.7931	0.7815 0.9435	0.7151 0.8477	0.8884
		MultiVERSE	0.8148 0.8171	0.8719 0.9909	0.8775 0.8776	0.9103

Table 6 .

 6 ROC

				LAZEGA				HOMO	MOL
	Operators	Method	CKM (95%)	(95%)	C.ELE (95%) ARXIV (5%) DIS (2,5%)	(2,5%)	(2.5%)
		node2vec-av 0.6764	0.9174	0.4526	0.8207	0.5578	0.7599	0.2989
		deepwalk-av 0.6564	0.9351	0.4416	0.7886	0.5486	0.7636	0.3164
		LINE-av	1.0	0.9924	0.8924	0.8204	0.4955	0.5191	0.4006
	Hadamard	Ohmnet	0.7842	0.8334	0.5329	0.9156	0.4811	0.6979	0.2591
		MNE	0.9505	0.9094	0.2728	0.7891	0.4218	0.3641	0.1316
		Multi-node-2vec	0.8352	0.8811	0.6875	0.8605	0.6063	0.7584	0.3123
		MultiVERSE 0.9687	0.9695	0.7436	0.9015	0.6734	0.8729	0.3674
		node2vec-av 0.5923	0.9494	0.5129	0.6922	0.5859	0.8123	0.3194
		deepwalk-av 0.5791	0.9784	0.4896	0.6878	0.5921	0.7984	0.3206
		LINE-av	0.9985	0.9953	0.9229	0.7837	0.4921	0.6839	0.3586
	Weighted-L1	Ohmnet	0.7355	0.8581	0.5785	0.8771	0.6025	0.8019	0.3769
		MNE	0.9926	0.975	0.4722	0.8593	0.4377	0.5241	0.1861
		Multi-node-2vec	0.8636	0.9235	0.7379	0.7684	0.6356	0.7649	0.2671
		MultiVERSE 0.8545	0.9638	0.7444	0.8705	0.6678	0.7913	0.3559
		node2vec-av 0.5886	0.9436	0.5097	0.6983	0.5953	0.8193	0.352
		deepwalk-av 0.5829	0.9672	0.5146	0.6877	0.5857	0.805	0.3233
		LINE-av	1.0	0.9962	0.9367	0.7749	0.4945	0.6697	0.392
	Weighted-L2	Ohmnet	0.7418	0.8687	0.5724	0.8694	0.6209	0.8143	0.3701
		MNE	0.9926	0.9764	0.4646	0.8818	0.4351	0.5529	0.176
		Multi-node-2vec	0.8644	0.93	0.7568	0.7548	0.6361	0.7896	0.2922
		MultiVERSE 0.8653	0.969	0.754	0.8776	0.6784	0.7876	0.3701
		node2vec-av 0.8408	0.917	0.4817	0.889	0.5587	0.6809	0.2686
		deepwalk-av 0.8331	0.9379	0.501	0.8853	0.5318	0.6714	0.2795
		LINE-av	0.9855	0.9382	0.7103	0.8725	0.5093	0.5677	0.3244
	Average	Ohmnet	0.9412	0.8287	0.5825	0.906	0.4989	0.6551	0.2887
		MNE	0.9179	0.9151	0.2966	0.7146	0.4175	0.352	0.1444
		Multi-node-2vec	0.9767	0.8937	0.6726	0.9498	0.6243	0.6216	0.2901
		MultiVERSE 0.978	0.9059	0.5326	0.9758	0.6316	0.7204	0.4143
		node2vec-av 0.5103	0.4936	0.18	0.2537	0.1825	0.116	0.0441
		deepwalk-av 0.4807	0.4776	0.1741	0.2835	0.1854	0.1036	0.0462
		LINE-av	0.3291	0.4974	0.1867	0.2638	0.2384	0.1476	0.0454
	Cosine	Ohmnet	0.5696	0.509	0.1718	0.2655	0.1984	0.1311	0.044
		MNE	0.3169	0.4536	0.1768	0.2445	0.1957	0.1667	0.044
		Multi-node-2vec	0.5127	0.52	0.186	0.273	0.195	0.1032	0.0461
		MultiVERSE 0.6395	0.5026	0.1818	0.254	0.1983	0.1522	0.0474
	Operators		Gene-disease bipartite		Drug-target bipartite
	Hadamard		0.95			0.9701
	Weighted-L1		0.7962			0.8057
	Weighted-L2		0.7951			0.8055
	Average		0.9603			0.9703
	Cosine		0.7765			0.8338

-AUC scores for link prediction using MultiVERSE on 2 multiplex-heterogeneous 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.

Scientific Reports | (2021) 11:8794 | https://doi.org/10.1038/s41598-021-87987-1

Scientific Reports | (2021) 11:8794 | https://doi.org/10.1038/s41598-021-87987-1

Lenormand, M. et al. Multiscale socio-ecological networks in the age of information. PloS ONE 13, e0206672 (2018).

Acknowledgements

We thank Pr. Alfonso Valencia and the Computational Biology Life Sciences Group at Barcelona Supercomputing Center for the warm welcome and discussions and the use of the MareNostrum supercomputer.

Funding

L.P-L is the recipient of a Short Term Collaboration Grant for HPC 2019 from the Eurolab4HPC consortium. This project has received funding from the Excellence Initiative of Aix-Marseille University-A*Midex, a French 'Investissements d' Avenir' program.

Author contributions

Conceptualization: A.B., L.P-L, A.V.; Methodology: L.P-L., A.V.; Software: L.P-L., L.T., A.V.; Formal analysis and investigation: L.P-L.; Writing-Original draft: L.P-L.; Review and editing: A.B., L.P-L., E.R., L.T., A.V.; Visualization: A.B., L.P-L., L.T.; Supervision: A.B.; Project administration: A.B., E.R.; Funding acquisition: A.B., L.P-L., E.R.

Competing interests

The authors declare no competing interests.