Network science constitutes now a fundamental framework for studying complex systems and modeling the ever-growing data deluge that occurs in virtually all fields of knowledge. Over the last decade, multilayer networks contributed to revolutionize the study of systems characterized by multiple scales or levels, i.e. ,layers of interactions. These models shed a new light on interconnected systems by exhibiting unexpected topological correlations, robustness and synchronization properties, just to name a few. Historically, multilayer networks leveraged their comparative advantages over classical networks, focusing on the interactions of nodes both within and between layers. Although recent studies started to characterize layer properties per se, the layer characterization and its interplay with the classical node characterization still remain to be explored.

In this manuscript, we propose the node-layer duality, a unifying framework to study the structural properties of nodes and layers. We show that both node and layer properties provide complementary information when considering the second moment of their distribution ,i.e. , the variance. We extensively study these complementarity by deriving a stochastic rewiring model that selectively rewires links according to either node, or layer or both dimensions. Using this rewiring approach, we analytically derive and numerically validate a formal duality relationship between the node and layer dimensions. In particular, we derive and validate a closed-form of this relationship in the case of random multiplex and multilayer networks.

Based on local connectivity, we show that the complementarity of the node-layer duality clusters realworld multiplex networks coming from di erent contexts (social, biological, infrastructure...) into networks that have a spatial connotation and others that do not have one. Moreover, we provide a method to e ciently visualise the connectivity of multiplex networks in the node and layer dimensions.

In this e ort to characterize real-world systems, we focus our analysis on multilayer brain networks. The brain is perhaps the most striking example of complex systems that has benefited from the multilayer thrust. Actually the brain at the large-scale can be recorded using many neuroimaging technics (MRI,EEG,MEG,fMRI...) and its functional properties can be investigated along di erent modes (temporal, frequential). The properties of these brain networks are fundamental to uncover the normal brain functioning as well as new e ective biomarkers to prevent, track or even cure brain diseases.

In particular, cross-frequency coupling, i.e., interactions between di erent brain frequencies has been shown to be a major component of information transfer within the brain across spatiotemporal scales. Although modeling cross-frequency coupling was already suggested to better characterize diseases, multilayer brain networks that integrate cross-frequency coupling are still relatively unexplored. In parallel, the frequency-based multilayer networks, i.e., layers representing the frequency of brain activity showed promising results in characterizing Alzheimer's disease. Therefore, we apply the node-layer duality framework to characterize simultaneously the regions (nodes) and frequencies (layers) of multi-frequency brain networks including cross-frequency coupling.

Based on the correlation between local and global measures of connectivity and cognitive decline, we find that Alzheimer's disease seems to be characterized not only by the disruption of connectivity in specific brain regions, but more importantly by the aberrant coordination between frequencies. In particular, we recover the importance of the connectivity disruption in the alpha band ( -Hz) which is a well-kown feature of Alzheimer's disease.

We conclude on the opportunity of systematically adopting a dual characterization in the study of the structure of multilayer networks by exploiting its characteristics which are shared among duality relationships in general.

Introduction

List of Figures

. Multiscale brain organization. The di erent organizational aspects of the brain system are represented over a multidimensional manifold. Three type of dimensions, or levels, are illustrated here, i.e. , time, space, and topology. From the top to the bottom of the manifold, the scales of each organizational level go from micro to macro. Image credit: Thibault Rolland.

Adapted b) The quantity Θ , ( ) represent the information about the structure in layer carried by the properties of nodes in layer . The ratio of the two aforementioned properties give a measure of similarity between layers. If the ratio is one, the property has the same level of information for the structure level as the property for the layer . Finally, to obtain the similarity function the ratio is symmetrized as show in the picture. brain sites according to the used neuroimaging technique. Links are estimated by measuring the FC between the activity of brain nodes; this information is contained in a connectivity matrix. By means of filtering procedures, based on thresholds, only the most important links constitute the brain network. The topology of the brain network is quantified by di erent network measures that can be represented as numbers (e.g., the colored bars). These measures can be the input to statistical analysis to look for significant di erences between populations/conditions (e.g., red points correspond to brain network measures of diseased patients or tasks, blue points stand for healthy subjects or resting states. progression from cell level to entire brain. K is a noninteracting collection of neurons. KI corresponds to a cortical column with su cient functional connection density. KII represents a collection of excitatory and inhibitory populations. KIII is formed by the interaction of several KII sets and simulates the known dynamics of sensory areas with /f spectra (see inset b). KIV is formed by the interaction of three KIII sets that models the genesis of simple forms of intentional behaviors. Panel b) Schematic view of major components involved in thalamocortical interactions. Di erent shading patterns code for di erent zones of the system, i.e. , from micro (relay nuclei, thalamic reticular nuclei (TRN)) to macro scales (cortex). As indicated by key, all connections shown are excitatory except for the connection from the reticular cell to the relay cell, which is inhibitory. ) between the coreness disruption index ( ) and the memory impairment of AD patients as measured by the free recall (FR) test. Panel c) Boxplots show the values of coreness disruption index ( ) obtained by progressively removing the links preferentially connected to the multiplex periphery of the HC group. The blue (x-axis HC) and red (x-axis AD) boxplots illustrate respectively the values for the HC and AD groups. . Topological duality in multilayer networks. The unitary element of a multilayer network is the duplet ( , ) corresponding to the node and the layer . Because of the intrinsic node-layer isomorphism, the system can be equivalently represented by nodes plunged into layers or by layers plunged into nodes. Depending on the representation side, two complementary descriptions of the same system can be obtained. In the primal nodewise, connectivity is integrated across the layers and the topological properties of the nodes are provided. In the dual layerwise, connectivity is integrated across the nodes and the topological properties of the layers are obtained instead. In the networks, solid lines intralayer links, grey lines interlayer links, dashed lines replica links. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. Nodewise and layerwise distributions are non-constrained from each other. The MC distributions are normalized in order to have the same mean. a) Nodewise (blue) and layerwise (orange) MC distribution of the EuAir multiplex network [ ]. The normalized nodewise standard deviation is X = 0.73. The normalized layerwise standard deviation is Y = 0.45.b) EuAir multiplex network is rewired in a way that layerwise degree distribution stays fixed and the nodewise one is modified using the stochastic rewiring algorithm ( = 1, = 1, see Sec. . . ). The normalized standard deviations become X = 0.10, Y = 0.45.c) EuAir multiplex network is rewired in a way that nodewise degree distribution stays fixed and the layerwise one is modified using the stochastic rewiring algorithm ( = 0.66, = 1, see Sec. . . ).The normalized standard deviations become X = 0.73, Y = 0.15. Kernel density estimation are obtained using the seaborn Python library. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. Schematic illustration of the di erent types of link rewiring. The type of rewiring is determined by the combination of three parameters. Specifically, determines the probability to make a link displacement while keeping the layers unchanged.

determines the probability to make a link displacement while keeping the nodes unchanged.

determines the probability to make a link displacement while changing both nodes and layers. Dotted lines old position, solid lines new position. . E ect of intralayer link additions on node and layer multidegree centralities. In general, adding one intralayer link in the nodewise description increases the multidegree centrality in the layerwise description by a factor of (right side), while the node multidegree centrality only increases by a factor of one (left side). Because in multiplex networks only changes within layers are allowed, this e ect increments more by construction the variance of the layer multidegree centrality and its distance as compared to the nodewise counterpart. . . . . . .

. Invisibility zones for multidegree centrality distances. Nodewise distances ( X ) increase linearly with (x-axis) and (white diagonals) but they cannot see link displacements that keep nodes unchanged ( → 1). Layerwise distances ( Y ) increase linearly with (y-axis) and

(white diagonals) but they are blind to link displacements that keep layers unchanged ( → 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. Linear relation between layerwise and nodewisee multidegree centrality distances in random multilayer networks. Lower slopes (higher X ) are obtained for > . Higher slopes (higher Y ) are obtained for > ). Solid lines correspond to Eqs. . , . in the case of random multilayer networks with = = 200, = 0.0005, by varying the entire range of possible rewiring (color line). Scattered points correspond to synthetic random networks simulated with the same parameters (one synthetic random network di erent for each point).

. Linear relation between layerwise and nodewisee multidegree centrality distances in random multiplex networks. Lower slopes (higher X ) are obtained for > . Higher slopes (higher Y ) are obtained for > ). Solid lines correspond to Eqs. . , . in the case of random multiplex networks with = = 200, = 0.01 by varying the entire range of possible rewiring (color line). Scattered points correspond to synthetic random mutiplex networks simulated with the same parameters (one synthetic random network di erent for each point). . In multilayer networks, nodewise and layerwise distances equally reach the highest value when the number of nodes N is equal to the number of layers M (black lines). However, in multiplex networks, the maximum is reached when there are more nodes than layers. In addition, layerwise distances (orange) are by construction higher than nodewise distances (blue). Solid lines correspond to the theoretical formulas for random networks with connection density = 0.0005, rewiring ratio = 0.5, and uniform rewiring probability (see Sec. In bold, we show the frequencies whose significance is lower than the mean false discovery rate (mFDR) for independent tests ( < 0.0253). Nine frequencies in the alpha band are close to the statistical significance after correction with the Benjamini-Hochberg method. . . . . . . .

List of Tables

At the dawn of the 21 century, network science appeared as a major paradigm for studying complex systems, i.e. , systems that exhibit emergent phenomena, rooted in the interactions between their components. Epidemics, climate and opinion formation are a few examples of complex systems whose understanding is crucial for our societies.

In network science, the components of a system are modeled by nodes and their interaction by links, which together form a network. Such abstraction is now a common language of an interdisciplinary endeavor ranging from physics to sociology passing through computer science to tame complexity as well as the ongoing data deluge present in virtually all fields of knowledge.

Among many other fields, the concepts of network science have found a prolific echo in neuroscience. Actually, the brain is a complex system par excellence exhibiting a wide repertoire of emergent phenomena from the molecular to the behavioral scale which makes it amenable to network modeling. In network neuroscience, nodes can correspond to di erent brain sites, such as neurons, neuronal ensembles or even large-scale areas but also to electrical or optical sensors. The links correspond to either anatomical/structural connections or functional/dynamical interactions between nodes. This holistic approach opens up novel perspectives on the brain, e.g. , by demontrating, that brain networks are globally characterized by nodes that form an optimal balance between being clustered and being integrated in the whole network at the same time. Network properties are fundamental for quantifying and understanding normal and abnormal brains and ultimately for unveiling e ective biomarkers to prevent, track or even cure brain diseases.

Nonetheless, networks and particularly brain networks su er from a major limitation: complex systems, such as the brain, often interact through multiple types of interactions along with di erent spatiotemporal scales. For instance, in a social network, individuals may exchange information through face-to-to face interactions, phone, email or online social networks. Each of these layer of interactions has a structure that is constrained by its nature and its characteristic temporal scale. Therefore, over the last decade, multilayer networks were put forward to model these interacting layers in an integrated framework. Multilayer networks exhibit unexpected properties arising from the topological interplay between layers that impact systemwide synchronization, robustness and controllability as compared to equivalent single-layer networks. Actually, multilayer network analysis has leveraged its comparative strengths, with respect to classical networks, focusing on the interactions of nodes both within and between layers.

However, a complementary characterization could in principle be derived based on how the layers interact across nodes. Recent studies ), 'Interlayer connectivity reconstruction for multilayer brain networks using phase oscillator models', started to characterize structural layer properties such as their clustering or ranking. Nonetheless, the peculiarities and the interplay between such nodewise and layerwise characterization remain to be elucidated. To bridge this gap, our contribution is to propose the nodelayer duality, an unifying framework to study the structure of nodewise and layerwise representation in multilayer networks. We show that properties of both representations are not trivially related when considering the second moment of their distribution ,i.e. , the variance. Based on the local connectivity, we adopt a perturbative approach to track how each link toggle is reflected in one or both of the dimensions. This approach allows us to establish a formal duality between the two dimensions based on a distance between a multilayer network and its perturbed version. In return, it quantifies the situations in which one dimension better discriminates two multilayer networks compared to the other one. Specifically, we derive a complete closed-form of the formal node-layer duality in the case of random multilayer and multiplex networks.

Next, we provide a first characterization of di erent real-world multiplexes using the node-layer duality based on the aforementioned perturbation approach. We find that the real-world multiplex connectivity structure in both dimensions seems to be characterized by the spatial nature of each network. Also, results suggest that classes of multiplexes could be characterized by their peculiar interplay between the nodewise and the layerwise dimension. Inspired by [ ], we also propose a visualization approach to represent the connectivity of multilayer networks in the nodewise and layerwise dimensions.

In this e ort to characterize real-world systems, we then focus our work on multilayer brain networks. Among real-world systems, multilayer brain networks rapidly emerged as an interesting new candidate for multiscale modeling of the brain (Fig. . ). In particular, the formalism allows to study the complex interplay between structural interactions and functional interactions within the brain. It also allows to understand how the di erent modes of brain activity are reorgnized during a cognitive task or a disease and how they are related to behavior in an unprecedented way. For instance, multifrequency brain networks recently have shown promising results to characterize Alzheimer's disease [ , ]. Also, the formalism seems well suited to study the cross-frequency coupling, i.e. , interactions between brain frequencies. Although integrating cross-frequency coupling was already suggested to better characterize diseases [ ], multilayer brain networks integrating cross-frequency coupling are still relatively unexplored [ , ].

In this perspective, we consider that the node-layer duality is a natural framework to study such systems by analyzing simultaneously regions (nodes) and frequencies (layers) to characterize Alzheimer's disease based on local connectivity. Relying on local and global measures on the connectivity, we find that the layerwise dimension provides better descriptors than the nodewise one in terms of the severity of the cognitive decline. In particular, we recover the connectivity disruption in the alpha band ( -Hz) which is a well-kown feature of Alzheimer's disease. All in all, it validates the need to systematically characterize the layerwise dimension in multi-frequency brain networks. The manuscript is structured as follows:

In Chapter , we briefly introduce the basic concepts and tools in network science that are of relevance for this manuscript. In particular we briefly review network comparison methods which provide an important background for group-wise comparison between normal and abnormal brain networks. Then, we describe how concepts and tools of network science have been generalized to multilayer networks. We also briefly review the literature that started to characterize layers in multilayer networks. In Chapter , we first present an overview of the neuroimaging techniques that constitute the basics for the contruction of brain network. Then we describe the di erent ways to define nodes and links in brain networks. In particular we highlight the peculiarities of defining either structural or functional links. Finally we provide a brief overview on the main results related to brain networks. Specifically, we review the network characterization of Alzheimer's disease which lies at the core of our data-driven work.

In Chapter , we extensively review the literature on multilayer brain networks. We first present a brief overview of the multiscale modeling of the brain before introducing multilayer networks. Then, we present the peculiarities of multilayer brain networks as compared to classical brain networks in terms of structure and filtering techniques. Next, we detail the associated literature under the following pilar questions: the relation between structure and functions, the balance between integration and segregation and the behavioral adaption to tasks or learning. Finally we detail the contribution of multilayer networks to the characterization of brain diseases. In particular, we put the emphasis on the multilayer characterization of Alzheimer's disease. In Chapter , we present our theoretical contribution, i.e. , the node-layer duality framework. We present its principles and where it is relevant for a dual structural analysis of multilayer networks. Then, we show that the dual analysis of multilayer network retains complementary information about the position and the number of links in the system which is the basis of our perturbation approach. To apply it, we present a stochastic rewiring model that is tailored to switch links in a way that emphasizes the connectivity di erences between the two dimensions. We present the analytical derivation of the expected multidegree centrality sequence after a rewiring by the model. We validate the model by comparing it with an equivalent numerical rewiring. In the limit of large networks, we derive a formal node-layer duality based on the notion of distance between a network and a particular perturbed version of it. We prove and numerically validate a closed-from of this formal duality relationship in the case of random multiplex and multilayer networks. In particular, we show how this relationship scales with the size of the network. We finally discuss the potential implications of the node-layer duality framework for multilayer network analysis. In Chapter , we present an application of the node-layer duality framework. Especially we test the hypothesis that the node and layer multidegree centrality distributions uniquely characterize classes of real-world multiplex networks ranging from transportation networks to social ones passing through neural networks. For each network, we compute the distance between it and a perturbed version of it, normalized by its equivalent random null counterpart to avoid size e ects. Using k-means clustering we show that the multiplex systems are well separated into two clusters that seem to distinguish between networks with a spatial connotation and others without. Aside, we provide a method to visualize the organization of connectivity in both dimensionswithin multiplex networks based on singular value decomposition. Finally, we discuss the potential of the nodelayer duality framework for characterizing classes of multiplex networks. In Chapter , we apply the node-layer duality framework to the study of cross-frequency coupling in Alzheimer's disease. We show that the nodewise and layerwise distances between the Alzheimer's disease patients' networks and the averaged healthy control have the same discriminative power between the two conditions. However as the frequency resolution becomes coarser the layerwise dimension discriminates patients better. We show that indices of global connectivity disruption correlate significantly with the severity of the cognitive decline in the layerwise dimension but not in the nodewise one. Similarly the disruption of local connectivity of specific regions and frequencies correlates with the cognitive decline although a higher number of strongest correlations were found in the frequency dimension. Notably, the strongest correlations include the caudal anterior cingulate gyri and the frequencies of the alpha band. Finally, we discuss the observed comparative advantage of the layerwise dimension since Alzheimer's disease seems to be characterized not only by the connectivity disruption in specific brain areas, but more importantly by the aberrant coordination between frequencies.

In Chapter , we draw a conclusion and a perspective on the significance of our work in the broad picture of multilayer networks in general and multilayer brain networks in particular. We also relate our work to the general concept of duality.

: In mathematics, V and E are often called the vertex set and the layer set respectively : When = 0 the graph is said to have no self-loops.

Network Science

. Classical networks

In this chapter, we introduce the basic concepts and tools of network science. Then, we show how these concepts and tools generalize in the context of multilayer networks. Specifically, we put the focus on the principle of network comparison and layer characterization which we will be used in chapter to .

. . Basics and concepts

A network is the wiring map of any system whose components interact through pairwise interactions. The mathematical representation of a network is a graph G = ( , ), where = { | ∈ [1, ]} is the nodes set and E is the links set with ⊆ { |( , ) ∈ × }. . The element of E represents the interaction between node and node . The number of components of the network, i.e. , the number of nodes is denoted N | | and its number of interactions, i.e. , the number of links, L | |.

The definition of a graph can be refined in two ways. First, if the graph is combined with a function : → ℝ that associates to each pair of nodes , a value , the graph is defined to be weighted, otherwise it is defined to be unweighted. Weighted graphs take into account the di erence of intensity in the interactions between nodes. Second, if it exists elements of E such that ≠ , meaning that the interaction from nodes to is di erent from the interaction from nodes to the graph is defined to be directed, otherwise it is defined to be undirected. Directed graphs take into account the directionality of the interactions between nodes. In the case of weighted directed graph, it translates to ≠ .

Any graph G = ( , ), can be represented by its adjacency matrix. Let G being a weigthed directed graph with N nodes. Its adjacency matrix is a × matrix whose elements read:

= if there is a link from node i to j 0 if there is no link between i and j ( . )

In the case of undirected graph, is symmetric since = . For undirected graph, = 1 if there is a link between i and j.

The adjacency matrix o ers an operational and completely determined representation of a graph. It allows to use the mathematical apparatus of linear algebra on graphs like the Perron-Frobenius theorem or the spectral theorem. In an undirected graph, the number of links L read:

2 = , ( . )
where the factor comes from the fact that each link is counted twice since = . In the case of directed network = , because the presence of a link from to does not necessarily imply a link from to . A fundamental aspect of network modeling is to correctly define the nodes and links which compose the system of interest.

. . Measures on networks

A major part of network science relies on deriving network measures in order to characterize network properties. In the following, we expose measures that are classical in the field. For the sake of simplicity, we use systematically the term networks even to designate graph.

Density.

The density of a network is the ratio of the total number of links over the maximum number of possible links in the network . For undirected graph the density read:

= = 2 ( -1) , ( . )
where = 2 , i.e. , the number of possible pair with N elements, in the undirected case. In the directed case, the maximum number of possible links is egal to ( -1), the number of o -diagonal elements in the adjacency matrix. Therefore = ( -1) .

For = 1 the network is said to be complete. For → 0, the network is said to be sparse.

Node degree. The node degree is the sum of the number of links connected to a given node. For undirected and unweighted network the node degree read: = ≠ ( . )

In the undirected weighted case, the previous formula is called the node strength, the sum of weights connected to a given node. For undirected weighted graph, the degree read = ≠ Θ( ), where Θ is the Heaviside function.

In directed unweighted network, we distinguish between the sum of links that terminate to a node, i.e., the in-degree and the sum of links that originate from a node, i.e. , the out-degree:

= ≠ = ≠ ( . )
In the rest of the manuscript, we will focus on the properties of undirected unweighted networks. ), 'Persistent patterns of interconnection in time-varying cortical networks estimated from high-resolution EEG recordings in humans during a simple motor act', Node degree distribution. The node degree distribution is the probability distribution to find a node of degree in the network (Fig. . ):

( = ) = , ( . )
where is the number of nodes of degree equals . Note that ( ) = 1. A node that has a high degree is called a hub. A classical result in network science is that a vast majority of real-world networks are scale-free, i.e. , they follow a power-law distribution ( ( ) ∼ -) [ ]. Scale-free networks were reported in technological, social and biological networks just to name a few [ , ]. Such universality of the power law distribution prompted a considerable attention on the generative mechanisms [ ] and the advantages of scale-free networks (e.g., the robustness against random attacks [ ]). However, the claim of universality is debated due to the lack of statistical testing on early studies [ ] and the rarity of scale-free networks when such tests are applied [ ]. On the other side, statistical tests could not detect a power-law distribution due to finite-size scaling of the latter [ ]. All in all, the minimal baseline in this controversy could be that the "knowledge of whether or not a distribution is heavy-tailed is far more important than whether it can be fit using a power law" [ ].

Network motif. Network motifs are recurrent connection patterns involving few nodes and are therefore easily interpretable (Fig. . ). They constitute the basic building blocks of a complex system architecture which can code for essential biological functions such as autoregulation, cascades and feed-forward loops [ -] in biological networks. To ensure suitable comparisons, frequencies of appearance of network motifs are often compared with the ones of matched null models (see Sec. . . ).

PageRank centrality.

The node degree is a centrality measure of the first order as it takes into account only the immediate vicinity of a node. The PageRank centrality is a higher-order centrality measure which avoids to nodes related to a highly central nodes to be quite central whereas they are not by other regards. If is the PageRank where , , are a centrality term added to every node, the balance between the added centrality and the one coming from the structure of the network and the out-degree of node , respectively. For undirected networks becomes . Interestingly, PageRank centrality of node can be interpreted (and computed) as the probability of a random walker to be present at node at equilibrium. In this directed version, PageRank is at the core of the Google web search engine [ ].

Clustering coe cient. In social networks (among others), we often observe that individuals form dense clusters of relationships. The local clustering coe cient measures the tendency for two neighbors of a given node to be also connected to each other. More accurately, it measures the ratio of closed triangles containing a node and two of its neighbors over the total possible pairs of neighbors.

For a node it read:

= 2 ( -1) , ( . )
where , are the number of closed triangles from and the node degree of node , respectively.

To obtain the global tendency of a network to form triangles, we can average the previous metrics to obtain the average clustering coe cient:

= 1 ( . )
Shortest path. An important question in network science is to know how far each node is separated from each other. In other words, it is important to know how convenient it is to navigate across the network by following path between nodes. The shortest path is defined as the smallest number of links that separates node and node . The characteristic path length is the average shortest path for all possible pairs of nodes in the network:

= 1 ( -1) ; ≠ ( . )
In practice, to find the distance between nodes, we use the Breadth-Search algorithm whose time complexity is linear in the number of Betweenness centrality. Based on shortest paths, the betweenness centrality of a node is its tendency to connect structurally connected distant parts of the network [ ]. More precisely, the betweenness centrality of node i is the proportion of shortest paths between any nodes j and h that pass through i: ) where ℎ ( ) is the number of shortest paths between h and j that pass through i and ℎ is the total number of shortest paths between h and j.

= 1 ( -1)( -2) ℎ≠ ,ℎ≠ , ≠ ℎ ( ) ℎ , ( .
Interestingly, the ranking given by the node betweenness centrality can be di erent from the one given by the node degree. It means that the higher degree node is not necessarily the most central in terms of betweenness [ ]. Under the hypothesis that information circulates on shortest paths on networks, the betweenness centrality informs on nodes that could control the passage of information or that can be a bottleneck for it [ ].

Local and global e ciency.

The local e ciency can be seen as a generalization of the clustering coe cient:

( ) = 1 ( -1) ,ℎ∈ 1 , ( . )
where is the subgraph including all nodes that are in the immediate topological vicinity of node . The measure is normalized to take into account all the possible pairs around node . Unlike the clustering coe cient, it takes into account direct and indirect paths that can link nodes in the vicinity of node . In this sense, it measures how well the neighborhood of node is integrated and resilient under the disruption of the aforementioned node. The global e ciency measures how well the information exchange works when passing by shortest paths in the network at the global scale [ ]. It read:

= 1 ( -1) ≠ 1 ( . )
Communities and modularity. In real-world systems , it is common to find communities, i.e., a group of nodes that are more densely connected to each other than to the rest of the network (Fig. . ). Detecting these communities, i.e. , detecting the number and the composition of the relevant partition in a network is an important task in network science since their properties can supposedly deviate from the average properties of nodes. In particular, communities can possess information processing or functional properties in a network [ ].

A widespread method to quantify the quality of the community detection is the modularity index. Its principle relies upon the hypothesis that networks where links are attributed at random do not have a ), 'Detecting richclub ordering in complex networks', community structure. Therefore, it quantifies how far the community structure of a network is from the one that would be expected by chance if the network would be random. More precisely, in its most popular declination, the modularity compares the network of interest to a network that has the same degree distribution but is completely random otherwise (see Sec. . . ). In this case the modularity read:

= 1 2 ( - 2 ) ( , ), ( . )
where ( , ) is the Kronecker delta function that is equal to one if and belongs to the same community, i.e., = and zero otherwise. is inferior or equals one. When it is close to , it indicates that the community partition is maximally di erent as compared to the one of a random network.

The modularity maximization is a community detection method that is based on the idea that the more maximized is the modularity the better is the associated partition of the network. Since finding the global optimal partition is a NP-hard problem , modularity maximization relies on greedy algorithms whose most popular is the Louvain algorithm [ ] associated with a time complexity that scales linearly in the number of links in the network ( ).

It is worth noting that modularity maximization su ers from some limitations like, for instance, a resolution limit under which a community is Rich-club phenomenon. The rich-club phenomenon is tendency for nodes with the highest degree to be more interconnected than expected by chance [ ]. To quantify this phenomenon, the rich-club coe cient at degree is the density of links of nodes with degree at least , > over the number of nodes with degree at least , > :

( ) = 2 > > ( > -1) (
( ) = ( ) ( ) , ( . )
where ( ) is the average rich club coe cient of an ensemble of random networks with the same degree distribution as the the one of ( ) (see Sec. . . ). ( ) > 1 indicates that the rich-club phenomenon is higher that what could be expected by chance.

Laplacian matrix. By modeling linear di usion processes on networks, the Laplacian matrix, naturally arises as the graph equivalent of the Laplace operator in calculus. It read:

L = -, ( . )
where is the degree matrix, i.e. , a diagonal matrix where theℎ diagonal entry is the node degree of node . L has many implications from denoising to low-dimensional embedding [ ]. In particular, the eigenvalues of L are informative on the way to partition the network into two parts or into communities [ ]. By construction L, the smallest eigenvalue is equal to zero. The second smallest eigenvalue is called the algebraic connectivity and relates to many properties of the network. For instance, its inverse is proportional to the time required to reach equilibrium in linear di usion processes and to reach synchronization for linearly and di usely coupled oscillators [ ].

.

. Null models in network science

Null models. Many measures on networks scale with the number of nodes and the number of links . For instance, higher clustering coe cients arise on networks with a high number of links. Actually, it increases the probability of making triangles by chance, without being informative on the clustering properties of the system.

The purpose of null models is to provide a reference model to the system under interest to disentangle what is an essential property from what is a mere consequence of contingent properties in the system [ ]. Concretely, null models provide rules to build an ensemble of synthetic networks that satisfies some desired properties. Ideally, the null model should match the properties of the system except the one under investigation which is not possible in reality. Therefore, only basic properties are often match like the number of nodes, the (expected) density or the degree distribution. Then, the comparison between real-world systems and the synthetic networks allows to assess whether the real-world properties can be explained by the null models or not. Since null models are also generative models of networks, they can provide explicability on the arising of patterns or properties in real-world networks. and the probability ∈ [0, 1] for any nodes to connect to each other. By construction, the degree distribution of the , is binomial:

Erdös

( = ) = -1 (1 -) -1- ( . )
For large and sparse networks which are closer to real-world systems, the degree distribution can be approximated by a Poisson distribution:

( = ) ≈ ( ) - ! ( . )
In any case, the node degrees are located at a single-scale around the mean with few fluctuations. In terms of generative model for real-world network, it limits its applicability since a vast majority of empirical degree distributions are right-skewed (see Sec. . . ).

However, its simplicity allows to derive analytically many of its properties such as its clustering coe cient or its characteristic path length [ ]. Furthermore, it constitutes a good minimal null model to investigate the dependence in the number of nodes and density of a property of interest. We can use the expected properties of an ensemble of Erdös-Renyi network to normalize a measure (e.g., the network motifs) on real-world networks. If the normalized value is close to one, it means that the property under interest in the real-world network is not di erent from what would be expected by chance. If the normalization is greatly di erent than one, the investigated property is deemed to be a characteristic of the system. Configuration model. The configuration model is a null model that allows to completely randomize a network while preserving its degree sequence, i.e., the exact degree of each node in the network. This enables to study the implications of particular degree distributions (e.g. , scale-free) on the properties of interest. The Maslov-Sneppen algorithm is a rewiring procedure that allows to generate the configuration ensemble from an initial network. It consists in successively selecting two pairs at random and then in rewiring these two between the four implicated nodes while keeping their node degree sequence fixed (Fig.

. . Network comparison

A major challenge in network science is to quantify the similarity or dissimilarity between two networks (or group of networks). Applications are numerous in biology such as comparing between networks of diseased proteins with healthy ones which can give insights on how to cure the disease or comparing genetics networks between species to give clues on the di erent evolutions' dynamics [ ]. In neuroscience, the statistical comparison of a population of diseased brain networks and a population of healthy ones can deepen our understanding on a given disease [ , ]. On the theoretical level, comparisons between real-world networks and synthetic networks generated with null models are crucial to select the potential generative mechanisms upon which empirical networks are built.

Theoretically and historically, network comparison relies upon the graph isomorphism problem [ ] which consists in finding whether two networks are isomorphic or not . To study real-world networks, the graph isomorphism problem su ers from two theoretical drawbacks. First, it is likely that there exists no method to solve it in a polynomial time [ ]. Second, real-world networks can be reconstructed from noisy, incomplete data which makes it likely that two instances of the same network would be deemed to be not isomorphic with no mean to quantify this noise-induced di erence.

Therefore, from the graph isomorphism problem, network comparison evolved towards defining heuristics which give a measure of distance that quantify the dissimilarity between graphs by focusing on its relevant aspects. In this framework, a distance equals to zero between two networks indicates that they are the same, regarding the property measured by the heuristics. Given two networks 1 and 2 , a measure 'M' having two arguments, denoted by ( 1 , 2 ), performing a mapping

from 1 × 2 → ℝ, is called a 'bivariate network measure'. ( 1 , 2 )
is either the similarity or the dissimilarity (distance) between two networks, hence, we use the symbol ( 1 , 2 ) to indicate the similarity and ( 1 , 2 ) for the distance. The definition of heuristics is a trade-o between their time-complexity and the interpretation and the relevance of their results regarding the domain where they apply [ ].

The comparison methods largely fall into two categories depending on the scale of network features: global and local. Global methods assess a global property of a network (e.g. , the average clustering coe cient, the average shortest path or the degree distribution). On the contrary, local methods assess property at the scale of subgraph or graphlet [ , ] or at the node scale (e.g., with the degree). Then, comparison methods can be classified according to di erent classes of distances based on global or local methods.

In the following, we focus on classes of comparison methods. In particular we focus on heuristics that requires networks to have the same number of nodes. Should the reader be interested in delving into network comparison methods, we refer him/her to the following review [ ]. ), 'Temporal networks', Euclidean distance. Let 1 and 2 be two networks to be compared associated with their adjacency matrices 1 and 2 respectively. One of the most naive comparison methods is to take the Frobenius norm of the di erence between 1 and 2 :

( 1 , 2 ) = || 1 -2 || = (( 1 -2 )( 1 -2 ) ) ( . )
Thought simple, this distance was shown to be highly discriminative as compared to much sophisticated one in a social and brain networks context [ ].

Structural descriptor-based distance. This class of distance is based on structural properties (e.g. , degree, clustering coe cients. . . ) that are mapped one-to-one between networks. One of the simplest measures of this kind is the degree-based distance. Let 1 and 2 be the degree sequence of each node of the networks 1 and 2 respectively, the distance read:

( 1 , 2 ) = || 1 -2 || 2 ( . )
Though easily computable, this distance reduces to one simple aspect of the structure of each network. Therefore, very similar networks according to this measure can be actually very di erent. For example, the previous distance between a network and its reshu ed counterpart with the configuration model will be equal to zero albeit they can be otherwise very di erent. Therefore, structural descriptor-based distances are usually composite distances that take into account several, hopefully irreducible aspect of the networks [ ].

Graph kernels. Graph kernels are functions : × → ℝ, that are positive semi-definite and therefore can be defined as an inner product ( ) where S is a vector of network features. Like the traditional inner product, they are measures of similarity, i.e., networks are maximally di erent where ( , ) = 0. In most of graph kernels, is defined as the frequency of subgraph patterns in the networks (e.g., it counts the di erent subgraphs of node size ). The most important di erence between graph kernels and previous methods is that ( ) can be defined explicitly with having an arbitrary dimensionality. Therefore, networks can be compared by algorithms-made features instead of human-made features [ ].

. Multilayer networks Quite rapidly in network science emerged the need to model systems by structures of the higher order compared to networks [ , ]. Actually, in many systems that can be represented by networks, there exists multiple types of interactions between their constituents. Within a social network, people have the opportunity to share information by means of direct personal interactions, telephone conversations, written correspondence, or various online social platforms. Each of these layers of interaction possesses a distinct framework, which is influenced by its inherent nature and the specific time frame it operates within. The conceptual advance of multilayer networks is to consider these layers in interaction which makes emerge new properties that are crucial for the understanding of systems.

In the next section, we define the di erent types of multilayer networks and how network measures and null models generalize to them. Finally, we insist on the characterization of layers, the emerging dimension of multilayer networks compared to single-layer networks.

. . Types of multilayer networks The need to investigate complex systems with multiple types of connectivity has emerged, almost independently, from di erent disciplines including social science, engineering, and computer science [ -]. More recently, the physics community also produced pioneering works on various notions such as networks of networks [ , ], node-colored networks [ , ] interdependent networks [ , ] or multilayer networks [ , ]. As a consequence, di erent terms have been introduced and adopted, thus producing a lack of a consensus set of terminology and mathematical formulation. Only in the last decade, we have eventually witnessed the dawning of general frameworks compatible with tools from complex systems and network science [ , ] is the set of nodes of and the set of its links, with ⊆ × . The set of links between the nodes of the graphs at di erent layers and is denoted by

C = { ⊆ × | ≠ }.
Similar to classical single-layer networks, a convenient and operational representation of a multilayer network is given by the so-called supraadjacency matrix. Let Gbe a binary undirected multilayer network with nodes and layers. Its supra-adjacency matrix is a × block matrix whose elements read:

=
1 if there is a link from node i in layer to node j in layer 0 otherwise ( . ) Hence, results in a matrix with blocks on the main diagonal, accounting for the connections within layers, and ( -1) o -diagonal blocks describing the links between di erent layers.

The above definitions are quite general and allow to describe complex systems exhibiting di erent number of nodes in each layer or scale, directed or undirected interactions, as well as weighted or unweighted connectivity. We consider here multilayer networks composed of replica nodes, i.e. , where all the layers have the same number of nodes representing the same units of the system across di erent layers. Note that this is a very particular condition, which however matches the nature of data presented in the large part of the studies so far. In this configuration = , ∈ {1, ..., } and only the connectivity within and between layers is allowed to change. In the following, we will refer to these general configurations as to multilayer networks (Fig. . . . where is the × identity matrix.

Multiplex networks are the most studied structure so far and pervade the field of complex systems [ , ]. As an example, they can model the di erent types of interactions in a social network whether it is characterized by means of communications (phone, email, face-toface...) or di erent social ties (friendship, enmity...) [ , ].

Based on the above configurations, many types of multiscale interconnected systems (e.g., spatial, temporal, multimodal) can be represented and investigated. For example, temporal networks are represented by a particular type of multiplex, where only replica nodes between adjacent layers are interconnected, and the blocks after the first diagonals in Eq. . become null matrices (Fig. . c). We notice that in general, the information contained in multilayer networks can be obtained neither from equivalent aggregated versions (e.g., where links are averaged across layers) nor from standard network metrics and tools [ , ]. For this reason, it is crucial to derive new concepts and methods to quantify the higher-order topological properties emerging from multilayer networks. In the next subsection, we introduce some of the metrics and tools that have been developed so far. For the sake of simplicity, we will focus on unweighted and undirected multilayer networks.

.

. Multilayer network measures

In the following, we briefly present some of the multilayer methods that are classical in the field, whith a particular focus on those that have been most frequently used in neuroscience (see Ch. ).

Overlapping degree. The overlapping degree is the sum of links connected to a node within layers and across all layers, it read:

= ≠ ( . )
For weighted multilayer networks, this measure translates into the overlapping strength.

Multilinks and multi-degree. Layers in real-world networks tend to be significantly correlated between each other, i.e., a given pair of nodes are connected through multiples layers. For instance, in a social multiplex network where layers are modalities and nodes are individuals, people that connect in the real-life interactions will likely to connect by phone or social media as well. To take into account these correlations, multilinks were defined as the -uplet ( being the number of layers) where a given pair of node and are connected. For instance, for M , if and are connected by multilink( , ), it means that they : The multi-degree, with a dash, should not be confused with the node and layer multidegree centrality that are introduced later in the manuscript (see Ch. The multi-degree

- →
is the number of multilinks -→ incident to node . Multiplex participation coe cient. Nodes of similar overlapping degree can have dramatically di erent roles across layers. Actually, one node can be a hub in one layer and almost isolated in the other ones whereas the links of other nodes can be equally distributed across layers. Furthermore, the degree distribution of each layer is heterogeneous. To take into account this heterogeneity and the contribution of nodes to each layer, the multiplex participation coe cient is defined as [ ]:

= -1 [1 - ( ) 2 ], ( . )
where is the degree of node at layer . When = 0, the node's links are all concentrated in one layer; when = 1, they are uniformly distributed across layers. ) which takes into account the possibility to form triangles by means of links belonging to two di erent layers. ), 'Multiplex core-periphery organization of the human connectome',

Multiplex clustering coe cient. A relatively straigtforward extension of the clustering coe cient is the multiplex clustering coe cient

[ ]: = ≠ ≠ , ≠ ( -1) ≠ , ≠ , ( .

Overlapping betweenness centrality and multiplex

= 1 ( -1)( -2) , ≠ , ≠ ( ) , ( . )
where ( ) is the number of shortest paths from node to passing through node in layer and is the total number of shortest paths between node and in layer .

PageRank centrality is another centrality measure that was adapted to multiplex networks [ ]. In multiplex networks, random walkers have the possibility to jump to adjacent nodes and teleport to nodes in other layers, according to a modified version of the transition probability [ ].

Multiplex networks motifs. When dealing with multiplex networks, motifs can be formed by links belonging to di erent layers [ ]. Hence, the total number of possible configurations depends on the number of layers but also on the type of interaction, (e.g., negative or positive). In these cases, -scores are typically used to determine the statistical abundance of a multiplex motif according to the following formula:

( ) = ( ) -¯ ( ) ( ) , ( . )
where is the occurrence frequency of a given multiplex motif, while ¯ ( ) and ( ) are respectively the mean frequency and its standard deviation obtained from a set of equivalent random multiplex networks R.

Multiplex modularity.

In the case of multiplex networks, the definition of modularity incorporates the relation between di erent layers and partitions all the layers simultaneously [ ]:

= 1 2 [( - 2 ) + ] , , ( . )
where is the total number of links in the multilayer, sets the granularity of the community structure in each layer, is the total number of links in layer , is a parameter that tunes the consistency of communities across layers and , = 1 when node in layer and node in layer belong to the same community, and zero otherwise. Maximization of is finally obtained by similar heuristic as in the singe-layer one. In temporal networks, nodal metrics reflecting mesoscale network properties can be defined by measuring, for example, the node flexibility, i.e., the average number of times that a node changes community assignment across layers [ ].

Multiplex core-periphery. A peculiar network partition consists in separating the network in a core of tightly connected nodes, and a periphery made by the remaining weakly connected nodes [ ].

Similarly to a rich-club, the presence of a core is crucial for the e cient integration of information between remote parts of the network [ , ].

[ ] introduced a fast core-periphery detection algorithm for multiplex ), 'Small-world behavior in time-varying graphs', networks .Based on local information [ ], the method first defines a multiplex richness of a node by combining its degrees in each layer. Nodes are then ranked according to their multiplex richness values and the core-periphery separation is given by the following optimal rank [ ]: * = ( + ) , ( .

)
where + is the richness obtained when considering only the links of the node ranked in position r towards nodes with higher ranks.

In the case of weighted multiplexes, the coreness of a node is given by the number of times it belongs to the core after filtering and binarizing the network's links with a range of di erent threshold values.

Overlapping global and local e ciency. Based on topological distances, one can also quantify the global tendency of a multiplex network to form highly clustered and e cient groups via the overlapping local-e ciency [ ]:

= 1 ( -1) , ≠ ∈ 1 ( -1) 1 ( , ) , ( . )
where is a sub-graph containing the neighbors of node i and ( , ) is the length of the shortest path between node i and j at layer . The multiplex extension of global e ciency consists in computing the shortest paths also across layers.

Characteristic temporal path length.[ ] extended the concept of topological distance to temporal networks by introducing the characteristic temporal path length , which measures the formation of shortest paths across consecutive layers.

Same authors also introduced a metric to quantify the probability that the neighbor set of a node that is present at time is also present at time + 1. By averaging over all the nodes, they eventually defined the temporal-correlation coe cient C as:

= 1 ( -1) =1 -1 =1 +1 ( )( +1 ) ( . )
Together, the last two global metrics measure how the system information is respectively integrated and segregated over time and can be used to assess the small-world properties of time-varying networks.

Supra-Laplacian matrix..In a full multilayer network, 2 is calculated from the associated supra-Laplacian matrix, whose elements are defined as: ,where is the is the total number of links incoming to node i at layer .

L = - , if i j, = - , otherwise ( . )
In multilayer networks, 2 is sensitive to the amount of intra-and inter-layer connectivity, and typically quantifies the integration/segregation balance among layers from a dynamical perspective [ , ]. Notably, 2 exhibits a phase transition when increasing the interlayer connection intensity, from layers being independent/segregated to a high overall dependence/integration [ ].

.

. Null models of multilayer networks

As for single-layers networks, it is necessary to compare the properties of real-world multilayer networks (e.g., multiplex motifs or dynamical properties) to random multilayer network ensembles. An important aspect of multiplex networks is that layers can be correlated. Therefore, we can distinguish between null models that preserve the layer structural correlations from the ones that do not. The configuration model (see Sec. . . ) has been adapted to multiplex networks along these two lines.

Configuration model with independent layers.. Since the layers are not required to be correlated, one of the simplest ways to generate multiplex networks with fixed layerwise degree sequence (i.e., fixed overlapping degree sequence) is to apply a configuration model on each layer of the multiplex network of interest. This model was shown to generate a multiplex network ensemble where layers are independent, i.e. , they do not show a significant overlap [ ]. Like for single-layer model, null models of multilayer networks can fix several other properties of the networks of interest (e.g., the community structure). The multilayer networks of the null ensemble are either built to strictly respect those constraints or only on average.

Configuration model with multi

. . Layer characterization

Historically, multilayer network analysis focused on the emergent properties of nodes when considering the di erent modes of connectivity of a system in interaction. Although still poorly explored, part of the literature started to characterize the layers, i.e., the emergent object in multilayer network as compared with single layer ones. In particular, studies developed along two main lines of forces. How similar layers are, i.e. , how do they cluster ? how central do layers are, i.e. , how to rank them ? [ ] developed a stochastic block model that cluster layers into groups called 'strata'. More specifically, their model co-cluster nodes to community and layers to strata. They show that if the multiplex networks are built from strata of layers, their model is more e cient than stochastic block models applied on the all network or on layer to layer to recover strata. Another approach consists in clustering layers through their redundancy as assessed with an entropy-based measure [ ]. More specifically, [ ] used Jensen-Shannon (JS) divergence [ ] between couples of layers and used their measure to hierarchically cluster layers in a human microbiome multiplex network, with similar results as in the literature. Similarly, [ ] built a measure of layer similarity based on the entropy of network ensembles constrained by some properties (Fig. . ). Through this approach, they discern how the organization of knowledge in physics manifests itself within the structure of collaboration networks. Finally, [ ] considered a measure of link overlap between layers. From that , they build a single-layer network whose interacting agents are the layers linked by their similarity measure. Then, they form cluster using regular tools of node community detection in single layer networks. For example, their approach grouped airlines that are based in the same regions into layer communities in an airline network. ), 'Tunable Eigenvector-Based Centralities for Multiplex and Temporal Networks', introduced a version of the eigenvector centrality that ranks nodes and layers of undirected multiplex networks in a coupled fashion based on multi-homogeneous maps . More precisely, they show that a multi-homogeneous map is the unique existing solution of a nonlinear equation where the centrality of a node depends on a product of the centrality of its neighbors times the centrality of the layer in which the nodes are connected. Similarly, the centrality of the layers depends on the product of the centrality of each pair interacting in this layer. These centrality measures provide other ranking as compared with classical eigenvector centrality measures on real-world multiplex. Elaborating on the principle of centrality by association, [ ] presented the MultiRank algorithm to rank both nodes and layers by the mean of coupled equations. The centrality is estimated through PageRank algorithm where the probability to hop to a neighbor node of a given layer depends on the centrality of the layer. Similarly, the centrality of a layer is determined by the centrality of nodes that are connected in it. The parameters of the model make it versatile to evaluate di erent types of centralities. For instance, it is able to consider or not consider the total weight of layers to determine their influence. Similarly, it can favor (or not) layers with high number of highly central nodes. Finally, it applies to weighted, directed multiplex networks.

Regarding the development of centrality measures for layers, [ ]

In the same fashion, [ ] recently proposed Co-Rank algorithm, based on a random walk, to estimate simultaneously the centrality of nodes, layers and times in a temporal multilayer network.

Finally, [ ] developed the supracentrality matrix, a block matrix whose diagonal blocks are the centrality matrices related to each layer (e.g. the adjacency matrices of each layer in the case of the computation of eigenvector centrality). The o -diagonals blocks represents the diagonal coupling that are tuned by a parameter. Once the vector of centrality is computed over the node-layers, it is then summed over correct dimensions to obtain either the centrality of layers or the centrality of nodes. Authors focus on the e ect of interlayer coupling on the node and layer centrality. They found that a weak coupling yields a layer decoupling and a strong coupling yields a layer aggregation e ect in terms of centrality. It emphasizes on the value of the interlayer coupling to determine nodes and layer centrality. In particular, for the European Airline multiplex [ ], they found that nodes that are important for the weak and strong coupling receive a boost of centrality at intermediate regimes compared to the ones that are highly central in only one of these limits. with a strong magnetic field. Then, radiofrequency waves are applied in short pulse of microseconds and their energy are absorbed by nuclei which deviate their spin from the former magnetic field. Next, nuclei released this energy up to thermal equilibrium in the so-called relaxation. The time it takes for nuclei to release their energy to thermal equilibrium, i.e., the relaxation and the amount of energy released, is informative on the chemical environment of the nuclei. : White matter are bundles of myelin that connect grey matter areas.

Network Neuroscience

: Grey matter regroups the brain areas that are mostly composed of neurons.

[ ]: Lazar ( ), 'Mapping brain anatomical connectivity using white matter tractography', brain anatomy while T is better suited to identify lesions (Fig. . ) Their high spatial resolution allows to distinguish the cerebrospinal fluid (long relaxation time, appears darker on the image) from white or grey matter (short relaxation time appears brighter on the image). The typical voxel size for a T -weighted MRI is about 1 3 .

Di usion weighted imaging: Di usion weighted imaging is based on the three-dimensional di usion (Brownian motion) of water molecules in the brain. Without constraints, water molecules di use in every direction, i.e., in isotropic directions. However, in certain brain tissues, such as white matter, the water di usivity is constrained into a preferred directions, i.e., in anisotropic directions that influences the nuclei properties observable by MRI. In particular, water di uses preferentially in the direction of the white matter fiber. Therefore, by measuring the di usion profile voxel-by-voxel (a process called Di usion Tensor Imaging, DTI), it allows to obtain the pathways of white matter fibers in the brain [ ] (Fig. . ).

Functional MRI: Neurons need oxygen to function. Therefore, at the capillari level, there is an increase of the cerebral blood flow and a change in the oxygen concentration when neurons fire. Yet, the NMR properties of hemoglobin molecules carrying oxygen are di erent from the ones which do not carry it. This is called the blood oxygen level dependent (BOLD) contrast and it allows an indirect measure of the brain functioning by targeting regions that are more in demand of oxygen as compared to the In the literature, frequency bands are associated with di erent tasks or behavior. band ( . -Hz) and ( -Hz) are observed when subjects are in deep or light sleeps, respectively. band ( -Hz) is generally associated with resting states (with eyes closed). band ( -Hz) is observed during motor tasks for instance. Finally, band ( Hz) is detected during specific highly demanding tasks.

Di erences between EEG and MEG are diverse. First, EEG is more sensitive to radial currents, i.e., activity located at the gyrus level whereas MEG is more sensitive to tangential currents ,i.e., activity located at the sulci level (Fig. . ).

Second, the EEG signal is highly deformed by crossing the di erent layers of materials (e.g., the skull) between the sources and the sensors : MEG signal is 10 9 weaker than the Earth magnetic field. ), 'Small-world and scale-free organization of voxel-based resting-state functional connectivity in the human brain', whereas MEG signal is less sensitive to this feature. Third, the MEG signal has a slightly better spatial resolution (around cm) than the EEG (around -cm). The drawbacks of EEG are compensated by the fact that it is a cheap and portable device as it mainly consists in placing electrodes on the scalp (Fig. . ). On the contrary, MEG requires to place the subject under a cumbersome and costly machine made of magnetometers and gradiometers (which measure the gradient of magnetic field) in a magnetically shielded room .

Unlike MRI data, M/EEG signals are not associated with a source in the brain but they are associated to a mix of signals recorded at the level of the scalp sensors. In order to reconstruct the brain sources from the recorded data, one has to solve the direct and inverse problems [ ].

The direct problem is about defining a physical model of the source and a model of how the source-generated electromagnetic field is recorded at the scalp level. The aim is to model the electromagnetic field produced by a cerebral source with known characteristics. For that purpose, it is necessary to consider both a physical model of the sources, and a model that predicts the way that these sources will generate electromagnetic fields at the scalp level. Realistic models are geometrical reconstructions of the head materials (cerebrospinal fluid, skull. . . ) obtained directly from the MRI of the subject. Then each material is associated with a electric/magnetic conductivity.

The inverse problem consists in estimating the number of sources.

Since it has an infinite number of solutions, it requires further assumptions to be solved. For instance, the dipole modeling method sets an a priori number of sources and associate each of them to dipoles. The source reconstruction problem is made challenging by volume conduction e ects, i.e., the di usion of the electromagnetic field through the head tissues.

When no source reconstruction is performed, the M/EEG analysis is said to take place in the sensor space. Otherwise, the M/EEG analysis takes place in the source space.

. Brain networks reconstruction . . How to define nodes ?

Building a network requires crucially to determine what are the nodes, i.e., what interacts in the systems of interest. While at the micro-scale, neurons constitute the straightforward unit to be modelled as nodes, the macro-scale is left with a multiplicity of brain atlases to define the regions of interests (ROIs), i.e., the nodes at the whole-brain scale.

Following [ ], the ideal node should encompass similar cell populations, i.e., should be a coarse-graining of the neuronal architecture at the lower scale, i.e. , the mesoscale.

A natural starting point is to define nodes as the thinnest possible unit that we can extract from data. It would consist in taking voxels as nodes in MRI data or sensors for M/EEG data [ , ]. However, voxels or sensors data do not necessarily correspond to cell population. Even worse, it can separate a given cell population into two nodes which would lead to a spurious increase of connectivity between those nodes.

Another approach consists in defining the nodes directly as regions that share a similar cell population and organization under the microscope, i.e., a similar cytoarchitectonic architecture. Since the seminal work of Broadmann [ ], who mapped the human brain in regions according to their cytoarchitectonic architecture other related methods using the myelin architecture instead were tried but not on the entire brain. The major problem lies in the interindividual variability of the cytoarchitectonic architecture that can be up to a factor of in the volume of brain regions [ ].

An alternative relies on delineating the ROIs based on macroscopic landmarks recordable with MRI such as gyri and sulci. For instance, a popular parcellation scheme is the Automated Anatomical labeling (AAL) [ ] thought it cannot account for interindividual variability since it is based on a single brain. To take into account interindividual variability, the Desikan-Killiany atlas [ ] is generated from a probabilistic map (each point in the brain is assigned a probability to belong to a ROI) based on a population of brains that are diverse in terms of sex, age and condition. Again, these methods do not reflect the underlying cytoarchitectonic architecture areas that cut across sulci and gyri. Furthermore, nodes can represent regions with huge size di erences which can bias the connectivity estimation in favor of the larger regions. Indeed, these variations can be biologically grounded (e.g., when considering a hub area) but there is no guarantee that the defined nodes are meaningful.

Finally, one can start from connectivity to define nodes by postulating that each functional area possesses a unique connection pattern. Therefore, a measure of similarity (e.g., Pearson correlation. . . ) is defined for every pair of voxels that is then cluster to obtain functional ROIs. Though largely used for functional connectivity data, it is still limited in structural data by the capacity of MRI to track long-range pathways. One of the interests of these methods is that it can be computed individually on each scan of each patient, leading to consistent parcellation for subsequent connectivity analysis [ , ]. Since there exists no gold standard, the choice of nodes definition should be guided by the scientific question and cross-validation with di erent schemes are necessary to provide rigorous claims on the human brain.

. . Structural connectivity

Once nodes have been defined, the next step to build brain networks is to define and compute the connectivity estimator between nodes. From DWI-data, the potential trajectories of axonal fibers can be reconstructed using tractography algorithms . Such algorithms define a local model of the direction of propagation of the fiber at the voxel level and then they use this model to build fiber tracts ,i.e., D trajectories of putative fibers. The most used model is the Di usion Tensor Imaging which seeds a voxel and then propagates a fiber tract along ), 'Imaging structural co-variance between human brain regions', the principal direction of water di usion (Fig. . ). A probabilistic version of this model exists in order to estimate the uncertainty relative to the fiber tract direction. DTI is limited since it cannot capture complex fiber entanglements like fiber branching or crossing. Models capturing multiple fiber orientations are, for instance, the multiple tensor model [ ], -ball imaging [ ], or the di usion spectrum imaging [ ]. Regions of grey matter are the most interesting to track since they are the most concentrated region in terms of neurons. This would make those regions interesting to seed. However, the poor water di usion in grey matter renders di cult the fiber tract in those areas and necessitates appropriate seeding strategies [ , ].

Once the fiber tracts have been generated between all pairs of nodes, it is possible to derive estimated measures of connectivity strength.

The most common structural measures consist in counting the number of streamlines between any given node and node to construct the adjacency matrix. Since it is not possible to recover the directionality of the fibers, the adjacency measures from fiber tract count is symmetric, i.e, the associate network is undirected. This measures postulates that there is a direct correspondence between the number of axonal fibers and the number of fiber tracts. However, the fiber tracts count is sensitive to many factors such as the presence of myelin or noise (head motion, physiological noise. . . .

. Functional connectivity

Functional connectivity estimator Once nodes have been defined, M/EEG and fMRI techniques provide brain functional data under the form of time series associated to each node. Estimating the functional connectivity consists in computing a measure of statistical dependency between each pair of nodes and from the time series. The magnitude of this dependency is the weight between nodes and . Functional connectivity estimators fall into several nested categories: linear or non-linear dependencies between the data. directed or undirected which informs if they return a di erent weight for the interaction from node to node compared to the interaction from node to node . Can be of pure statistical nature (Pearson correlation. . . ) or taking advantages of the oscillatory nature of the signal (amplitude, phase, power. . . ) or based on an information theory measure.

One of the simplest linear connectivity estimators is the Pearson correlation of the time series. It provides undirected weights between -1 and 1 where -1 means that the two nodes are perfectly anticorrelated and 1 that they are perfectly correlated. A common practice Other linear measures can be for instance coherence, the equivalent of Pearson correlation in the frequency domain that captures amplitude correlations between two signals at the same frequency but at di erent nodes [ ]. In the frequency domain, one can quantify the pure phase synchronization at a given frequency for a pair of nodes through the phase locking value [ ]. FC estimators like the imaginary coherence [ ], or the phase lag index [ ] can even take into account phase lag correlations which are thought to describe underlying partial synchronization or traveling wave phenomena at the cell population scale [ ]. Regarding directed FC estimators, Granger causality is the most widely used one and it is defined both in the frequency and time domain. In essence, it quantifies how much the past signal of node is more informative on the signal of node than the past signal of node alone [ ]. The quantification of non-linear measures can be made throught the information theoretical measures of mutual information [ ] or synchronization likelihood [ ]. An important requirement of almost all the FC estimator is the stationarity of the signals, i.e., signal properties like mean and variance do not change over the recorded period. To reach stationarity, a widely used method is the sliding window analysis. It consists in estimating the connectivity in a small window of the recorded time (e.g., on the first time points over ones) and then to slide the window and again estimate the connectivity (e.g, on the time points to ).

The choice of the FC estimator should be grounded in the scientific question that the model is meant to solve. In [ ], Authors make the parallel with measuring the level of acquaintance between people in a social network by either measuring the number of mails or the number of time people see each other. Yet, depending on the context, some people may exchange mail while never meet each other and viceversa. For instance, a directed measure would be more appropriate to capture the spread of the epileptic synchronization than and undirected one.

Thresholding and binarization

By construction, functional networks are complete ,i.e., each pair of nodes in the graph has an attributed weight. However, some pairs are likely to be induced by noise either coming from the parcellation scheme, head motion or volume conduction e ects or simply the number of nodes and the length of the signal [ ] among other sources. Even if some noise reduction can be performed a priori by appropriate FC estimators [ ] it is a very common procedure to filter the network by keeping only the most significant weights. A widely used procedure consists in fixing a threshold below which any link weight is cancelled. Fixing the threshold can rely on statistical methods ,e.g., by fixing a statistical threshold on a null distribution obtained with surrogate data. The choice of the threshold can also be driven by the subsequent network analysis since many network measures are As an example, the characteristic path length decreases as the density increases which can render di cult the comparison on di erent subjects or conditions if filtering methods do not lead to the same density for each considered graph. Furthermore, it is important to mention that network measures were conceived for networks where nodes were connected with a measure of physical distance and not for similarity measure. Therefore, we need to be cautious when interpreting the measures like the characteristic path length on similarity networks.

Since most network measures were validated in the context of unweighted network, the significant weights after filtering can be binarized ,i.e., their weights set to to obtain an unweighted network.

All in all, the FC network reconstruction and analysis is a combination of choices where it exists no gold standard. [ ] report an extensive literature that has shown that FC analysis is strongly influenced by these choices which makes tedious any claim of generalities on one side but allows claims of specificity, i.e. , analysis whose choices were made toward a specific scientific question, in the other side. A summary of the functional brain networks analysis pipeline is given in Fig. . . Nodes correspond to specific brain sites according to the used neuroimaging technique. Links are estimated by measuring the FC between the activity of brain nodes; this information is contained in a connectivity matrix. By means of filtering procedures, based on thresholds, only the most important links constitute the brain network. The topology of the brain network is quantified by di erent network measures that can be represented as numbers (e.g., the colored bars). These measures can be the input to statistical analysis to look for significant di erences between populations/conditions (e.g., red points correspond to brain network measures of diseased patients or tasks, blue points stand for healthy subjects or resting states.

. Characterization of brain networks

In regards of the diversity of brain networks reconstruction and experimental conditions, no unique characterization of the node degree seems possible to be found to characterize brain networks [ ]. Nonetheless, like a vast majority of real-world systems, brain structural and functional networks were found to exhibit a broad node degree distribution, i.e., that brain networks exhibit the presence of hubs. While some studies reported a power-law degree distribution [ , ] possibly corrected by an exponential decay at its tail [ ], others reported distribution with an exponential decay [ , ]. Remarkably, it was shown that brain networks are more resilient under attacks than equivalent random or scale-free networks which could be explained by the previously reported tail exponential decay [ ].

Functional and structural networks at di erent scales all exhibit a small-world behavior [ ]. It means that they exhibit both a qualitatively high average clustering and a low characteristic path length or equivalently a high global and local e ciency. This balance have been related to the neuroscience concept of balance between information integration and segregation. It states that the brain functioning is rendered possible by the treatment of the information through specialized units integrated through large and e cient communication pathways. However, it was shown that the concept of small-world should be invoke with caution especially for functional brain networks due to pitfalls such as their reconstruction [ ]. Focusing on segregation, the high clustering was shown to emerge from the presence of communities of densely interconnected brain regions which were shown to be part of the same functional system [ ]. Furthermore, modules were shown to exhibit a nested hierarchical structure ,i.e., modules into modules that are typical of information processing networks [ , ]. More recent researchs have acknowledged the role of weights and in particular of weak connections in the understanding of brain functioning [ ].

Regarding integration, hubs are the natural candidates to interrelate modules within the brain. Although, they can be identified as the nodes with the higher degree, other measures can identify them more directly as the linking nodes between communities (e.g, through betweenness centrality). In structural brain networks, medial and superior parietal cortex as well as regions from the orbitofrontal, superior frontal, and lateral prefrontal cortex were identified as hubs [ , ].This claim is corroborated by their diverse and widespread functional connectivity across many di erent tasks [ ]. Hubs regions are also associated with regions with higher white matter concentration [ ].

Hubs seem to be connected between each other forming a rich-club of densely connected regions. In structural brain networks, regions including the superior frontal cortex, superior parietal cortex, and the precuneus, in addition to several subcortical regions including the thalamus, hippocampus, and part of the basal ganglia were shown to form a rich-club [ ]. Furthermore, almost of shortest paths between regions non included in the rich core pass through it. This supports the vision of a core-periphery structure of the brain where The brain is spatially embedded in space which induces constraints on the number of structural connections but also energetical cost to build long pathways across regions and maintain the electrochemical gradients for electrical signaling along them [ ]. Actually, modules of structural brain networks are often made of nodes that share a physical proximity and therefore short links [ ]. However, the brain exhibit also e cient integration properties that is characterized in network analysis in terms of high global e ciency. For instance, higher global e ciency of structural and functional network were associated with higher IQ [ , ] supporting the need of long-range pathways for e cient executive function. Therefore, "brain networks negotiate an economical trade-o between minimizing wiring cost and maximizing expensive but advantageous topological properties such as e ciency"[ ]. Based on this principle, a simple generative model of the human functional brain networks was able to successfully reproduce properties of empirical functional brain networks [ ].

Brain network changes in connectivity were related to brain diseases and their severity. The common trait of any brain diseases is a "dysconnectivity syndrome" ,i.e., an alteration of the network properties of the brain compared to healthy ones. These alterations can be spatially distributed as in Alzheimer's disease, or localized, as in stroke or traumatic injuries [ ]. Globally, brain diseases are characterized by a loss of node centrality, a less optimal small-world organization, related to loss of global e ciency and detrimental community structure reorganization [ ]. The quantification of the changes of network properties could represent an advance to monitor the disease and select new protocols to cure diseases.

. Network characterization of Alzheimer's disease Alzheimer's disease (AD) is the most prevalent type of dementia and a neurodegenerative condition. It manifests with initial mild memory impairments that progressively worsen, resulting in severe cognitive ), 'Functional neural network analysis in frontotemporal dementia and Alzheimer's disease using EEG and graph theory', ), 'Activity dependent degeneration explains hub vulnerability in Alzheimer's disease', deficits and ultimately leading to death. As of , approximately million individuals worldwide were a ected by AD and other dementias, and this number is expected to rise due to increased life expectancy [ ]. At the cellular level, AD is characterized by the gradual accumulation of tangles and -amyloid plaques, which contribute to the degeneration and loss of neurons and synapses. Consequently, brain atrophy occurs from subcortical regions to primary sensory areas, passing through association areas (Fig. . ). To characterize AD using network analysis, structural or functional networks are statistically compared with healthy controls (HC) that ideally match the AD group (in terms of age, gender, . . . ). AD-related changes in connectivity is a conundrum since studies alternatively report significant decrease or increase of functional connectivity in AD patients compared to HC [ ]. However, decrease are more often reported in M/EEG, MRI based networks while increase are reported for fMRI ones which suggests the major role played by the type of studied connectivity [ ]. In most of functional studies, AD patients undergo a loss of clustering coefficient [ -] Similar to connectivity, results on the characteristic path length are contradictory showing alternatively an increase or a decrease. Nonetheless, all studies report that AD brain networks are small-world, thought a majority report that AD brain networks are closer to equivalent random networks than HC ones [ -]. Linking behavior to observed networks, a study reported a positive correlation between lower characteristic path length and cognitive abilities [ ].

The more consensual result on AD-related e ects on brain networks seems to be the decrease of centrality of hubs compared to HC, in particular in association areas, such as the temporal lobe, medial parietal, posterior and anterior cingulate, and medial frontal areas [ -] which are areas known to undergo atrophy in the disease. [ ] show that the propensity of these hubs to form a rich-club is also a ected in the M/EEG alpha band. Furthermore, loss of centrality was positively associated with cognitive declined [ ]. From these considerations, [ ] show that selectively attacking hubs in HC brain networks make them resemble to AD brain networks in regard of their network properties. These results support the hypothesis of activity dependent-degeneration in AD [ ] stating that the more active regions are the most a ected by AD lesions. Altogether, the characterization of brain networks in AD seems dependent on the recording modalities and the various choice of functional brain networks reconstruction options. More recent studies attempt to benchmark protocols to maximize the di erence between AD patients and HC brain networks [ , ].Also, recent research endeavors are made towards more homogeneous patients cohorts and better match with HC [ , ]. Alternatively, AD-related brain networks studies, evolve towards integration of temporal or multiscale higher order structure.

Multilayer Brain Networks

. Multiscale brain modeling

The goal of multiscale modeling is to describe system's behavior by simultaneously considering multiple features, or mechanisms, taking place and interacting on di erent levels of information. Since most real-life phenomena involve a broad range of spatial or temporal scales, as well as the interaction between di erent processes, multiscale modeling has been widely adopted in several disciplines, from material science and algorithmics, to biology and engineering [ ].

In neuroscience, multiscale modeling has historically considered multiple levels ranging from microscopic single neuron activity to macroscopic behavior of collective dynamics. This is achieved by bridging biophysical mechanistic models of neuron dynamics and experimental neuroimaging data [ ]. Such "bottom-up" approach allows to predict macroscopic observables by integrating information at smaller scales, typically under the assumption of mean-field approximations [ -]. This means that neuronal ensembles' dynamics are progressively "averaged" across scales leading to a characteristic hierarchical nested structure where multiple units at finer-grained levels map into a new entity at coarser-grained ones (Fig. With the advent of new technologies and tools that allow gathering more precise experimental data and e cient processing, multiscale brain modeling has witnessed a significant transformation in the last decade. Increasingly more sophisticated and accurate models have been proposed including, among others, large-scale anatomical and functional brain connectivity (see Sec. . . ) [ , ]. However, to fully understand a multiscale system, models at di erent scales must be coupled together to produce integrated models across multiple levels. Indeed, global brain dynamics are strongly dependent on the interaction of several interconnected subnetworks that di erently contribute to generate them. Thus, the study of how intra-scale and inter-scale The use of "top-down" approaches, which start with the observation of biological characteristics in the intact system and then construct theories that would explain the observed behaviors, o ers complementary solutions. In particular, data-driven methods based on statistical signal/image processing of neuroimaging data allow to infer network representations of the brain at both anatomical and functional levels. Thus, while multiscale modeling in neuroscience has historically had a strong spatial connotation, it nowadays spans disparate levels of information, from structure/function to multiple oscillatory regimes and temporal evolution. Top-down approaches can be therefore used to generate richer and more realistic models reproducing real brain connectivity schemes and not just simulated ones [ ].

However, richer information and more accurate models also mean higher complexity and harder interpretation. These are both typical characteristics of multiscale problems that require the use of e cient algorithms to simulate the fully integrated model and appropriate In all this turmoil, questions like: how to model within-level and between-level relationships, how to characterize the resulting higherorder network properties, what are the critical phenomena emerging from the interaction of multiple levels, appear to be essential for advancing multiscale models. These questions and associated notions motivate the construction of a theory that explicitly builds on the capability to simultaneously characterize intralayer and interlayer connectivity.

In the next section, we will introduce the methodological framework of multilayer network theory, which is at the basis of recent developments in multiscale modeling of neural functioning.

. Common types of multilayer brain networks

. . Multiplex brain networks

Up to date multilayer brain networks have been mostly derived from experimental neuroimaging data in humans, with nodes representing the same entities, i.e. , brain areas across layers. Multiplex networks (see Sec. . . ) represent the easiest way to bridge brain connectivity at di erent levels, as one does not have to explicitly infer interlayer connections. In this situation, interlayer links only virtually connect the replica nodes and the associated meaning is basically the one of identity between the same nodes across layers (Fig.

. b) [ ].
This type of representation has been largely used to describe multimodal brain networks, whose di erent layers may contain structural and functional connectivity [ , , ], as well as interactions at different signal frequencies [ , , ]. A common situation when dealing with multimodal networks is that the nodes might not correspond to the same entity in their native space. This is for example the case of brain networks derived from fMRI and EEG signals, where nodes correspond respectively to image voxels and scalp sensors . To overcome this issue, advanced image and signal processing tools are used beforehand for projecting the native signals into the nodes of a common anatomical brain space, typically extracted from the structural MRIs of a subject's head

. . Temporal brain networks

Multiplex networks have been also adopted to describe temporal brain networks, i.e. networks whose topology is changing over time (see Sec. . . ) [ , , ]. In this case, each layer corresponds to a specific point, or instance, in time and only the replica nodes of temporally adjacent layers are interconnected according to a Markovian rule (Fig. . c). Unlike multimodal brain networks, the layers of a time-varying brain network do not correspond to di erent spatial or temporal/frequency scales, but they typically capture the dynamic network evolution within a fixed time resolution. This is typically in the order of milliseconds for motor behavior, minutes/hours for human learning, or years for aging as well as for neurodegenerative diseases.

. . Full multilayer brain networks

Full multilayer network representations, containing both intra-layer and inter-layer nontrivial connectivity, have been mostly adopted to characterize brain signal interactions within and between di erent oscillation frequencies (Fig. . Particular aspects of multilayer brain networks . . Multilayer brain networks are more than the sum of their layers

Multilayer networks give richer description than standard network approaches, but do they really represent a step forward into the modeling of brain organization? Why aggregating layers is not enough? Are all layers necessary to capture the main organizational properties? [ ] addressed these questions by introducing a structural reducibility approach to maximize the quantity of non-redundant topological information between the layers of a multiplex network with respect to its aggregated counterpart (Fig. . a). For a large spectrum of networks, from protein-protein interactions to social networks, structural reducibility showed that the best configuration in terms of distinguishability is not necessarily the one with the highest number of layers [

]. On the contrary, [ ] showed that multifrequency brain networks derived from fMRI signals were not easily reducible since all the layers brought some non-redundant topological information

. . Filtering spurious links in multilayer brain networks

It's important to remind that brain connectivity networks are estimated from experimental data. This necessarily implies the presence of spurious connections, often among the weakest ones, due to the statistical uncertainty associated with the connectivity estimator and/or the presence of signal artifacts during the experiment [ , ]. As a result of the construction process, multilayer brain networks are also influenced by such noise, which might alter the true association between the multiscale brain network organizational properties and the subject's characteristics and behavior.

To mitigate the presence of unwanted alterations in the estimated links, two main strategies have been so far adopted following what has been done in standard network analysis. The first approach consists in manipulating the brain signals, while the second one operates directly on the connectivity matrices. [ ] used the first approach to silence the e ects of head motion on the recorded brain signals and in turn on the estimated brain network. They tested di erent signal denoising strategies, mainly based on regression and source separation techniques [ ] on temporal brain networks constructed from fMRI data. Specifically, they evaluated their ability in attenuating the nuisance e ects on several network metrics, such as multiplex modularity and node flexibility (see Sec. . . ). Despite some variability, the obtained results suggested that regression-based approaches outperform source separation-based techniques, possibly due to their ability to explicitly incorporate the nuisance variables in the denoising process [

].

The second approach consists in filtering the network's links. This is typically achieved by fixing a threshold either on the percentage of the strongest links to retain or on their weights. Depending on the threshold value the resulting networks might have di erent densities and/or intensity.

[ ] evaluated the impact of network filtering on several topological properties such as multiplex PageRank (see Sec. . . ), multiplex modularity (Eq. . ) and participation coe cient (Eq. . ). Specifically, they tested several filtering criteria, e.g., spanning tree (MST) [ ], e ciency cost optimization (ECO) [ ], singular value decomposition (SVD) [ ] applied to each single layer separately, or adapted to the whole multiplex.

By using both synthetic and neuroimaging-derived multiplex networks, results indicated that SVD techniques lead to multilayer network properties that are quite robust to changes in connection density/intensity. MST and ECO techniques were instead e ective only when filtering each layer separately, and therefore useful when dealing with multimodal brain networks, where layers are estimated from di erent type of data and the nature of the interlayer links cannot be straightforwardly established.

Note however, that these results have been obtained for multiplexes and the extension to full multilayer networks still remains to be investigated. anatomical connection is likely to induce a synchronized activity between the corresponding brain regions [ ]. However, other significant configurations were reported including the presence of triangles in the functional layer with no support in the structural one. Overall these results indicated that intrinsic functional organization of the brain is non-trivially constrained by the underlying anatomical network [ ] and cannot be solely explained by it. Down the line, [ ] investigated the multilayer modularity of DTI-fMRI multiplex networks. Main results showed that the structural layer is mostly dominating the community structure of the multiplex over a broad range of topological scales explored by varying the granularity parameter (Eq. . ). Notably, the communities of the structural layer tended to spatially overlap with the cytoarchitectonic brain organization and were highly consistent across individuals. Instead, the communities of the functional layer were more heterogeneously distributed and less consistent across subjects, reflecting the dynamic repertoire of the brain functions [ , ].

By looking at DTI-fMRI multiplex networks, [ ] measured to what extent nodes with similar overlapping degrees tended to "wire" together, a property often referred as to assortativity. Results indicated that multimodal brain networks have a propensity to be assortative, which translates into an overall ability to facilitate system dynamics and resilience to random attacks (e.g., node removal) [ ]. This evidence resolved the assortative/disassortative dichotomy previously observed with single-layer analysis of structural/functional brain networks. Notably, such multilayer assortativity resulted from a nontrivial structure/function interplay and pointed out a novel organizational mechanism optimally balancing the resilience to damages and restrainability of their e ects.

Modeling the emergence of large-scale brain dynamics from microscale neuronal interactions is crucial for a mechanistic understanding of neural multiscale organization. An early study by [ ] proposed a computational model based on the structural connectome of the cat cortex. By parametrizing the coupling between several Fitzhugh-Nagumo oscillators according to the available connectome, they simulated the ongoing activity in each region, and estimated the interareal functional connections via Pearson's correlation [ ]. By means of this simple model, the Authors showed that a weak coupling parameter was su cient to generate biologically plausible macroscale activity, with functional connectivity patterns mostly overlapping the modular organization of the structural network.

[ ] used a similar approach based on the structural connectome of a macaque cortex and Wilson-Cowan neuronal models [ ]. More relevant to this chapter, they analyzed the behavior of multiplex clustering patterns (such as in Eq. . ) in the structural-functional networks as a function of two model parameters, i.e. , one tuning the input to excitatory neurons, and the other one modulating the input to the inhibitory ones. Specifically, they defined multiplex clustering indices to quantify the presence of functional links associated with common drivers in the structural layer. Main results showed that such quantities were maximal at the boundaries of the phase transition, from steady-state

b).

Di erently from previous results on single-layer analysis, this nontrivial behavior suggested that the system criticality does not only depend on the structure-functional interplay of the brain network, but also on the type of ongoing dynamics.

At the level of single neuron, [ ] proposed a multiplex approach to represent synaptic connections (structural) as well as extrasynaptic signaling interactions (functional) inferred from gene expression data of the C. Elegans worm. Despite the low degree of overlap between the synaptic and extrasynaptic connectomes, Authors found highly significant multiplex motifs (similar to the ones in Sec. . . ), pinpointing locations in the network where aminergic and neuropeptide signalling modulate synaptic activity. The presence of directed monoamine interactions and reciprocal synaptic connections was particularly significant among specific neurons implicated in learning, memory and motor functions. These results support the evidence that the structural/functional interplay is crucial to better understand the communication pathways between di erent parts of the C. Elegans nervous system.

In this direction, [ ] identified the shortest paths from touch sensory neurons to motor neurons allowing information flowing across di erent type of neurotransmitters and neuropeptides layers. By applying a time-delayed feedback control on the identified neurons, the Authors could eventually reproduce the typical C. Elegans locomotion, and characterize the neuromuscular multilayer connectivity mechanisms associated with the central pattern generator [ , ].

Multilayer network theory has just started to provide new tools and insights into the complex interplay of the brain structure and function.

Several issues remain to be explored such as how to establish interlayer connections [ ] or incorporate multilayer network mechanisms in the laws modeling the large-scale neuronal dynamics [ ].

.

. Information segregation and integration

Clustering and shortest paths are general concepts in complex systems that are both essential for e cient organization of many real-world networks [ , ]. These concepts reconcile two long-standing opposed views of the brain functioning. On one hand, phrenology-based theories, which associated di erent cognitive tasks with segregated brain regions [ ]. On the other hand, global workspace theories, which instead hypothesize the necessity of interareal integration of information to realize the very same tasks [ ]. Notably, network science has provided the tools to quantify network segregation and integration by demonstrating respectively the presence of many clustered connections and few shortest paths between areas. More recently, integration in the brain has been revisited and hypothesized to be determined by the presence of few core hubs in the network, and not directly by shortest paths [ , ]. By considering multilayer brain networks, segregation and integration become a joint property of both nodes and layers, thus providing information about higher-order phenomena such as cross-frequency coupling [ ], multimodal information [ ] [ ] investigated information segregation and integration in MEG full multifrequency brain networks. They first observed the presence of strong dependencies between intra-and interlayer connectivity. By decomposing the multilayer networks into representative connectivity structures, or "eigenmodes", they demonstrated that the overall amount of interlayer connectivity was associated with the second eigenmode, containing specific fronto-occipital network components common to all frequencies. In addition, they compared the empirical MEG multifrequency networks with those obtained from large-scale signals simulated with a thalamo-cortical model [ , ]. By increasing the model structural coupling parameter, the Authors reported a progressive increase in the resulting functional interlayer connectivity. Notably, real MEG multilayer networks maximally fit the model at the transition point of such increment, suggesting an optimal balance between segregation and integration of information between di erent frequency bands.

As for multimodal connectivity, [ ] investigated the associated integration properties by evaluating the core-periphery structure of DTI-fMRI multiplex networks. They specifically calculated the multiplex coreness (see Sec. . . ), which integrates information from di erent layers and provide a possibly more accurate characterization of the mesoscale brain network properties. Compared to single-layer analysis, results identified new core areas in the sensorimotor region of the brain that are key components of the so-called default mode network (DMN), i.e. , a set of brain regions that is active when a person is not focused on the outside world [ ]. Besides, results excluded previously established areas in the frontal region, whose belonging to the core system was still debated [ ]. By including structural (DTI) and functional (fMRI) network information, these findings o ered a

On longer time scales, [

] investigated how brain segregation changes with age by using longitudinal fMRI data acquired over a 4 years timespan. By computing the multiplex modularity (Eq. . ), they showed that the global flexibility, i.e. , the average node flexibility, is significantly higher in healthy elderly as compared to a temporal null model, where the brain network layers are randomly shu ed [ , ]. Results also demonstrated that people with more segregated temporal networks tended to be more resistant to transient changes in modular allegiance [ -]. Notably, older age was related to higher temporal variability in modular organization. However, no correlations were found with cognitive behavior, such as processing speed and memory encoding. Since flexibility is in general a good predictor of Taken together, these findings provided some concrete examples on how concepts such as segregation/integration of information can be broaden to capture multilayer brain mechanisms and provide complementary information about the system's behavior. While most of the studies have focused on undirected connectivity, future research will be crucial to include directed links and better inform on communication pathways in neuronal systems [ ].

. . Brain organizational properties of human behavior

The results presented in the previous paragraphs aimed to quantify the intrinsic structural and functional brain organization, with no reference to any specific mental state or behavior. Nonetheless, the brain is an extremely flexible and adaptive system, capable of altering its organization depending on endogenous and exogenous stimuli coming from the external environment (a property often referred to as plasticity). In this paragraph, we present some of the most recent results showing how multilayer brain network properties change according to specific behaviors, and how those higher-order topological changes are associated with inter-subject variability. Main results showed that, with appropriate model parameter values, the di usive interaction between the two layers prevents runaway growth of synaptic strength, both during ongoing activity and during learning. These findings suggested a previously unappreciated role for fast dynamic glial transport of metabolites in the feedback control stabilization of slow neural network dynamics during learning (Fig. .

b).

Notice that this is so far one of the few examples where multilayer network theory is used to model microscale neural organization across multiple temporal scales.

At larger spatial scales, [ ] used a multilayer network approach to characterize human learning during a simple motor task. In particular, they built temporal brain networks from fMRI signals across consecutive experimental sessions. They used the multilplex modularity (Eq. . ) to find long-lasting modules and found that community organization changed smoothly with time, displaying coherent temporal dependence, as in complex long-memory dynamical systems [ ]. Results also showed that the network flexibility changed during learningfirst increasing and then decreasing-demonstrating a meaningful biological process. In particular, the nodal flexibility (see Sec. . . ) was stronger in frontal, posterior parietal and occipital regions. Also, it predicted the relative amount of learning from one session to the following one (Fig. . ). These predictions could not be obtained via conventional task-related fMRI activation or standard network analysis, and confirmed the relation between network flexibility and cognitive performance. Indeed, network flexibility has been found to correlate with several mental states, such as working memory and planning [ , ], but also with mental fatigue [ ] and sleep deprivation [ ]. At this stage, it would be interesting to elucidate whether network flexibility is an aspecific predictor of cognitive performance or it can also distinguish between di erent dynamic brain states.

[ ] further study the cognitive load during attentional tasks in a EEG frequency-based multiplex framework. Based on betweenness centrality (see Sec. . . ), they observed an outflow of shortest paths from low frequencies toward high frequencies in the fronto-parietal regions. These findings suggest that cross-frequency integration of information is not only an intrinsic characteristic of the brain functioning [ ], but it is also modulated by attentional tasks as well as by drowsiness [ ].

In a recent study, [ ] investigated how the brain supports expressive language function by looking at MEG multifrequency brain networks. In particular, they aimed to identify the brain regions that are important for successful execution of expressive language in typically developing adolescents. To this end, Authors first identified the multifrequency hubs by means of a modified version of the multilayer PageRank centrality and then reranked them according to their importance in fostering interlayer communication. Compared to standard singlelayer analysis, this two-step procedure allowed to capture nonlinear interactions and resolve the task-related brain areas with a higher spatial resolution. These regions mostly lied in the left hemisphere and represented possible conduits for interfrequency communication between action and perception systems that are crucial for language expression [ ].

Planning and executing motor acts is accompanied by changes in brain activity and connectivity on very short time scales of the order of milliseconds [ , ]. [ ] used an EEG temporal network approach to characterize such fast brain functional organization during a simple foot movement task. Compared to network sequences with randomly shu ed layers, brain networks showed a higher temporal clustering and a similar characteristic temporal path length (see Sec. . . ). Put di erently, dynamic brain networks exhibited a temporal small-world propensity, supporting both segregation and integration of information through time. While single-layer analysis had previously unveiled that segregation/integration properties fluctuate and adapt over the di erent phases of the movement [ ], these findings provided new evidence on the intrinsic global temporal properties of motor-related brain networks.

. Multilayer network-based biomarkers of brain disease

Like any other complex system the brain can exhibit anomalous connectivity, which in turn may lead to abnormal behavior and clinical ), 'Abnormal neural oscillations and synchrony in schizophrenia', symptoms. Those brain connectivity changes can be spatially distributed, such as in schizophrenia or Alzheimer's disease, or localized such as in stroke or traumatic injuries [ ]. Looking at the network organization in both healthy and diseased conditions appears therefore fundamental to understand the resilience and vulnerabilities of the brain [ ]. From a medicine perspective, network-based biomarkers would represent advanced tools to monitor the disease progression and inform new therapeutics to mitigate or counteract the e ects of the disease. In the last decade, standard network analysis has accumulated evidence documenting general reorganizational properties such as departure from optimal small-world configurations, aberrant modular reorganization, as well as significant loss of node centrality [ ]. Thus far these network changes have remained associated with a particular aspect, or layer, of information. Since brain pathologies typically result from multifactor processes at di erent scales and levels, multilayer brain networks naturally constitute a more appropriate integrative modeling approach. In the following, we present some of the most recent results obtained for di erent brain diseases, that provide new perspectives on the impacted multiscale network properties and can be used to improve diagnosis and prediction.

. . Neuropsychiatric disorders Among neuropsychiatric disorders, schizophrenia is certainly one of the most studied ones due to its large population incidence. In , over million people su ered from schizophrenia worlwide [ ]. Typical clinical symptoms include hallucinations, emotional blunting and disorganized speech and thoughts. The biological causes of schizophrenia are still poorly understood and many hypotheses are currently being investigated based for example on neurotransmitter dysregulation [ ], myelin reduction [ ] as well as oxidative stress [ ]. At large spatial scales, low and high frequency neuronal oscillations, as well as their interactions, have been widely documented as a core feature of the neuropathology underlying schizophrenia [ ]. Functional connectivity changes, within and between frequency bands, have been reported in schizophrenic patients [ ] and associated with persistent symptoms leading to disorganization of visuomotor mental functions [ ].

By using a multiplex approach, [ ] provided a first integrated characterization of the topological changes in schizophrenia from resting state fMRI-derived multifrequency networks. In particular, they evaluated the multiplex PageRank centrality (cf. Sec. . . ) and showed a substantial reorganization of the most important multifrequency hubs of the brain, such as the precuneus cortex, a key region for the basic physiological brain organization [ ]. When injected into a random forest classifier, multiplex PageRank centrality metrics led to a classification accuracy of 80%, which is higher than standard network approaches, but comparable with otherwise much more sophisticated machine learning techniques. At cellular levels, schizophrenia has been hypothesized to result from excitatory-inhibitory neuronal dysfunction, with a consequent abnormal temporal coordination between large-scale macro areas of the cerebral cortex [ , ]. By investigat- From a pure classification perspective multilayer brain networks have been also used as alternative multidimensional features to better discriminate between schizophrenic and healthy subjects. [ ] considered a working memory fMRI experiment and built a -layers multiplex brain network where each layer contained a di erent type of nonlinear functional connectivity. For each layer they extracted standard nodal centrality metrics (i.e., strength, betweenness, clustering, and PageRank) and used them as classification features. Compared to a single-layer networks, built from simple linear correlations, they ), 'Removing high-frequency oscillations', achieved a significantly higher classification (≈ 90% vs. ≈ 70%) for di erent types of working memory tasks. Following the same goal, [ ] considered resting state fMRI data over a group of healthy individuals and a group of patients with schizophrenia. Originally, they built for the two groups a multiplex brain network, where each layer represents the functional network of a specific individual. By extending the popular node vec unsupervised network embedding procedure [ ], they learned continuous node feature representations from multilayer networks based on random walkers which are allowed to move across layers. The resulting embeddings revealed a higher variability for the similarity between the nodes in the default mode network and salience subnetwork, suggesting a less stable within-module brain organization in the schizophrenic group. While the overall classification accuracy did not outperform state-of-the-art performance, learning the features in an unsupervised approach might be nevertheless important for future applications in automatic diagnosis.

Major depressive disorder (MDD) is clinically characterized by severe fatigue, aphasia, di culty to focus and suicidal thoughts in extreme cases. Symptoms are diverse and their severity largely di ers among patients. Since e ective treatments are currently available, scientific research mostly focuses on identifying predictive biomarkers to enable a more personalized therapeutics. Previous studies suggested that MDD leads to several brain signal alterations a ecting functional connectivity within but also between di erent frequency bands [ , ]. To fully exploit this multifrequency information, [ ] proposed a full multilayer approach to improve the diagnosis of MDD. Specifically, they developed a convolutional neural network that directly takes as input the full multilayer brain networks to learn and extract the most discriminant features. The resulting classification accuracy (≈ 97%) was comparable to state-of-the-art methods based on specific frequency bands. While promising, these findings suggested that machine learning algorithms for multilayer brain networks still remain to be finetuned in view of their concrete implication in the identification of the best intervention strategy to cure or alleviate MDD-related symptoms.

. . Other neurological diseases

Epilepsy is a group of neurological disorders characterized by seizures, which may vary in time and intensity, from short mild awareness loss to long vigorous convulsions. Epileptic seizures are underlied by excessive synchronized neuronal activity in the entire cerebral cortex or in parts of it. In , about 27 million people su ered from epilepsy [ ] among which 30% are not curable with drug treatments [ ]. Clinical research mostly aims at identifying predictive neural markers of the seizures to allow preventive treatments or to localize the origin of the seizure to inform precise surgery [ ].

Recent evidence has showed that epilepsy seizures are characterized by brain functional connectivity changes within, but also between, di erent brain signal frequencies [ -]. From a topological per- ), 'The JFK Coma Recovery Scale-Revised: Measurement characteristics and diagnostic utility', spective, decrements of network e ciency have been reported between low-high frequency bands, before the seizure onset, and were associated to sensorial disturbance and mild loss of consciousness [ ].

The intrinsic relationship between structural and functional layers can also unveil hidden connectivity structures characterizing di erent types of epilepsy. In this direction, [ ] used a DTI-fMRI multiplex approach to classify between epileptic seizures originating in di erent zones of the brain, namely the frontal and temporal lobe. In particular, Authors extended the concept of multiplex motifs to include subgraphs with more than 3 nodes (cf. Sec. . . ). The most frequent multiplex patterns consisted of links from both structural and functional layers that were spatially localized. Notably, the structural components were quite stable across conditions and involved regions belonging to the DMN system ,i.e., cuneus, precuneus, and peripheral cortex [ ]. Instead, the functional counterparts of the multiplex patterns, were highly variable and mostly involved regions concentrated in the respective epileptogenic zones, i.e., temporal and frontal lobes. Eventually, Authors demonstrated the superiority of these multiplex connectivity patterns to discriminate between epileptic patients and healthy controls (72-82% classification accuracy), as compared to equivalent single-layer metrics or di erent multiplex metrics such as multiplex PageRank or algebraic connectivity (see Sec. . . ). These results are in line with the one-to-many relationships between structural and functional brain networks [ ], and can be used to finetune the research of predictive biomarkers in epilepsy.

Consciousness disorders regroup a variety of symptoms which go from complete loss of awareness and wakefulness, such as coma, to minimal or inconsistent awareness [ ]. The di erential diagnosis between the di erent types of disorders of consciousness is paramount to identify the best medical therapeutics. Recent results suggest that frequency-dependent functional brain connectivity is crucial to characterize impairments of consciousness, as well as to predict possible recovery processes [ -]. In an e ort to provide an unified picture on the role of brain connectivity within and between frequency bands, [ ], adopted a multilayer network approach. By investigating brain networks derived from source-reconstructed EEG signals, Authors aimed to distinguish between patients su ering from unresponsive wakefulness syndrome (UWS) and minimally conscious state (MCS), which often present similar symptoms [ ]. Results showed that several nodal multiplex metrics, including overlapping clustering, betweenness and multiplex participation coe cient (see Sec. . . ), were significantly lower in UWS as compared to MCS patients. This was particularly evident in the frontoparietal regions of the brain whose relative loss of multiplex centrality is associated with the behavioral responsiveness of the patients quantified by the coma recovery scale [ ]. By adopting a full multilayer network approach, Authors eventually reported a significantly lower interlayer connection intensity in the UWS group and could spot out those patients who regained consciousness one year after the experiment. Notably, the discrimination between UWS and MCS patients was not observed when looking separately at frequency-specific network layers. Although very pre- While the consequences of these changes on large-scale brain networks have been widely investigated, the accumulated results are often discordant and depend on the considered spatial or temporal scale [ , ]. Multilayer networks represent an interesting approach to get an integrated, potentially more informative picture of the disease. Of note, these losses of functional hubs could not be observed when looking at individual frequency layers. Such multilayer hub disruptions correlated with the accumulation of -amyloid plaques in the cerebrospinal fluid, but also with the cognitive impairment of patients, demonstrating a potential clinical relevance. By using the multiplex participation coe cient (Eq. . ), results indicated that most vulnerable hub regions in patients with AD also lost their ability to foster communication across frequencies compared to healthy control subjects. Similar results were obtained independently by [ ], showing a significant loss of multifrequency hubs in DMN regions and a strong association with memory impairment. By using a classification analysis, they eventually showed that integrating multiparticipation coe cient values with equivalent single-layer network metrics leads to improved distinguishability of AD and healthy subjects. More recently, [ ] showed that AD patients could be predicted by the lower values of algebraic connectivity 2 in resting-state MEG multifrequency networks. These results spot out new network mechanisms that hinder information load from flowing through di erent frequency bands and eventually impair the cognitive abilities of AD patients.

[

], addressed similar questions in EEG multifrequency brain networks. They showed that both multiplex clustering (Eq. . ) and multiparticipation coe cients (Eq. . ) presented significant decrements with respect to healthy controls in the posterior areas of the brain. These results confirmed a general tendency in AD patients to loose segregation and integration of information across signal frequencies. Yet, few observed increases in frontal areas suggested the presence of some compensatory mechanisms to be further elucidated [ ]. In the same study, Authors also investigated the dynamic aspects of EEG brain networks in AD from a purely temporal perspective (Fig. . a). By using the aforementioned multilayer metrics, they showed that AD temporal segregation was mostly impacted by AD in frontal and occipital areas, while temporal integration properties were less a ected as compared to healthy subjects, mainly because of its higher variability across nodes. However, when combined together, nodal values of temporal segregation and integration led to a very high discrimination between AD and healthy subjects (> 90% accuracy), suggesting that spatial heterogeneity of temporal integration may also be related to progression of the disease (Fig. . c).

To integrate and disentagle the role of di erent neuroimaging modalities in AD, [ ] built multiplex networks composed of di erent connectivity types derived from DWI, fMRI and MEG data. This represents so far the most complete type of multiplex brain network merging together structural and functional information (Fig. . a). By focusing on the mesoscale properties (see Sec. . . ), Authors showed a selective reduction of multiplex coreness in the AD population, mainly involving temporal and parietal hub nodes of the DMN that are typically impacted by the anatomical atrophy and -amyloid plaque deposition [ ]. Such significant loss was mainly driven by few layers notably DWI, fMRI and MEG in the alpha (7 -10 Hz) frequency range, and could be explained by a simple model reproducing the progressive random disconnection of the multilayer network via the preferential attacks of its core hubs (Fig. . b). From a clinical perspective, Authors eventually reported that patients with larger coreness disruption tended to have more severe memory and cognitive impairments, in line with the general tendency observed in other previously described studies [ ] (Fig. .

c). Recently, [

] built -layer multimodal networks from gray matter atrophy and amyloid deposition across di erent stages of AD in humans. Within a rigorous controlled study, they provided very specific results, not obtainable with traditional approaches from single imaging modalities. Notably, multiplex modularity (Eq. . ) revealed a characteristic module in the temporal brain area, likely to reflect the transition to AD dementia. Decreased values of multiplex participation coe cients (Eq. . ) in atrophy-related hub regions were also found in the later AD stage as compared to healthy controls. This study shed light on the non-trivial interplay between -amyloid level and grey matter atrophy and its clinical relevance for AD.

Taken together, these results indicate that AD is characterized by a previously unappreciated multimodal and temporal dysconnection mechanism that primarly a ects regions impacted by the atrophy process. Future research will be crucial to elucidate whether such a disruption tendency is compensated by other multilayer mechanisms, possibly involving more intact cortical systems (e.g., sensory motor) [ , , ]. . Principles of node-layer duality

. . Primal and dual description of multilayer networks

Unlike single-layer networks, the basic interacting unit of a multilayer network is the duplet ( , ) which quasi systematically refers to node interacting across layer . In this agent-centric perspective, links are neutral and only account for the presence (and the intensity in case of weighted multilayer network) of an interaction between node-layer duplets, the latter carrying the nature of the interaction. Let us define , the supra-adjacency matrix of this system whose entry contains the connectivity between the node at layer and the node at layer [ ].

Since the ordering of the node-layer duplet ( , ) is arbitrary we can rewrite it as ( , ) if we carefully transpose all the related indices in our mathematical representation [ ]. Now, the duplet ( , ) refers to layer (considered here as the interacting agent ,i.e., 'nodes' ) interacting across node (considered here as interaction spaces ,i.e., 'layers' -Fig. . ). Let us define , the supra-adjacency matrix of this transposition whose entry contains the connectivity between the layer at node and the layer at node .

From now, we call the nodewise description the one where nodes are the units of interest (matrix ) which we refer to the symbol X. Similarly, we call the layerwise description, the one where layers are the units of interest (matrix ) which we refer to the symbol Y, in the rest of this manuscript.

At the node-layer level, and are trivially related by the formula = , which implies that = where is a permutation matrix. Therefore, the structural properties of the node-layer duplet are description-invariant ,e.g., a clustering coe cient of ( , ) is the same as the one of ( , ). In the next section, we explore where the relevance of this dual representation of multilayer networks is by exploring the di erent scales in a multilayer system. The unitary element of a multilayer network is the duplet ( , ) corresponding to the node and the layer . Because of the intrinsic node-layer isomorphism, the system can be equivalently represented by nodes plunged into layers or by layers plunged into nodes. Depending on the representation side, two complementary descriptions of the same system can be obtained. In the primal nodewise, connectivity is integrated across the layers and the topological properties of the nodes are provided. In the dual layerwise, connectivity is integrated across the nodes and the topological properties of the layers are obtained instead. In the networks, solid lines intralayer links, grey lines interlayer links, dashed lines replica links. Any structural measure associates a property (e.g., number of links or number of triangles) to the unit of interest. In the nodewise description, the units of interest could be the node-layers, the nodes, or the whole multilayer network, depending on the scale. In the layerwise description, layers have the role of nodes. A direct consequence of the framework is that, since layers are the interacting units in the layerwise description, we can conceptually attribute nodal measures on layers and viceversa. For instance, we can conceptually define the clustering of layers as well as less evident measures such as the structural reducibility of nodes [ ].

[ ]: Artime et al. ( ), Multilayer Network Science, : In real-world multilayer networks, the number of layers are small. Hence the characterization of the layer distribution is irrelevant in a pure statistical perspective. However, even for = 2, the standard deviation of the layer distribution is relevant in characterizing the system as we saw in the example (see also next section).

In the nodewise description, node properties are obtained by integrating the information of the node-layer level across the layer dimensions. For instance, the node multidegree centrality (node MC) is given by X = ,i.e., the sum of links adjacent to node across all layers [ ]. In the layerwise description, layer properties are obtained by integrating the information of the node-layer level across all nodes (Fig. . ). For instance, the layer multidegree centrality (layer MC) is given by Y = which is the sum of links adjacent to layer across all nodes.

Nodes and layers properties are, in general, not trivially related since connectivity changes a ecting one dimension are not necessarily visible in the other one. As an example, let us explore the first two moments of the node and layer MC distribution. By definition, X = 1 X

, where is the number of nodes and

Y = 1 Y
where is the number of layers. It implies that the first moment of the node and layer MC distributions are related by the formula: X = Y . Therefore, the first moment of the node and layer MC distribution is description invariant. In extenso, the properties of the whole multilayer networks are also representation invariant (e.g. , 2 = X = Y ) since information on the global scale is needed to compute the first moment.

However, the second moments or the standard deviation of the distributions are not trivially related (Fig. . ). For instance, if we have a system with = 2 layers and nodes, and we fix the node MC distribution, we can either imagine a system where the layer MC is evenly distributed across the two layers, i.e. , the second moment of the layer MC distribution is zero. Or, we can imagine a system with still the same node MC distribution but with 90% of links in one layer and 10% in the other one. It would lead to a large standard deviation in this case .

Therefore, studying the properties of nodes and layers at the 'second order' (referring here to the moment) is the level at which nodes and layers' properties emerge as the two degrees of freedom of multilayer systems, i.e. , the results become description-dependent (Fig. . ). Based on these considerations, we can build synthetic systems that have completely uncorrelated node and layer structural properties 'at the second moment'. Therefore, the potential relationship between these two descriptions' properties could give us precious insights on the nature of real-world systems. In the next section, we propose a "perturbation" approach in order to establish a formal relationship between both descriptions using node-layer duality.

. . Primal and dual descriptions retain complementary information on the system

As we saw in the previous section, we can think about two systems that share similar nodewise properties but completely di erent layerwise properties and viceversa (Fig. . ). In the meantime, multilayer systems are uniquely defined by the conjunction of nodewise and layerwise properties. For example, the multidegree centrality (MC) is associated with the number and position of links in a system. If we change one link by at least one layer indices while keeping the nodes indices fixed, the node MC X stays fixed while the layer MC Y changes.

Similarly, if we change a link by at least one of its node indices while keeping the layer indices fixed, Y stays fixed and X changes. Finally, any other changes that a ects at least one node and layer index will imply both X and Y to change. Therefore, any perturbation in a system is followed by a change of the MC sequence X and/or Y (and the second moment of their related distribution). In this perspective, X and Y uniquely characterize a multilayer system in terms of number and position of links. Indeed,fixing X and Y define a network ensemble whose elements have the same X and Y but which are di erent regarding another network property. Nonetheless, any perturbation of this other property can be track by either the nodewise and/or the layerwise description. As another example, the node and layer clustering coe cient sequences uniquely characterize the position and number of triangles in a multilayer system.

From these considerations, we would like to quantify the e ects of the local connection perturbation into the nodewise and layerwise representation. In particular we ask, what are the perturbations that are visible in the nodewise description only, in the layerwise description only ? For which perturbations one description better discriminates the initial and the perturbed networks ? Are they related by formal duality ? Does it depend on the initial system ?

Next, to answer all these questions, we elaborate a stochastic rewiring model which allows us to control the amount and the types of perturbations in a multilayer system. The normalized nodewise standard deviation is X = 0.73. The normalized layerwise standard deviation is Y = 0.45.b) Eu-Air multiplex network is rewired in a way that layerwise degree distribution stays fixed and the nodewise one is modified using the stochastic rewiring algorithm ( = 1, = 1, see Sec. . . ). The normalized standard deviations become X = 0.10, Y = 0.45.c) EuAir multiplex network is rewired in a way that nodewise degree distribution stays fixed and the layerwise one is modified using the stochastic rewiring algorithm ( = 0.66, = 1, see Sec.

. . ).The normalized standard deviations become X = 0.73, Y = 0.15. Kernel density estimation are obtained using the seaborn Python library.

Figure . : Schematic illustration of the di erent types of link rewiring.

The type of rewiring is determined by the combination of three parameters. Specifically, determines the probability to make a link displacement while keeping the layers unchanged.

determines the probability to make a link displacement while keeping the nodes unchanged.

determines the probability to make a link displacement while changing both nodes and layers. Dotted lines old position, solid lines new position.

+ + = 1.
. All the indices are drawn unformly at random. At least one node index and one layer index must be di erent from the one of the selected links.

For example, = 0.3, = 0.6 and = 0.1 implies that about 30% of the selected links will be rewired without altering the initial connected layers, 60% without altering the nodes, and 10% by altering both layers and nodes. In essence, the model selectively randomizes part of the node degree sequence, the layer degree sequence or both. Note that when = 1, the layer MC sequence is entirely preserved. Similarly, when = 1, the node MC sequence is conserved.

The algorithmic implementation of this model is available at the Appendix and in a GitHub repository and its computations are further compared with the model for validation of the latter (see Sec. . . ).

.

. Parameters for a complete uniform rewiring

To make a connection to existing reshu ing models, we show here the combination of parameters of the stochastic rewiring model that makes it a uniform rewiring of all the links within networks. We consider two cases:

Multiplex network: links are present within layers and are only allowed between replica nodes across layers. No self-loops are included.

Multilayer networks: links are allowed within layers, and between any nodes across layers. Self-loops are included.

The complete uniform rewiring reduces to find the probability parameters , , for which the e ective probability for a link to be rewired anywhere else is 1 -1 , i.e. , uniform across the network. Therefore, the probabilities depend on the size of the system with = -1 -1 and = -1 -1 with and , the total number of unique links in a (inter)layer and the total number of unique links for a given pair, respectively (see Table of Sec. . . ). The -1 term ensures that we do not rewire the link at its place.

For multilayer networks, the probability parameter for complete uniform rewiring read:

[ ]: Newman ( ), Networks: An In- troduction, = 2 -2 ( ) 2 -2 ( . ) = 2( 2 -1) ( ) 2 -2 ( . )
For >> 1, >> 1, ∝ -2 and ∝ -2 . Therefore, = 1 reduces to a complete uniform rewiring in the case of large multilayer networks.

For multiplex networks, the probability parameter for complete uniform rewiring read:

= ( -1) -2 ( -1) -2 ( . ) = 2( -1) ( -1) -2 ( . )
For >> 1, >> 1, ∝ -2 and ∝ -1 . Therefore, = 1 reduces to a complete uniform rewiring in the case of large multiplex networks. With most real-world multiplex networks having ∼ 1, has often a non-negligible contribution in the complete uniform rewiring.

. . Preliminaries and hypothesis

We make the hypothesis that multilinks are allowed in the rewiring process, i.e., links can be rewired where links exist already [ ]. This hypothesis ensures that the rewiring process is markovian and so, analytically tractable. In other words, it ensures that the transition probability does not depend on every previous rewiring step. We also impose the constraint that each link that is rewired cannot be rewired anymore in the network. This ensures the number of rewiring to be equal to the number of unique links in the system at the maximum. Furthermore, we reduce our analysis to undirected unweighted multilayer networks which exhibit the same set of nodes across each layer.

Our objective is to find ( ), the value of link connecting ( , ) to ( , ) averaged over all possible configurations of the rewiring model after rewiring. If we denote ( ( ) = ), the probability for the latter link to have the value at the -th rewiring, we can write ( ) = ( ( ) = ). Since we deal with undirected and unweighted networks (links have value or ) we obtain that ( ) = ( ( ) = 1). Since we are ultimately interested in the expected MC sequences, we will progressively integrate ( ) in the following derivation.

. . Π( ) does not depend on the rewiring step By construction, each link that is rewired at step cannot be rewired anymore at subsequent steps. Also, it means that the population of links to be rewired progressively decreases as the network is rewired. Furthermore, Π ( ), the probability of selecting the link connecting (a, ) and (b, ) at the -th rewiring depends only on what happens on the ( -1)th rewiring step: either the link was selected or it was not selected at a previous step. If it is selected, the probability of selecting a link is zero for all subsequent steps. Let us denote the conditional probability Π ( | -1) for any link to be selected at timestep knowing it was not at step -1. It implies that:

Π ( ) = Π ( | -1)×Π ( -1| -2)×.......×Π (2|1)Π (1), ( . )
where bars, represent the complementary probability of links not being selected at a given timestep. At any timestep , the link is drawn uniformly at random among the population of links that is still rewirable. If we suppose that the link connecting (a, ) to (b, ) is not rewired at time -1, it implies that:

Π ( | -1) = (0) 1 - 1 -( -1) = (0) - -+ 1 , ( . )
where (0) ensures that the binary link is present in the initial network before any rewiring. Since it is binary, we make it appear only once in related expressions. Using Eqs. . , . , it read:

Π ( ) = (0) -( -1) Π -1 =1 - -+ 1 = (0) -( -1) -( -1) = (0) ( . )
Therefore, the probability to select a given link stays constant for any rewiring step.

.

. Derivation of the expected multidegree centrality sequences

The node MC degree derivation is not shown since it follows the exact same canvas.

Following the notations, we derive the following rate equation for ( + 1) = ( ( + 1) = 1). Since we derive only expected values, we lighten the notations by not considering overlines in the rest of

We define (0) and (0) as the total number of unique links present in the initial multilayer network in the (inter)layer , and the total number of unique links present in the initial multilayer network in the pair i,j, respectively. Simplifying the previous equation leads to:

( + 1) - ( ) = (0) - (0)
-1 ( . )

+ (0) - (0) -1 + (( - (0) -(0)) -( -Δ + 1) (0)
-Δ

We note that the right-handside is a constant that depends only on the structure of the initial multilayer network and not on the rewiring step. Therefore ( ) follows a simple arithmetic sequence that we can explicit:

( ) = (0) + (0) - (0) -1 .
)

+ (0) - (0) -1 + ( -(0) - (0) -( -Δ + 1) (0)
-Δ

Since we are interested in knowing the expected node and layer MC of a network that underwent the rewiring process, we need to sum the above expression. Here, we show the derivation for the layer MC by summing over the node component i,j. Note that > 1 = .We set = with ∈ [0, 1]:

( ) = (0) + (0) - (0) -1 ( . ) + - (0) -1 + ( -1) - (0)( -) -Δ
The term that multiplies equals zero, which is expected by construction of the stochastic rewiring model. By summing over one of Node-layer duality of multilayer networks the layer components, we obtain the layer MC:

Y ( ) = Y (0) + - Y (0) -1 ( . ) + ( -1) -Y (0)( -) -Δ ,
where is a coe cient that depends on the multilayer type (multilayer or multiplex-see definitions in Sec. . . ). For the sake of readability, the derivation for the node MC is not shown since it follows the exact same canvas. The node MC read:

X ( ) = X (0) + - X (0) -1 ( . ) + ( -1) -X (0)( -) -Δ ,
where is a coe cient that depends on the multilayer type (multiplex, full multilayer).

Depending on the type of multilayer networks , , , have di erent values. We distinguish two main types of multilayer networks: , we have that Y = 2 . For multiplex networks, the sum reduces to one term so = 2 and for multilayer networks it reduces to terms so = 2 . counts the number of terms needed to go from to X . Since is the number of unique links from to , X = . For multiplex networks, the sum reduces to -1 (no self-loops) term so = -1 and for multilayer networks it reduces to terms so = .

.

. Fit between analytical derivation and numerical computations

First, we validate Eqs. . , . by numerically implementing them to ensure that they correspond to a rewiring process, i.e. , the sum of the total node and layer MC is conserved for every combination of , , . Then, the main hypothesis that supports the analytical derivation of the model is that multilinks are allowed in the rewiring process, i.e., links can be rewired where links exist already.

Nonetheless, real-world multilayer networks do not exhibit multilinks in general. It implies that the transition probability is zero for a link that exists already in a realistic rewiring. Consequently, the transition probability to rewire a link in any available position is comparatively higher in the realistic rewiring than in the analytical rewiring. Therefore, the higher is the density, the higher is the discrepancy between the model and the realistic rewiring. Since real-world networks are sparse (low density), the discrepancy should stay low in practice.

Indeed, we show that the analytical model and the realistic rewiring (implemented in the algorithmic version of the model, see Appendix ) match for sparse multilayer networks and sparse multiplex networks (Fig. . ). To do so we define MC-based Euclidean distances that compare the MC sequence of the initial multilayer network to the expected one for both types of rewiring.

In general, for two MC sequences 1 , 2 the Euclidean distances read:

( 1 X , 2 X ) = ( 1 ,X -2 ,X
) 2 ( . )

( 1 Y , 2 Y ) = ( 1 , Y -2 , Y
) 2 , ( . )

We will use these MC-based distances extensively in the rest of this manuscript.

. . Expression of Y ( ) and the distances in the thermodynamic limit

In the previous section, we obtained the expression of the expected node and layer MC after rewiring according to the model (see Eqs. . , . ). Note that these equations are valid for any and . In the limit of large networks ( >> 1 and >> 1), the previous formulae reduce to simpler, more interpretable expressions:

lim , →∞ ( ) = (0) + (1 - )( - (0)) lim , →∞ Y ( ) = Y (0) + (1 - )( Y -Y (0)) ( . )
where Y is the ensemble average of the layer MC. To obtain this result, we remark that = , by definition of and and also that 1 -= + . Similarly, the expression of the node MC in the thermodynamical limit read:

lim , →∞ X ( ) = X (0) + (1 - )( X -X (0)), ( .
)
where X is the ensemble average of the node MC.

These results are consistent because it implies that the MCs are driven by the di erence between the initial MC and the ensemble average layer MC. Therefore, the higher / and the higher , the higher the rewired network tend to be to an equivalent random multilayer network where each layer/node MC is close to Y and perfectly equals to it on average for / = 0 and = 1.

We finally compute the distance between the initial network and the expected layer MC. Distances are the euclidean distances between the components of the MC sequences of two multilayer networks (Eqs. . , . .) With the MC formula, distance between initial and expected perturbed network read:

( Y , Y ( )) = (1 - )( Y -Y (0)) 2 ( . ) ( Y , Y ( )) = (1 - ) √ ( Y ), ( . )
where ( Y ) is the standard deviation of the initial layer MC distribution. Similarly, the node MC distance read:

( X , X ( )) = (1 - ) √ ( X ), ( . )
where ( X ) is the standard deviation of the initial node MC distribution. Note that Eqs. . , . , . are valid for any type of multilayer network composed of the same set of nodes across layers. In general, both distances increase with the system size, the MC heterogeneity (i.e., the variance), as well as with the type and amount of changes. Specifically, X increases with (and/or ), while Y increases with (and/or ). A node-layer duality relationship then read:

Y = 1 - 1 - ( Y ) ( X ) X ( . )
. . Closed-form of the standard deviation of random multilayer networks

We need to derive the standard deviation of the node and layer MC for random multilayer networks to obtain a closed-form of the distance and the formal duality relationship. Since we deal with random multilayer networks, their number of links follows a binomial law. Let us define , , , the maximum number of interlayer, intralayer and replica links in the network, per interacting unit X, Y (node or layers), respectively. The probability of a unit to have intralayer links, interlayer links and read:

( = ) = (1 -) - ( = ) = (1 -) - ( . ) ( = ) = (1 -) - ,
where is the probability parameter. The contribution of each link to the MCs depends on its type.

Actually, each intralayer links contributes to two units of the layer MC, Y = 2 and each interlayer link and each replica link contribute to one unit of the layer MC, Y = , Y = (Fig. . ). Each replica link contributes to two units of the node MC, X = 2 and intralayer and interlayer links contribute both to one unit for the node MC, X = and X = . The following table makes explicit , and :

X ( -1) 2 ( -1)( -1) 2 ( -1) 2 Y ( -1) 2 2 ( -1)( -1) 2 ( -1)
Since the number of links follows a binomial law, the MC variance has the typical form 2 ( ) = (1 -), with the parameter of the random model.

Figure . : E ect of intralayer link additions on node and layer multidegree centralities.

In general, adding one intralayer link in the nodewise description increases the multidegree centrality in the layerwise description by a factor of (right side), while the node multidegree centrality only increases by a factor of one (left side). Because in multiplex networks only changes within layers are allowed, this e ect increments more by construction the variance of the layer multidegree centrality and its distance as compared to the nodewise counterpart.

. . Closed-form of the distance for random multilayer and multiplex networks in the thermodynamic limit Multilayer case Using Eq. . , the layer MC variance for random multilayer networks read:

2 ( Y ) = 2 (2 + + ) ( . ) = 4 2 ( ) + 2 ( ) + 2 ( ) = ( -1)( -1) (1 -)(1 + 2 -1 + 1 -1 ) Similarly, 2 ( X ) = 2 ( -1) (1 -)(1 + 2 -1 + 1 -1
). In the thermodynamic limit ( >> 1, >> 1), we saw that the nodewise and layerwise distances depend on the standard deviation of their respective MC distributions. Therefore, the nodewise standard deviation scales as

( X ) = √
and the layerwise one scales as ( Y ) = √ in the thermodynamic limit. Interestingly, the standard deviation of the random multilayer networks are trivially related X = Y . Therefore, the random multilayer network is an example where the first and second moment are trivially related across descriptions. However, the comparative reorganization of the initial random multilayer networks provide a non trivial relationship between distances. Using the previous equations and table we find that for random multilayer networks, the distances read: with , the probability parameter. Notably, changes that involve different nodes with layer being unchanged ( → 1) do not a ect layerwise distances, and viceversa changes involving di erent layers with nodes being unchanged ( → 1) do not influence nodewise distances (Fig. . ). Both distances increase linearly in and . The above formula implies that:

( X , X ( )) = (1 - ) (1 -) ( . ) ( Y , Y ( )) = (1 - ) (1 -), (
( ) = 1 - 1 - ( ) ∝ ( ) ( . )
The formula shows that the ordering between distances is only influenced by the proportion of changes visible in one or the other description. We validate this equation by computing distances between random multilayer networks and a numerical rewired version of them (Fig. . ). Although the formula are computed for expected MC sequence, the fit is decent even with only one instance of a rewired MC sequence. Taken together, Eq. . and Eq. . demonstrated the complementarity of nodewise and layerwise distances and established a formal node-layer duality in multilayer networks.

Multiplex case

In the case of multiplex networks, distances as defined in Eq. . apply as well. For realistic purposes, we consider the case where the multiplex network has all its replica links that are not rewireable, i.e. , the term = 0 (see previous Table ). Since we deal with and ( Y ) = 2 ( -1) (1 -). In the thermodynamical limit, the distances for random multiplex networks read:

( X , X ( )) = (1 - ) √ (1 -) ( . ) ( Y , Y ( )) = (1 - ) √ 2 (1 -) ( . ) Y = √ 2 1 - 1 - X , ( . )
We validate these equations by computing distances between a random multiplex network and a numerical rewired version of it (Fig. . ).

Although the formulae are computed for expected MC sequence, the fit is decent even with only one instance of a rewired MC sequence. The layerwise distance is √ 2 times higher than the nodewise one for the same condition as compared with random multilayer networks. Notably in multiplexes ∝ √ 2 by construction, so that their dual layerwise description a-priori better emphasizes di erences with respect to the primal nodewise one (Fig. . ). . . Scaling of distances for random multilayer and multiplex networks in the thermodynamical limit

We want to know where the closed-form of the distances for multiplex networks reach their maximum in function of the size , of the network. In general, the higher is and the higher are those distances (for fixed probability parameter). Therefore, here we fix + = , with c a constant to limit the size of the network. In the thermodynamical limit, X and Y have the same value for every combination of , in random multilayer network (see Eqs. . , . ). Morever they trivially reach their maximum when = (Fig. . ).

For random multiplex networks, X and Y have the same , dependence, i.e., they reach their maximum for the same combination of and (Fig. . ). Therefore, we use only X to derive the maximum of both distances (Eq. . ). We set = -, as the variable of interest. Therefore Eq. . read:

X = (1 - ) ( + ) √ - 2 √ 2 ( . )
X = (1 - ) -3 2 √ - ( . )
By setting X = 0, it implies that the maximum of X is * = 3 leading to the position of the maximum when :

( -) * = + 3 ( . )
By construction, layerwise distances are always higher than nodewise one for multiplex networks (Fig. . ). However, they reach their maximum at the same non-trivial point when is higher than (Eq. . ).

. Discussion

Multilayer network science has provided a common language and methodology for understanding and analyzing the complex structure and dynamics of multilevel systems from a node-centric perspective. However, in a multilayer network the layers constitute also a basic component and the properties of the system can be similarly studied from a layer-centric viewpoint [ ]. Recent studies have for example inspected the layer properties to quantify the amount of redundancy in the networks [ ], or to define notions of clustering [ , ], or ranking [ , , ] on layers (see Sec. . . ). While these studies provided intuitions on how to look at the system connectivity from a layerwise perspective, a formal characterization of its properties and potential benefits is still lacking. By focusing on the multidegree centrality (MC), i.e., the total number of connections that a unit has within and between di erent levels [ ], we provided a first basic characterization of the layerwise description of multilayer networks. We showed that the node-layer duality is of interest when characterizing the structural properties of nodes and layers "at the second order" in terms of moments. Based on a perturbation approach, we showed that local changes that do not alter the nodal connectivity, are only visible in the layerwise dimension. We derive a stochastic rewiring model which allows to disentangle link perturbations that are visible in one description but not in the other. We showed that, in general, the formal duality made by the MC-based distances in this perturbation approach is a complex interplay between the size of the network, the organizational di erence between the compared networks (quantified by and ) and the nodewise and layerwise heterogeneity (i.e., the variance) of the MC distributions. In the case of random multiplex and multilayer networks we derived closed-forms of the distances which allows to uncover a closed-form of the duality in this case. In random multiplex networks, we show that by construction the variability of the distance is higher in the layerwise description than the nodewise counterpart, making layerwise dimension a natural candidate for studying real-world systems.

One could argue that we can study directly the local connectivity of node-layers to characterize multilayer systems. Although interesting, we think this approach is limited. For nodes and layers, the MC sequence has a size while the conjunction of node and MC sequence has a size + while having the same characterizing power (in terms of number and position of links). In this perspective, this makes the use of node-layer duality more parsimonious to describe structural properties of the system.

Our study suggests to systematically consider the layerwise dimension in conjunction with the nodewise one to characterize multilayer system. Actually, characterizing only the node properties leads to overlooking valuable hidden information on the layerwise dimension in the multilayer system of interest. This hidden information can happen to be crucial in the network comparison problem (see Sec. . . ) where it discriminates better than the nodewise dimension normal and abnormal state of networks, for example (see Ch. ). Beyond quantifying the layer properties, the node-layer principle may unveil a potential interplay between some node and layer properties in real-world systems that is certainly richer than for random multilayer networks. According to the principle of emergence, interactions and correlations between nodes and layers can bring new information about a system [ ]. In particular, how the structural symmetries of layers a ect the nodes dynamics and viceversa can be a promising venue for nodelayer duality [ ]. At the modeling level, the construction principle of the stochastic rewiring model could be used to build generative or null models that take into account both dimensions. Finally, applying layer measures on network and viceversa may reveal new unexpected properties. For instance, operating a 'node reduction' in the sense of layer reducibility [ ] is an interesting perspective in this direction. .

Preliminaries and methods

To build on our previous results on node-layer duality, we exploit the complementarity of the nodewise and layerwise description to characterize real-world multiplex networks ranging from transportation and social networks to genetic and neuronal systems (see Data description in Sec. . ).

More specifically, we test the hypothesis that the node and layer multidegree centrality (MC) distribution better classify/characterize systems than the node MC distribution alone. To do so, we compute the node and layerwise distances between real-world multiplex networks and uniformly randomized version of them. Then we normalize these distances by those obtained from equivalently rewired random multiplex networks. Finally, we plot the results in the layerwise description in function of the one of its nodewise counterpart.

.

. Connectivity-based measure of dual characterization

Each real-world multiplex network is binarized and symmetrized before analysis. It restrains our analysis to the topology of multiplex networks but allows us to use the results of the stochastic rewiring model of the previous chapter. We extract the node and MC sequence of each multiplex network. Using Eqs. . , . we compute the expected node and layer MC sequence with and defined as in Sec. . . . Therefore, we obtain the expected MC sequence of each multiplex network where the links are reshu ed anywhere uniformly at random. We finally compute the euclidean distances X , Y between the MC sequences and the expected MC sequences.

To avoid size and link density e ects we normalized the actual distances by those obtained from equivalently rewired random multiplex networks [ ]. In other words, for each multiplex network, we build its equivalent random multiplex in terms of number of nodes, layers and density. More specifically, we generate numerically the MC sequences of such random multiplex networks based on the binomial law, since numerical computations are heavy for very big systems (e.g., Twitter events N , M

). Then, we compute the euclidean distances X , Y between the random MC sequences and the expected MC sequences (Eqs. . , . ). The normalization factor is finally obtained by averaging over nodewise and layerwise random MC sequences.

In the limit of large networks, the measures reduce to (Eqs. . , . ):

X X = X X = X √ (1 -) ( . ) Y Y = Y Y = Y √ 2 (1 -) , ( . )
In essence, we make a null model analysis where we compare the fluctuations to the mean of real-world MC to the one of random MC sequence. Therefore, we characterize real-world multiplex networks by how far their MC heterogeneity is from a random MC heterogeneity. It means that a value close to indicates that the real-world multiplex is not distinguishable from a random multiplex network.

. Results

. . Separation between real-world multiplex networks

We plot the logarithm of the normalized layerwise distance in function of the logarithm of the normalized nodewise distance (Fig. . ).

Results confirmed that adding the layerwise dimension allows for a better separation of multiplex networks, which would be otherwise indistinguishable by only looking at nodewise distances (e.g. , German Transport networks completely overlaps with Uganda village one and partially with Genetic networks).

Since MC-based distances are closely related to the standard deviation of the real-world MC sequences, this separation was mainly due to the higher heterogeneity of the layer MCs as compared to the node ones (Tab. . ). Furthermore, the distributions are not single-scaled in realworld networks , almost all the networks exhibiting a heterogeneity at least order of magnitude higher than random multiplex networks (x and y-axis di erent from zero on Fig. . ). .

. Real-world multiplex networks clustering in the dual space

Using a k-means clustering we found that networks tend to optimally separate into two subgroups (Fig. . ). Those with relatively low X and Y values were mostly associated with systems physically embedded in space (e.g., German transport, EUAir, C.Elegans connectome). Those with relatively high nodewise and layerwise distances were instead not characterized by a strong spatial connotation (e.g., PierreAuger, Arxiv, Twitter events). While there were few exceptions (e.g., Human microbiome), these results could be explained by the typical limited node degree heterogeneity of spatial networks due to environmental ) except for the FAO trade multiplex network. This indeterminacy concerning the latter could be due to the fact that we do not take into account weights and direction of interaction which we hypothesize can be detrimental in the characterization of this particular system.

. At a smaller scale, the clustering of Arxiv and PierreAuger is remarkable. Although they come from di erent datasets and have di erent sizes, and density (Fig. . ), they have a similar distance to their uniformly randomized expected MC sequence. This could suggest that these combinations of nodewise and layerwise heterogeneity is characteristics of multiplex networks of scientific collaboration. Indeed, this hypothesis requires a validation on more scientific collaboration datasets.

The German Transport classes represent urban transportation systems of german cities. Since it is the most populated class of our dataset, we decided to exploit the complementarity of both descriptions to investigate if their interplay is characteristics of these systems. In Fig. . , we see that the normalized layerwise distance is positively correlated with the normalized nodewise one (Pearson, = 0.68, = 9.6 -6). Furthermore, the regression line ( = 0.44, 2 = 0.45) shows that the normalized layerwise distance tends to evolve slower to normalized nodewise distance in terms of order of magnitude. In other words, the nodewise MC heterogeneity tends to evolve comparatively more than the layerwise one in these systems. This could be explained by the fact that most new bus,train and metro lines have a more or less standard number of stops and connections but that they all connect to few hubs, the latter driving the nodewise heterogeneity.

[ ]: Arenas et al. ( ), 'Optimal map of the structure of complex networks',

. Visualisation of multiplex networks

To help understanding how the nodes contribute to layers and viceversa, we adapted the method proposed by [ ] to visualize partitioned networks. To do so, we first computed the contribution matrix containing the number of links that a node share with a layer . We next considered a typical truncated singular value decomposition (tSVD) = Σ † , where Σ is a diagonal matrix containing the singular values, and and are respectively the left and right orthogonal matrices associated with the nodewise and layerwise dimension, respectively. To visualize the nodewise contribution we projected the space spanned by the first two left singular vectors, i.e., (1)-( 2). In this D space, the layers are represented as lines whose direction depends on their cohesiveness, i.e., layers that share many links tend to be represented along similar directions. The nodes are instead represented as points. The more the nodes contribute to a specific layer the more they tend to be aligned to its direction. The distance of each point from the origin is finally proportional to its MC. Similar visualization can be obtained for the layerwise contributions by projecting the space spanned by the first two right singular vectors, i.e., (1)-( 2).

In Fig. . , we show examples of projection plots using Singular Value Decomposition. First, we see that we can appreciate the heterogeneity of connectivity of nodes and layers by looking at how the markers are spread in the D space. In particular we recover that for all the multiplex networks in Fig. . , the layerwise heterogeneity is comparatively bigger in the layerwise than in the nodewise description. Second, we see that the C.Elegans multiplex connectome (Fig. . a) and Arxiv multiplex (Fig. . b) exihibit nodes and layers markers that lie in the middle of the layers and nodes direction respectively. This means that nodes are globally well integrated into each layer and that each layer is globally well integrated into nodes. On the contrary we see that Chemnitz (Fig. . c) and PierreAuger ( Fig. . d) multiplex networks exhibit nodes and layer markers that lie in orthogonal directions. This means that nodes are connected into some layers but with few or no connections in other layers and similarly for layers across nodes. The inset in Fig. . d) is a good summary of the illustrative power of the method. It shows that many nodes are well integrated into all the layers but one (the more horizontal one). This lasting layer has many nodes that have links that lie almost only it. Then in the middle we observe nodes that share about the same proportion of links with the two groups of layers, i.e., the more horizontal one and the others. .

Discussion

Network comparison and classification is an essential task of network science. (see Sec. . . ). In essence, it consists in deriving su cient statistics to quantify how similar or dissimilar the networks are. For instance, it can be used to check the relevance of a generative network model: the closer the model to the target, the better. On the other hand, if the similarity of networks is known a priori, it allows to evaluate the characterizing power of measures used for comparison. For instance, when comparing between normal and abnormal state of the brain we can assess whether or not the measure is predictive of the abnormality (see Sec. . ).

Based on the local connectivity of nodes and layers, the node-layer duality enabled us to better discriminate di erent multiplex real-world systems and capture intrinsic properties not included in the model such as their spatial nature. Furthermore, the relative clustering in space of the di erent classes of networks (e.g., Arvix and PierreAuger) suggests that real-world multiplex networks could be characterized by their MC heterogeneity pairs. At the level of a class of networks, we show that German Transport networks seem to be characterized by a particular interplay between their node and MC heterogeneity. Yet, we show in the previous chapter that node and layer MCs are theoretically uncorrelated. Overall, it suggests that the node-layer duality framework is well suited to characterize and potentially o er new insights for multiplex systems resulting from the non-trivial interplay between nodes and layers properties.

Binarization of weighted multiplex networks and symmetrization of the directed networks constitute the main limitation of this work. Indeed, these operations can have distorted the classification. For instance, by binarizing and symmetrizing FAO Trade, we put on the same scale millions of dollars exchange with thousand dollars ones and we disregard commercial unbalance between countries. This lay the space for future refinement in particular in adapting the stochastic rewiring model to weighted and directed multilayer networks.

Inspired by [ ] we show that the SVD-decomposition based projection greatly illustrates qualitative connectivity-related features of multiplex network. By its ability to represent the whole network in a synthetic way, we consider it can be a good complement for multiplex network visualisation compared to softwares that actually draw multiplex networks. [ ].

While these results were obtained using relatively simple network quantities and models, other features can be explored within the nodelayer duality framework in an e ort to improve the overall classification and result interpretation. For instance, we can classify multiplex networks using clustering coe cients in both descriptions. Another idea could be to extend the concept of small-world index [ , ] to the nodewise and layerwise classification. These investigations would also require the use of other null models such as the one that preserves the multi-degree sequence (see Sec. . . ).

. Data

The first eight datasets are available at Manlio De Domenico's repository. Uganda village dataset is available in Netzschleuder repository. German Transport is available at Kai Bergermann's GitHub repository. All real-world multiplex networks were symmetrized and binarized before any analysis. Their size and density are summarized in Fig. . : . Twitter events: layers, from to nodes. Nodes represent Twitter users, layers correspond to retweet, mentions and replies between users. Data are acquired througt the data API and in a time interval that correponds to a specific event and filtered by keywords related to it [ , ]. In this chapter, we build on these observations to test our dual approach in multilayer networks derived from experimental neuroimaging recordings. We evaluate the impact of Alzheimer's disease (AD) on local connectivity changes across di erent brain areas and across di erent frequencies of brain activity. To do so, we considered multilayer brain networks estimated from source-reconstructed magnetoencephalography (MEG) signals in a group of AD human subjects matched with a group of healthy controls HC (Methods in Sec. . ). Specifically, we used bispectral coherence to simultaneously infer weighted interactions among regions of interest (the nodes) within and between di erent signal frequencies (the layers) [ ].

. Results

. . Discriminative power of primal and dual descriptions

First, we wanted to investigate the nodewise and layerwise distances between the AD patients and the average HC. Since the related multilayer networks are complete and weighted, we compute the multistrength centrality for each AD patients and the average HC, i.e. , the sum of the weights adjacent to the interacting unit of interest. In this weighted case, distances captures simultaneously the di erence in weighted topology but also the mere di erence in global connectivity. Therefore, we check that the total weights of patients are roughly equal to the one of the average HC (Fig. . ).Results show that there is no significant di erence between nodewise and layerwise distance in the initial configuration (Fig. . ). However, as we decrease the number of layers the more significant was the di erence between Y and X in favour of Y . This behavior did not solely result from the reduction in network size, but instead originated from the greater consistency in connectivity changes within frequencies compared to between frequencies (Fig. . ). This suggests that, globally, the layerwise description better discriminates AD patients from HC than the nodewise one as the frequency resolution gets coarser. negative (Δ X = -0.44, 2 = 0.23) which means that highly connected regions of interests (ROIs) tend to be the most disrupted in AD patient whereas weakly connected ones are less disrupted or even enhanced. Similarly, the layer disruption index Δ Y is negative (Δ X = -0.38, 2 = 0.78) meaning that highly connected layers are the most disrupted in AD patients. Also, we note that the linear decay has quite good performances in describing the disruption of the layer multistrength centrality.

In order to relate this global reorganization of connectivity to behavior, we conduct a correlation analysis between the Mini-Mental State Examination (MMSE) score [ ] and the individual node and layer disruptions indices for each patient (Fig. . c-d). Individual disruption indices are defined as the slope of the line obtained by regressing the di erence between the AD patients and the average HC multistrength centrality. Results show that the individual node disruption indices do not significantly correlate with the MMSE score of AD patients (Spearman = -0.12, = 0.58). On the contrary, the individual layer disruption indices significantly and positively correlate with the MMSE score of the AD patient (Spearman, = 0.60, = 2.2 -3). In other words, the higher the layer disruption index, the higher is the MMSE score and so the less severe is the cognitive decline of patients. Of note, the closer is the disruption index from zero or slightly positive, the more patients resemble the average HC. Those results hold as we perform the analysis without the extreme outlier having Δ > 2, Δ < -4 visible in (Fig. . c-d) in the nodewise (Spearman, = -0.12, = 0.60) and the layerwise descriptions (Spearman, = 0.62, = 1.9 -3). Therefore, the layerwise description seems to provide a predictor of the severity of the disease where the nodewise one does not.

. . Dual characterization of the local disruption of connectivity

Next, we investigate the local reorganization of the node and layer multistrength centrality (Fig. . a-b). For each ROI or frequency, we statistically compare the multistrength centrality distribution of AD patients and HCs using a non-parametric Wilcoxon test. Compare to previous results we perform an all-to-all comparison between the AD patients and the HCs and not only AD patients against the average HC. Results show that each ROI and each frequency has a disrupted functional connectivity in AD patients as compared to HCs ( < 0-Fig. . a-b). Specifically, the paracentral lobules and left precuneus were among the most impacted areas. In the layerwise side, several brain frequencies within the alpha range (8 -13 Hz) in the AD group presented reduced multistrength centrality as compared to HC (Fig. . a-b). However, these local network changes were not statistically significant (| | < 1.96).

In order to relate the local changes of functional connectivity to the behavior of AD patients, we performed a correlation analysis between the MMSE score of patients and each ROI and each frequency's multistrength. centrality The most predictive areas were the caudal anterior cingulate cortex in both hemispheres, [ , ] (Fig. . c) but also the right precuneus, right inferiorparietal and right superiorparietal (Tab. . ). The latter areas survive the mean FDR test but not the test for multiple comparison. Notably, a higher number of strongest associations were found in the frequency-wise description. The most significant ones were all located within the alpha frequency range (Tab.

. ), e.g. , the Hz MEG frequency (Fig. . d).Most of the frequencies of the alpha band survive the mean FDR but they barely survive the correction for multiple comparison ( < 0.055, FDR corrected -Tab. . ).

Altogether, these results suggest that the layerwise dimension is more predictive that the nodewise one in terms of the severity of the disease, when considering the connectivity of the frequency with crossfrequency coupling computed with bicoherence.

. Discussion

The brain is a complex network exhibiting multiple interactions across di erent scales and levels [ ]. The interaction between neural oscillations at di erent frequencies is crucial for coordinating and integrating information across di erent brain regions and cognitive processes [ ]. For example, phase-amplitude coupling between theta (4 -8 Hz) and gamma (30 -100 Hz) oscillations is involved in working memory and attentional processes as well as sensory perception [ ]. Disruptions in cross-frequency coupling have been observed in several neurological and psychiatric disorders, including epilepsy, schizophrenia, and Parkinson's disease [ ]. Understanding these abnormalities may provide insights into the underlying mechanisms of these disorders and help identify biomarkers and novel treatment approaches.

Alzheimer's disease (AD) is the most common form of dementia accounting for a significant proportion of cases worldwide and placing a substantial burden on healthcare systems and economies [ ]. AD is characterized by a progressive neurodegeneration and anatomical atrophy followed by the disruption of multiple functional large-scale network properties [ , , , , , ]. These include the loss of key regional hubs in the parietal and temporal regions within alpha frequency ranges [ , , , ] but also between theta/alpha and gamma bands [ ]. Notably, disrupted theta-gamma coupling has been linked to deficits in working memory, which is commonly a ected in AD, and to a lesser extend to worse clinical outcomes [ ]. While these studies have provided fresh evidences on how AD brain networks change within and between frequencies, the existence of possibly higher-order e ects between multiple frequencies has been systemically overlooked. We first showed that the frequency domain has a better discriminative power than the region domain as the frequency resolution becomes coarser. This is due to more consistent changes in connectivity within frequencies compared to between frequencies. This could suggest that when considering specific frequencies associated with specific cross-frequency coupling (e.g. ,theta-gamma coupling [ ]) the layerwise description should be favored in terms of discriminative power between AD patients and HC.

We recovered the global disruption of functional connectivity in the frequency-space and the region-space. Results on the regionwise disruption of local connectivity are in line with literature where crossfrequency coupling are not integrated [ , ]. We found that individual layer disruption indices were associated with the cognitive decline of AD patients as measured with the mini-mental state examination (MMSE) [ ]. In particular, AD patients with disruption index close to or slightly above zero, i.e. , patients similar to controls have the best tested performance. We did not find such an association for individual disruption indices of regions. This suggests that the global reorganization of frequency coordination may be predictive of the clinical outcome of AD patients and should therefore be preferred to the regionwise one.

We finally found that both nodewise and layerwise multistrength cen- trality decrements were associated with more severe cognitive decline in AD patients as measured by the MMSE score (Tabs. . , . ). In the brain space, the most predictive area was the caudal anterior cingulate cortex in both hemispheres, a well-known hub of information processing in the brain that plays an essential role in AD pathophysiology [ , ]. Notably, a higher number of strongest associations were found in the frequency space. The most significant were all located within the alpha frequency range, considered one of the most reliable non-invasive functional predictors of AD-related cognitive symptoms [ , ] .

Node-layer duality o ers a natural framework for capturing these properties from multilayer brain networks. Here, we showed that AD is characterized not only by the connectivity disruption in specific brain areas, but more importantly by the aberrant coordination between frequencies. Specifically, the cognitive decline of AD patients was significantly associated with the reduced capacity of the alpha band to facilitate cross-frequency information integration.

Indeed, this work could be improved either by increasing the population size or exploring di erent network sizes (i.e., di erent parcellation schemes) or by using di erent network measures. Actually, given the unavoidable presence of noise and artefacts, it is often di cult to detect groupwise di erences based on local connectivity measures [ ]. Also, as discussed in Ch. , the brain network reconstruction supposes arbitrary choices such as the selection of the connectivity measure.

With bicoherence, we focus on non-linear phase coupling within and between frequencies. Although results can be consistent between di erent protocols and connectivity measures (e.g. , the disruption of the alpha band in AD), one should be aware that our results can be attached to a particular aspect of the brain activity. Staying in the same protocol, the results should be validated by keeping the asymmetry of the bicoherence links. Furthermore, the correlation analysis including the MMSE score should be corrected for age and level of education since it has been shown that they influence the outcome of the score [ ].

A thrilling perspective, would be to derive a mechanistic model of the functional neurodegeneration of Alzheimer's disease based on the node-layer duality. For instance, using the parameters of the stochastic rewiring model, we could test whether the disease is a multilevel randomization of normal brain functions at the node and layer level with a potential predictive power on the severity of the disease.

In conclusion, our approach provides the first proof-of-concept accounting for complementary aspects of multilevel brain systems with broader implications in basic and clinical neuroscience.

.

Material and Methods

For reproducibility, we provide an Owncloud repository (hyperlink included in the online version) containing the supra-adjacency matrices described in the following sections.

.

. Participants

The research was conducted in accordance with the Declaration of Helsinki. A written informed consent was obtained from subjects after explanation of the study, which was approved by the local ethics committee of the Pitie-Salpetriere Hospital. .

. Experiment

Magnetic resonance imaging (MRI) acquisitions were obtained using a T system (Siemens Trio, -channel system, with a -channel head coil). The MRI examination included a D T -weighted volumetric magnetization-prepared rapid gradient echo (MPRAGE) sequence with the following parameters: thickness mm isotropic, repetition time (TR) , ms, echo time (TE) . ms, inversion time (TI) ms, acquisition matrix

x .

The magnetoencephalography (MEG) experiments consisted in a resting state with eyes closed. Subjects were placed in magnetically shielding room and were asked to remain at rest. MEG signals were collected using a system that comprised magnetometers and planar gradiometers (Elekta Neuromag TRIUX MEG system) at a sampling rate of , Hz and online low-pass filtered at Hz. The ground electrode was located on the right shoulder blade. To favor physiological artifacts reduction, an electrocardiogram (EKG, Ag/AgCl electrodes) was placed on the left abdomen and a vertical electrooculogram (EOG) was simultaneously recorded. Four coils were fixed to the participant in order to monitor head position and to enable coregistration with the anatomical MRI. To digitize the physical landmarks (namely, the nasion, the left and right preauricular points) we used a Polhemus Fastrak digitizer (Polhemus, Colchester, VT). Three consecutive clean epochs of approximately min each were extracted for further analysis.

. . Processing

To remove external noise, a signal space separation was performed using MaxFilter. An in-house software was used to remove cardiac and ocular blink artifacts from MEG signals via a principal component analysis. MEG epochs that still presented contamination were removed after visual examination. At the end of the process, we obtained three Hz included for a total of M layers. The connectivity measure is the bicoherence [ ] which is defined to consistently estimate the connectivity within and across frequency layer and was proved to grasp non-linear similarities between ROIs. By definition, values of bicoherence lie between and . We estimate the connectivity between any pairs of nodes at any frequency with non-overlapping windows of length s. Since the bicoherence is asymmetric, we symmetrized the matrix by selecting the highest value of bicohence between any two symmetric pairs in order to obtain an undirected network. The average weight di erence between the symmetric and asymmetric version is .

in a favor of the symmetric version (patients and controls included) with a maximum di erence of . .

Conclusion

In this manuscript, we elaborated on the notion of topological duality in multilayer networks. In a broader sense, duality refers to a complementary, often antithetic representation of a given phenomenon or object. For example, the electric and magnetic fields are dual objects, always orthogonal and induced by the temporal variations of one with respect to the other. Actually, Maxwell's tour de force was to recognize that both fields are the complementary descriptions of a single physical object, the electromagnetic field.

Beyond its conceptual appeal, duality is also often a crucial operational concept. By switching descriptions, it often allows di cult problems to be mappped onto easier ones. One of the most famous examples of this is the Fourier transform which maps problems from the timedomain to the frequency domain. By doing so, it transforms any timedependent di erential equations that are potentially di cult to solve into a simple algebraic equation.

Another important aspect of some dualities is that they can map systems in which variables are quantitatively large to systems in which quantities are quantitatively small. This type of duality is exemplified by the Kramers-Wannier duality in statistical physics. In fact, this duality establishes a formal relationship between strong and weak spin coupling in an Ising model [ , ].

In network science, the existence of dual communities in planar graphs [ ], the duality between equilibrium and growing networks [ ], the dual relationship between structure and dynamics [ , ], are just a few examples witnessing the fertile research area in the field. Network duality is important not only from a fundamental perspective, but also because it allows to reuse tools belonging to one discipline in the dual one, and viceversa. For example, the geometric mapping between time series and networks unlocks the use of signal processing techniques in complex networks, as well as the use of graph theoretic methods for time series analysis and modeling [ , ].

While the dual properties of complex networks have recently begun to be explored, their definition and meaning in multilayer networks

has not yet been introduced. In this manuscript we proposed the node-layer duality to provide a first intuitive way to simultaneously characterize structural properties of multilayer networks in both their primal nodewise description and their dual layerwise counterpart.

We first defined the node-layer duality as two equivalent representations of a multilayer network where nodes are the interacting units on the one hand and layers are the interacting units on the other hand. We showed that the two representations are not trivially related with respect to the second moment of the distributions of the structural properties. Then, we demonstrated the complementarity of the two descriptions in characterizing one aspect of multilayer systems.

In particular, we showed that certain changes in local connectivity were visible in one representation and not in the other. Also, we showed that the wiring map and the number of connections of a multilayer network are constrained by both the node and layer multidegree centrality (MC) sequences. We derived a stochastic rewiring model that allowed us to systematically investigate the complementarity of the representations by perturbing the wiring of the connections. We analytically derived the expression of the perturbed MC sequences in both representations which we subsequently summarized in distances between the perturbed sequences and the initial ones. In the limit of large networks, we established a formal node-layer duality relationship using MC-based distances.We derived a closed-form of this duality in the case of random multilayer and multiplex networks. By focusing on the MC, i.e. , the total number of links that a unit has within and between di erent levels, we provided a first basic characterization of the node-layer duality of multilayer networks.

Then, we exploited this basic characterization of the node-layer duality in an attempt to characterize di erent classes of real-world multiplex networks. We showed that the complementarity of the nodewise and layerwise descriptions allows us to better discriminate between the real-world multiplex systems at hand. We showed that they form two clusters characterized by their spatial connotation which is not a priori included in our model of local connectivity. The results also suggest that di erent classes of systems could be characterized by a particular interplay between the two descriptions. We have exploited the complementarity of our characterization of the node-layer duality to provide a method for visualizing multilayer networks. In contrast to classical representations of multilayer networks, it allows to summarize basic properties related to the connectivity and to represent the whole multilayer network rather than just a part of it.

Next, we used our characterization of the node-layer duality to study the functional neurodegeneration in Alzheimer's disease in a multilayer network model that integrates cross-frequency coupling. We hypothesized that one of the descriptions may be more relevant to this problem than the other, a typical feature of dual representations. We showed that the layerwise (i.e. , frequencies) description has a better discriminative power between AD patients and HC compared to the nodewise one (the nodes). Then, using a global measure of the connectivity changes between the abnormal and healthy state, we found that the global disruption of frequency connectivity could be a good predictor of the cognitive decline of in patients. We did not find any predictive power on the nodewise description. Regarding the connectivity of each interacting unit, many regions and frequencies were associated with the cognitive decline in patients. However, the strongest associations were found in the layerwise description, especially in the alpha band which is one of the most consistent biomarkers of the functional AD neurodegeneration in the literature.

In conclusion, we have shown that the node-layer duality has the distinctive features of dualities (as discussed in the previous paragraphs) that can be further exploited. Through the stochastic rewiring model and the study of real-world multiplex networks we have shown the complementarity of the two descriptions. This could inspire generative models that act on both descriptions simultaneously. For instance we can model the functional neurodegeneration of Alzheimer's disease according to these two perspectives. In our work, we focused on the local connectivity but investigating other measures such as the clustering coe cient or the small world index could be a natural continuation to exploit the complementarity of the node-layer duality approach. .Emergent phenomena in complex systems are said to be rooted in the interaction of their components. Similarly, the correlation of nodes across layers are believed to be at the core of emergent phenomena in multilayer networks as compared to their single-layer counterpart [ ]. In this perspective, the node-layer duality paves the way to study the interaction between nodes and layers in multilayer networks. In fact, the non-trivial interplay between nodes and layers structural and dynamical properties can potentially uncover invariants that characterize multilayer network across and within a particular class (biological, infrastructure, social...).

The study on Alzheimer's disease showed an example of the strong relevance of the dual description compared to the weak relevance of the primal one. This suggests that both descriptions should be systematically explored to uncover which one is the most relevant for a given problem. Furthermore, it may prompt new hypotheses about why and how one or the other description is more relevant.

Since the role of nodes and layers are intertwined from one description to the other, we can reuse tools that were defined for the nodewise description on the layerwise description and viceversa. This can enrich our description of both nodes and layers.

Considering the broad scope of multilayer network models, we hope that the echo of node-layer duality will be felt universally. and conditions: return [ , , , ] return None def select_ inter(A,N,M,uplet,rng, com,trial): """Return the hole where 'uplet' is rewired according to the p_ layer event, which keeps the node indices fixed Input: A -supra-adjacency matrix N -number nodes M -number of layers uplet -tensor indices (layer,layer,node,node) of the link to rewire rng -random number generator com -dictionary with two keys ["be interlayer"] and ["replica"] trial -number of proposed new rewired position Output: [ , , , ] -Position where the link 'uplet' is rewired, it links node i in layer a to node j in layer b """ k, l, i, j uplet[ ], uplet[ ], uplet[ ], uplet[ ] the coordinates of the selected uplet rand_ M list(rng.integers(M, size [

, 2])) for a, b in rand_ M: conditions build_ instructions([ , , , ], com) avoid_ self_ loops not (a b and i j) condition to avoid rewiring in self-loops row, col to_ supra_ index(N, [ , , , ]) if A [row,col] and conditions and avoid_ self_ loops: return [ , , , ] return None def select_ tel(A,N,M,uplet,rng,com,trial): """Return the hole where 'uplet' is rewired according to the p_ tel event, which changes at least one of the layer index and one one the node index Input: A -supra-adjacency matrix N -number of nodes M -number of layers uplet -tensor indices (layer,layer,node,node) of the link to rewire rng -random number generator com -dictionary with two keys ["be interlayer"] and [row,col] and conditions: return [ , , , ] return None def stochastic_ rewiring(s_ A,multi_ degree,multi_ degree_ DS,event_proba,R,dic_ instructions,number_ of_ trial_ per_ event ): """Return the hole where 'uplet' is rewired according to the p_ tel event, which changes at least one of the layer index and one one the node index Input: s_ A -supra-adjacency matrix multi_ degree -node multidegree sequence multi_ degree_ DS -layer multidegree sequence event_ proba -array [ _ , _ , _ ] giving the probability of each related rewiring event R -proportion of links to rewire (from to ) dic_ instructions -dictionary with two keys ["be interlayer"] and ["replica"] number_ of_ trial_ per_ event -number of proposed new rewired position Output: s_ B -supra-adjacency matrix of the rewired multilayer network l_ K -node multidegree sequence of the rewired network l_ Ktranspose -layer multidegree sequence of the rewired network number_ of_ unrewire_ edge -number of links that are not rewired due to any problem """ R int(round((s_ A.nnz * R))) convert the proportion into number of links number_ of_ unrewire_ edge s_ A s_ A.tocsr() s_ A s_ A s_ A.T s_ A s_ A.tolil() s_ B s_ A.copy() save the original matrix, s_ B is the one that is rewired l_ K multi_ degree.copy() l_ Ktranspose multi_ degree_ DS.copy() N len(l_ K) number of nodes M len(l_ Ktranspose) number of layers for r in range(R): iteratively rewire R links row,col s_ B.nonzero() row and col of all existing link in the network random_ index rng.integers ( , len(row), size ) select one link uniformly at random among all exisitng link while s_ B[ while new_ uplet None: event choices(["intra", "inter", "tel"], event_ proba, k )[ ] randomly choose the rewiring event according to the probability in event_ proba if event "intra": p_ node rewiring event new_ uplet select_ intra(s_ B,N, uplet, rng, dic_ instructions, num-ber_ of_ trial_ per_ event) elif event "inter": p_ inter rewiring event new_ uplet select_ inter(s_ B,N,M, uplet, rng, dic_ instructions, num-ber_ of_ trial_ per_ event) elif event "tel": p_ tel rewiring event new_ uplet select_ tel(s_ B,N,M, uplet, rng, dic_ instructions, num-ber_ of_ trial_ per_ event) security if security : if the structure of the net makes it impossible to move the uplet so we keep it at its place new_ uplet uplet hoping it is not the case for a lot of nodes (normally not because sparse net) print("activate security") number_ of_ unrewire_ edge k, l, i, j

.

  from [ ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Illustration of a Barabási-Albert network (N ) and its node degree distribution. . . . . . The possible connected -nodes subnetworks with directed links. Adpated from [ ]. . . Communities and rich clubs. (A) A network composed of four communities that are linked by hub nodes (black). (B) With the addition of further inter-module connections hub nodes now form a densely interconnected rich club, consisting of nodes with a degree of or higher.Figure and caption adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Illustration of Maslov-Sneppen rewiring procedure. The links of two pairs of nodes are reshu ed in a way to keep the node degree sequence fixed. . . . . . . . . . . . . . . . . . . . . Main configurations of multilayer networks. Panel a) Full multilayer network. Both within-and between-layer connections are allowed with no specific restrictions. This configuration could typically be adopted to model multifrequency brain networks. Panel b) Multiplex network. Only interlayer connections between the replica nodes are allowed. No restrictions on connections within layers. This configuration is typically used to model real-world systems with multiple interactions such as multimodal brain networks. Panel c) Temporal network. Interlayer connections are allowed only between adjacent layers. No restrictions on connections within layers. This configuration is typically adopted to model time-varying networks. Image credit: Thibault Rolland. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Example of all possible multilinks in a multiplex network with M layers and N nodes. Figure and caption adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Schematic representation of the construction of a multiplex network with a given multidegree sequence. Stubs of di erent types are associated with nodes. Solid line, dashed line and dot-dashed line stubs refer to multilinks ( , ), multilinks ( , ) and multilinks ( , ) respectively. Stubs of the same type are randomly matched, forming the multiplex network. Figure and caption adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Layer clustering based on the entropy of network ensemble. a) The entropy , in layer with respect to a property is computed and then compared with a randomized entropy distribution with mean [ , ( ) ] and standard deviation [ , ( ) ]. Θ , ( ) measures how far is the computed entropy compared to the mean of the null distribution (it is a z-score).

  Adapted from [ ]. . . . . . -Tesla MRI of a normal human brain. Left) T image. Right) T image. Adapted from [ ]. . . . Schematic of the di usion of a water molecule within the fiber bundles. Adapted from [ ] . Illustation of the alternance of gyri and sulci at the surface of the cerebral cortex. . . . . . . EEG (up) and MEG (down) apparatus. Adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . Deterministic tractography. Seeds points are represented by white circles. They are the starting point to reconstruct, voxel-by-voxel, the principal direction of water di usion (white lines). Adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Pipeline for functional brain networks modelling and analysis. Nodes correspond to specific

  Figure and caption adapted from [ ]. . . . . . . . The di erent stages of Alzheimer's disease neurodegeneration. . . . . . . . . . . . . . . . . . . Bottom-up hierarchical modeling. Panel a) The so-called K-set hierarchy showing the model

..

  Pictures and captions adapted from [ ] (panel a) and [ ] (panel b). Inset b adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . Structural reducibility of multifrequency brain networks. Panel a) For each combination of layers a quality function measures the amount of new information added with respect to an equivalent single-layer model. Panel b) Median values of quality function obtained from fMRI multifrequency brain networks in healthy subjects. Shaded areas indicate the standard deviation around each value. Pictures and captions adapted from [ ]. . . . . . . . . . . . . . . Emergent properties in full multilayer brain networks Panel a) Intralayer and interlayer links in the multifrequency MEG network. -link between regions at the same frequency;link of the same area between di erent frequency bands; -link between di erent nodes at di erent frequency bands. Panel b) Algebraic connectivity 2 as a function of the total interlayer connectivity ( ). The vertical solid line corresponds to the actual value of interlayer connectivity, i.e., without modifying their weights. Pictures and captions adapted from [ ]. . . . Multiplex motif analysis of multimodal brain networks. Panel a) Structural-functional -layer brain network. Interlayer links between replica nodes are omitted for the sake of visibility. Five nontrivial multiplex motifs of two nodes are possible based on the type of connectivity in the DTI structural layer and in the fMRI functional layer. The Z-scores show the motifs that are overrepresented and underrepresented as compared to equivalent random networks. Panel b) Patterns of multiplex triangles comprising directed structural tuples (solid connections) closed by a functional link (dashed connections). The overall motif counts normalized by equivalent random multiplexes are illustrated as a function of basal activation parameters P and Q of the Wilson-Cowan model. Pictures and captions are adapted from [ ], (panel a), and [ ] (panel b). . Multiplex core-periphery structure of the human connectome. Scatter plot of multiplex coreness against single-layer corenesses obtained from structural (DTI) and functional (fMRI) layers. Labels indicate brain areas whose multiplex coreness cannot be predicted by looking at the coreness values in the respective structural and functional layer. Picture and caption adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temporal network flexibility correlates with brain perfomance. Panel a) An overview of network switching (or flexibility) in a temporal network. Red and blue colors identify the nodes belonging to two di erent communities according to the multilayer network modularity metric. Panel b) Brain maps of switching rate and dynamic fMRI connectivity. Values were normalized into z-scores to ensure both connectivity dynamics and switching values were scaled equally. Pictures and captions adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Stabilization of critical dynamics in multilayer glia-neuronal networks. Panel a) Left side: Glia cells redistribute metabolic resources from the bloodstream to neural synapses. Right side: Associated two-layer network model. Black arrows indicate neural synaptic interactions. Arrow thickness indicates synaptic strength which evolves according to spike time-dependent plasticity (STDP). Red arrows which terminate on black arrows represent the resource supply to the corresponding synapse. Panel b) Stability analysis of the two-layer STDP model. The largest eigenvalue of the neuronal network layer and the total resource R of all glia and synapses are illustrated as a function of time. The data plotted in black correspond to a 'baseline' condition. For the data plotted in red (labelled 'instability'), the initial evolution is the same as for the baseline data up until the di usion of resources between the glial cells is turned o (vertical arrow). Pictures and captions adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . Temporal network flexibility predicts future learning rate. Significant predictive Spearman correlations between flexibility in session and learning in session (black curve, ≈ 0.001) and between flexibility in session and learning in session ( ≈ 0.009). Each point corresponds to a subject. Note that relationships between learning and fMRI network flexibility in the same experimental sessions ( and ) were not significant; > 0.13 was obtained using permutation tests. Picture and caption adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temporal network flexibility as a clinical marker of shizophrenia genetic risk. Panel a) Significant increases in the mean dynamic reconfiguration of modular fMRI brain networks in una ected first-grade relatives (gray bar, REL) and patients with schizophrenia (black bar, SZ) in comparison to matched healthy controls (white bar, HC) [F( , ) . , P . ]. Bars indicate mean values, and whiskers represent standard error means (SEMs). Panel b) Significant increases in the mean dynamic reconfiguration of modular brain networks in healthy controls after application of dextrometorphan (DXM) [dark gray bars; repeated measures ANOVA placebo (PLA) versus DXM: F( , ) . , P . ] relative to PLA (light gray bars). Pictures and captions adapted from [ ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Multifrequency and temporal reorganization of brain networks in Alzheimer's disease. Panel a) Multiplex brain networks are constructed by layering di erent frequency-specific networks, while temporal networks were constructed by concatenating time-specific networks within frequency bands. Panel b) Top: hub disruption of MEG multifrequency networks in patients with Alzheimer's disease. Each point corresponds to a di erent brain area; k slope of the regressing line. Bottom: Brain regions with significant between-group di erence in overlapping weigthed degree. PCUN.R right precuneus; HIP.L left hippocampus; IPL.R right inferior parietal, but supramarginal and angular gyri; SPG.R right superior parietal gyrus; MOG.L left middle occipital gyrus; SOG.L left superior occipital gyrus; IOG.L left inferior occipital gyrus. Panel c) Scatter plot showing the Mahalanobis distance of each subject from AD or control group when combining of multiplex clustering coe cient (MCC) and participation coe cient (MPC) extracted from time-varying networks (gray line indicates equal distance). Pictures and captions adapted from [ ] (panel a, b) and from [ ] (panel b). . . . . . . . . . . . . . . . . . . . Multimodal brain networks reveal disrupted core-periphery structure in Alzheimer's disease. Panel a) Multimodal brain networks (multiplex) are constructed by layering DTI, fMRI, and several frequency-based MEG brain connectivity. Panel b) Spearman correlation (R . , p .

  Pictures and captions adapted from [ ] (Panel a) and [ ] (Panel b,c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . E ect of the density on the matching between analytical and numerical distances. Curves are the analytical distances as obtained with Eqs. . , . for the layer multidegre centrality (MC) and Eqs. . , . for the node multidegree centrality.Dots are the numerical distances obtained with the stochastic rewiring algorithm (see Appendix ). The numerical expected MC is obtained by rewiring of the initial networks. Error bars represents times the standard error of the mean. a) Nodewise distance in function of r for random multilayer networks ( = 20, = 20) at di erent densities ( = 1). b) Layerwise distance in function of r for random multilayer networks ( = 20, = 20) at di erent densities ( = 1). c) Nodewise distance in function of r for random multilayer networks ( = 20, = 20) at di erent densities ( = 1). d) Layerwise distance in function of r for random multilayer networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dependence of multidegree centrality distances on multilayer and multiplex networks' size.

.

  . . ). Di erent combination of number of nodes and layers are considered, with the condition + = 200. The figure boundaries are not accurate since the scaling is obtained in the limit of large networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dual characterization of real-world multiplex networks.Scatter plot of the nodewise (x-axis) and layerwise (y-axis) distances of real networks' MCs from uniformly rewired counterparts. To avoid network-size and density biases, all values are further divided by the distances obtained from equivalent random networks. The logarithm of the aforementioned quantities are shown in the plot. The optimal partition through k-mean clustering separates the network into two clusters (green and purple). Zones indicate in which cluster the networks belong to. . . . . . Optimal partition of the dual characterization of real-world multiplex networks (k-means clustering). a) Elbow method. Sum of the squared error in function of the number of clusters. The elbow point corresponds to the partition of space into two clusters which is considered as the optimum trade-o between error and the number of clusters (with knee-point detected using python package kneed). b) The silhouette coe cient in function of the number of clusters. The bigger is the silhouette coe cient, the more cohesive and separated the clusters. Its values ranges from -to . Here the maximum is obtained for two clusters, = 0.70. . . . . . . Clusters of real-world multiplex networks using hierarchical clustering. Colors of each network indicate in which cluster they belong to (green or purple). . . . . . . . . . . . . . . . . . . . . . . Projection plots of real-world multiplex networks. . In the nodewise description, layers are represented by lines and nodes are represented by markers (blue). In the layerwise description, layers are represented by markers (orange) and nodes by lines. Nodewise and layerwise plots are put on the same scale. a) Projection plots of the C.Elegans multiplex connectome ( nodes, layers). b) Projection plots of the Arxiv multiplex network ( nodes, layers). c)Projection plots of the Chemnitz multiplex network in the German Transport dataset ( nodes, layers). d) Projection plots of the PierreAuger multiplex network ( nodes, layers). Inset in the nodewise projection plot represents a zoom where markers are concentrated. . . Number of nodes in function of the number of layers for real-world multiplex networks. Size of the markers is proportional to the logarithm of the density of each multiplex network. Purple indicates networks with a strong spatial connotation and green, networks with no or weak spatial correlation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Distribution of the total multistrength centrality of AD patients compared to the groupaveraged HC. The total multistrength centralities of AD patients (orange bars) are compared to the one of the group-averaged HC(blue line). For every tested layer configuration, the ratio between the total multistrength centrality of the group-averaged AD patient over the one of the group-averaged HC is constant equal to 1.06. All AD patients except outliers have a total multistrength centrality that is very close to the one of the group-averaged HC. a) layers configuration. b) layers configuration. c) layers configuration. d) layers configuration. The layers (frequencies) are removed evenly on the range to Hz. . . . . . . . . . . . . . . . . . . Group-averaged nodewise and layerwise distances between AD and HC for di erent frequency resolutions (from M to M ). Layers are iteratively removed with a linear spacing to ensure no frequency layers are favored. The asterisk marks the number of layers ( <= 31) for which the di erence between the distances becomes significantly di erent ( < 0.05, FDR corrected). Vertical bars denote standard deviations. . . . . . . . . . . . . . . . . . . . . . . . . . . Distribution of connectivity changes in AD multifrequency brain networks across di erent number of layers. Connectivity changes are computed as the absolute value of the di erence between the links' weights of each AD network and those from the group-averaged HC brain network. Colored histograms correspond to di erences between replica (green), intralayer (blue) and interlayer links (yellow). The number of nodes (brain areas) stays the same, i.e., N . a) layers configuration. b) layers configuration. c) layers configuration. d) layers configuration. The layers (frequencies) are removed evenly on the range to Hz. . . . . . . Alzheimer's disease multilayer brain networks' global disruption in region-wise and frequencywise dimensions. a) Di erence between the average node multistrength centrality distribution of AD patients and HC in function of the average node multistrength centrality of HC. Regression line coe cient, i.e. , nodewise disruption index is Δ X = -0.44( 2 = 0.23). b) Di erence between the average layer multistrength centrality distribution of AD patients and HC in function of the average layer multistrength centrality of HC. Regression line coe cient, i.e., layerwise disruption index is Δ Y = -0.38( 2 = 0.78). c) Cognitive decline of AD patients (y-axis, MMSE) as a function of the individual nodewise disruption indices Δ X (x-axis). Spearman correlation test returns = -0.12, = 0.60 between Δ X and MMSE scores. d) Cognitive decline of AD patients (y-axis, MMSE) as a function of the individual layerwise disruption indices Δ Y (x-axis). Spearman correlation test returns = 0.63, = 1.9 -3 between Δ Y and MMSE scores. . . . . Alzheimer's disease multilayer brain networks' local disruption in region-wise and frequencywise dimensions. Multilayer brain networks are inferred from source-reconstructed MEG signals using cross-frequency coupling. In the primal dimension, nodes correspond to brain regions ( = 70) and layers to di erent frequency bins ( = 77). a) Statistical di erence (Wilcoxon test, Z-score) between node multistrength centralities of AD patients and healthy controls (HC) across di erent brain regions. b) Statistical di erence (Wilcoxon test, Z-score) between layer multistrength centralities of AD patients and healthy controls (HC) across di erent brain frequencies. c) Cognitive decline of AD patients (y-axis, MMSE) as a function of the multistrength centrality decrement in the left anterior caudal cingulate cortex (x-axis, X ). Spearman correlation test returns = 0.54, = 8 -3 between X and MMSE scores. d) Cognitive decline of AD patients (y-axis, MMSE) as a function of the multistrength centrality decrement at Hz (x-axis, Y ). Spearman correlation test returns = 0.61, = 2 -3 between Y and MMSE scores. Regressing curves resulting from a linear fit are shown for illustrative purposes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  and within band functional networks in resting-state MEG', [ ]: Tewarie et al. (

.

Figure . :

 . Figure . : Multiscale brain organization. The di erent organizational aspects of the brain system are represented over a multidimensional manifold. Three type of dimensions, or levels, are illustrated here, i.e. , time, space, and topology. From the top to the bottom of the manifold, the scales of each organizational level go from micro to macro. Image credit: Thibault Rolland. Adapted from [ ]

Figure

  Figure . : Illustration of a Barabási-Albert network (N ) and its node degree distribution.

Figure

  Figure . : The possible connected -nodes subnetworks with directed links. Adpated from [ ].

  [ ]: Blondel et al. ( ), 'Fast unfolding of communities in large networks', [ ]: Fortunato et al. ( ), 'Resolution limit in community detection', [ ]: Colizza et al. (

  undetectable [ ]. Community detection method based on bayesian inference were shown to overcome the limitations of modularity maximization [ ].

Figure . :

 . Figure . : Communities and rich clubs. (A) A network composed of four communities that are linked by hub nodes (black). (B) With the addition of further inter-module connections hub nodes now form a densely interconnected rich club, consisting of nodes with a degree of or higher.Figure and caption adapted from [ ].

  Figure . : Communities and rich clubs. (A) A network composed of four communities that are linked by hub nodes (black). (B) With the addition of further inter-module connections hub nodes now form a densely interconnected rich club, consisting of nodes with a degree of or higher.Figure and caption adapted from [ ].

Figure . :

 . Figure . : Illustration of Maslov-Sneppen rewiring procedure. The links of two pairs of nodes are reshu ed in a way to keep the node degree sequence fixed.

  Boccaletti et al. ( ), 'The structure and dynamics of multilayer networks', [ ]: Holme et al. (

Figure . :

 . Figure . : Main configurations of multilayer networks. Panel a) Full multilayer network. Both within-and between-layer connections are allowed with no specific restrictions. This configuration could typically be adopted to model multifrequency brain networks. Panel b) Multiplex network. Only interlayer connections between the replica nodes are allowed. No restrictions on connections within layers. This configuration is typically used to model real-world systems with multiple interactions such as multimodal brain networks. Panel c) Temporal network. Interlayer connections are allowed only between adjacent layers. No restrictions on connections within layers. This configuration is typically adopted to model time-varying networks. Image credit: Thibault Rolland.

[

  or based on tensorial formalisms [ ]. Formally, a multilayer network is defined as M = ( G, C) where G is a set of graphs and C a set of links connecting the nodes of the di erent graphs More precisely, G = { | ∈ ℕ} with = ( , ) being a

  . a). The supra-adjacency matrix of multilayer networks has the following form: = 11 12

  Yu et al. ( ), 'Selective impairment of hippocampus and posterior hub areas in Alzheimer's disease', are connected within the layers of the multiplex network (Fig. . ). In general the pair , is connected by a multilink -→ = ( 1 , 2 , ..., ).

Figure . :

 . Figure . : Example of all possible multilinks in a multiplex network with M layers and N nodes. Figure and caption adapted from [ ].

  PageRank centrality. The extension of betweenness centrality to multiplex networks is the so-called overlapping betweenness centrality which read [ ]: Rombach et al. ( ), 'Core-Periphery Structure in Networks', [ ]: Verma et al. ( ), 'Emergence of core-peripheries in networks', [ ]: Battiston et al. (

  based brain-computer interfaces: principles and applications', [ ]: Yu et al. ( ), 'Supplementary materials for Selective impairment of hippocampus and posterior hub areas in Alzheimer's disease: an MEG-based multiplex network study', [ ]: Tang et al. (

  -degree sequence.. In order to preserve the structural correlations across layers, one can fix the multidegree sequence [ ]. Similar to the Maslov-Sneppen algorithm, it consists in assigning -→ stubs to node . Next, pairs of stubs (of the same type -→ ) belonging to di erent nodes are matched. If nodes have matching stubs of type -→ , then they are connected with a multilink -→ (Fig. . ).

Figure

  Figure . : Schematic representation of the construction of a multiplex network with a given multi-degree sequence. Stubs of di erent types are associated with nodes. Solid line, dashed line and dot-dashed line stubs refer to multilinks ( , ), multilinks ( , ) and multilinks ( , ) respectively. Stubs of the same type are randomly matched, forming the multiplex network.Figure and caption adapted from [ ].

  Figure . : Schematic representation of the construction of a multiplex network with a given multi-degree sequence. Stubs of di erent types are associated with nodes. Solid line, dashed line and dot-dashed line stubs refer to multilinks ( , ), multilinks ( , ) and multilinks ( , ) respectively. Stubs of the same type are randomly matched, forming the multiplex network.Figure and caption adapted from [ ].

  model that produces communities by attributing a probability for nodes to randomly connect within and outside pre-defined communities. [ ]: De Domenico et al. ( ), 'Spectral Entropies as Information-Theoretic Tools for Complex Network Comparison', [ ]: Lamberti et al. ( ), 'Metric character of the quantum Jensen-Shannon divergence', [ ]: Iacovacci et al. ( ), 'Mesoscopic structures reveal the network between the layers of multiplex data sets', [ ]: Kao et al. ( ), 'Layer Communities in Multiplex Networks', [ ]: Tudisco et al. ( ), 'Node and Layer Eigenvector Centralities for Multiplex Networks',

Figure . :

 . Figure . : Layer clustering based on the entropy of network ensemble. a) The entropy

:

  A multi homogeneous map is a function f of several variables such that ( ) = ( ), with k the degree of homogeneity. [ ]: Rahmede et al. ( ), 'Centralities of nodes and influences of layers in large multiplex networks', [ ]: Zhang et al. ( ), 'Co-Ranking for nodes, layers and timestamps in multilayer temporal networks', [ ]: Taylor et al. (

  Since biological tissues are composed of high quantity of hydrogen nuclei, NMR is able to provide crucial information on their nature and organization in living beings. During a MRI experiment each D-volume of tissue is summarized in a small unit called voxel. The reconstruction of the studied objects(e.g., the brain) is made possible by repeating the NMR over each voxel. Diverse methods, each adapted to experiments or clinical conditions, are derived from the NMR principles among which T -and T -weighted MRI, di usion weighted imaging (DWI) and functional MRI (fMRI).

Figure

  Figure . : -Tesla MRI of a normal human brain. Left) T image. Right) T image. Adapted from [ ].

Figure . :

 . Figure . : Schematic of the di usion of a water molecule within the fiber bundles. Adapted from [ ]

Figure . :

 . Figure . : Illustation of the alternance of gyri and sulci at the surface of the cerebral cortex.

  others [ ]. fMRI has a poor temporal resolution ( ˜ . -second at best) and a voxel resolution of 3 3 . . . EEG and MEG Electroencephalography (EEG) and Magnetoencephalography (MEG), abbreviated by M/EEG, are recording technics of the brain electric and magnetic field, respectively. These fields are primarily rooted in intracellular and extracellular ionic currents existing at the neuron level. MEG signals arises more from intracellular currents whereas EEG arises from extracellular currents. More precisely, M/EEG signals record the oscillatory activity of a large population of synchronous neurons that are spatially localized. Oscillatory by nature, the M/EEG signals are characterized by their frequency, amplitude, phase, shape and their duration [ ].

Figure

  Figure . : EEG (up) and MEG (down) apparatus. Adapted from [ ].

  in ihren Prinzipien dargestellt auf Grund des Zellenbaues, : Composition of a biological tissue in terms of number of cells, their type and their spatial organization. [ ]: Amunts et al. ( ), 'Cytoarchitectonic mapping of the human amygdala, hippocampal region and entorhinal cortex', [ ]: Tzourio-Mazoyer et al. ( ), 'Automated Anatomical Labeling of Activations in SPM Using a Macroscopic Anatomical Parcellation of the MNI MRI Single-Subject Brain', [ ]: Desikan et al. ( ), 'An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest', [ ]: de Reus et al. ( ), 'The parcellation-based connectome: Limitations and extensions', [ ]: Thirion et al. ( ), 'Which fMRI clustering gives good brain parcellations?', : The most common deterministic tractography algorithms are FACT for Fiber Assigned by Continuous Tracking [ ] or the streamline tracking algorithm [].

  Figure . : Deterministic tractography. Seeds points are represented by white circles. They are the starting point to reconstruct, voxel-by-voxel, the principal direction of water di usion (white lines). Adapted from [ ].

  [ ]: Zalesky et al. ( ), 'On the use of correlation as a measure of network connectivity', [ ]: Smith et al. ( ), 'Network modelling methods for FMRI', [ ]: Bastos et al. ( ), 'A Tutorial Review of Functional Connectivity Analysis Methods and Their Interpretational Pitfalls', [ ]: Lachaux et al. ( ), 'Measuring phase synchrony in brain signals', [ ]: Nolte et al. ( ), 'Identifying true brain interaction from EEG data using the imaginary part of coherency', [ ]: Stam et al. ( ), 'Phase lag index: Assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources', [ ]: Muller et al. ( ), 'Cortical travelling waves: mechanisms and computational principles', [ ]: Ding et al. ( ), 'Granger Causality: Basic Theory and Application to Neuroscience', [ ]: Wang et al. ( ), 'Brain functional connectivity analysis using mutual information', [ ]: De Vico Fallani et al. ( ), 'Graph analysis of functional brain networks: practical issues in translational neuroscience', [ ]: Hero et al. ( ), 'Large-Scale Correlation Screening', [ ]: Palva et al. ( ), 'Ghost interactions in MEG/EEG source space: A note of caution on inter-areal coupling measures', [ ]: Kaiser ( ), 'A tutorial in connectome analysis: Topological and spatial features of brain networks', is to take the absolute value of negative correlation to obtain only positive weights on the network. Though it could be sensible to biases [ ], the Pearson correlation on fMRI signal was shown to perform well against other measures [ ].

  [ ]: Korhonen et al. ( ), 'Principles and open questions in functional brain network reconstruction', directly influenced by the density of the network [ ].

Figure . :

 . Figure .: Pipeline for functional brain networks modelling and analysis. Nodes correspond to specific brain sites according to the used neuroimaging technique. Links are estimated by measuring the FC between the activity of brain nodes; this information is contained in a connectivity matrix. By means of filtering procedures, based on thresholds, only the most important links constitute the brain network. The topology of the brain network is quantified by di erent network measures that can be represented as numbers (e.g., the colored bars). These measures can be the input to statistical analysis to look for significant di erences between populations/conditions (e.g., red points correspond to brain network measures of diseased patients or tasks, blue points stand for healthy subjects or resting states. Figure and caption adapted from [ ].

  Figure .: Pipeline for functional brain networks modelling and analysis. Nodes correspond to specific brain sites according to the used neuroimaging technique. Links are estimated by measuring the FC between the activity of brain nodes; this information is contained in a connectivity matrix. By means of filtering procedures, based on thresholds, only the most important links constitute the brain network. The topology of the brain network is quantified by di erent network measures that can be represented as numbers (e.g., the colored bars). These measures can be the input to statistical analysis to look for significant di erences between populations/conditions (e.g., red points correspond to brain network measures of diseased patients or tasks, blue points stand for healthy subjects or resting states. Figure and caption adapted from [ ].

  [ ]: Van Den Heuvel et al. ( ), 'Highcost, high-capacity backbone for global brain communication', [ ]: Bullmore et al. ( ), 'The economy of brain network organization', [ ]: Sporns ( ), 'Structure and function of complex brain networks', [ ]: Van Den Heuvel et al. ( ), 'Efficiency of Functional Brain Networks and Intellectual Performance', [ ]: Li et al. ( ), 'Brain Anatomical Network and Intelligence', [ ]: Bullmore et al. ( ), 'The economy of brain network organization', [ ]: Vértes et al. ( ), 'Simple models of human brain functional networks', [ ]: Hallett et al. ( ), 'Human brain connectivity: Clinical applications for clinical neurophysiology', [ ]: Stam ( ), 'Modern network science of neurological disorders', core regions integrate information coming from the peripheral ones [ ].

Figure

  Figure . : The di erent stages of Alzheimer's disease neurodegeneration.

  [ ]: Haan et al. ( ), 'Disruption of functional brain networks in Alzheimer's disease', [ ]: Stam et al. ( ), 'Graph theoretical analysis of magnetoencephalographic functional connectivity in Alzheimer's disease', [ ]: Haan et al. (

  . a)[ -].The thalamocortical model is perhaps one of the simplest examples that can reproduce disparate physiological and pathological conditions, from Parkinson's disease to epileptic seizures[ -]. In this model, both basic microscopic neurophysiology (e.g., synaptic and dendritic dynamics) and mesoscale brain anatomy (e.g., corticocortical and corticothalamic pathways) are progressively incorporated to predict large-scale brain electrical activity (Fig.. b).

  [ ]: Siettos et al. ( ), 'Multiscale modeling of brain dynamics', interactions give rise to collective behavior and to relationships with their environment is a central theme of modern multiscale brain modeling. Because of the substantial lack of biological evidence, especially concerning inter-scale connectivity, large parts of the studies have focused on analytical and numerical approaches [ ]. For example, intra-scale interactions have been simulated adopting cellular automata perspectives [ ], while inter-scale connectivity have been established using wavelet transformations [ ].

Figure . :

 . Figure . : Bottom-up hierarchical modeling. Panel a) The so-called K-set hierarchy showing the model progression from cell level to entire brain. K is a noninteracting collection of neurons. KI corresponds to a cortical column with su cient functional connection density. KII represents a collection of excitatory and inhibitory populations. KIII is formed by the interaction of several KII sets and simulates the known dynamics of sensory areas with /f spectra (see inset b). KIV is formed by the interaction of three KIII sets that models the genesis of simple forms of intentional behaviors. Panel b) Schematic view of major components involved in thalamocortical interactions. Di erent shading patterns code for di erent zones of the system, i.e. , from micro (relay nuclei, thalamic reticular nuclei (TRN)) to macro scales (cortex). As indicated by key, all connections shown are excitatory except for the connection from the reticular cell to the relay cell, which is inhibitory. Pictures and captions adapted from [ ] (panel a) and [ ] (panel b). Inset b adapted from [ ].

  data-driven multiscale modeling of brain circuits', [ ]: Sanz Leon et al. ( ), 'The Virtual Brain: a simulator of primate brain network dynamics', [ ]: Battiston et al. ( ), 'Structural measures for multiplex networks', [ ]: Battiston et al. ( ), 'Multiplex core-periphery organization of the human connectome', [ ]: Lim et al. ( ), 'Discordant attributes of structural and functional brain connectivity in a two-layer multiplex network', [ ]: Simas et al. ( ), 'An algebraic topological method for multimodal brain networks comparisons', [ ]: Guillon et al. ( ), 'Loss of brain inter-frequency hubs in Alzheimer's disease', [ ]: Yu et al. ( ), 'Selective impairment of hippocampus and posterior hub areas in Alzheimer's disease', [ ]: De Domenico et al. ( ), 'Mapping Multiplex Hubs in Human Functional Brain Networks', [ ]: Baillet et al. ( ), 'Electromagnetic brain mapping', [ ]: Michel et al. ( ), 'EEG source imaging', [ ]: Grech et al. (), 'Review on solving the inverse problem in EEG source analysis', ways of analyzing and interpreting them [ ]. This is actually one of the main challenges of big research projects supported by funding agencies around the world, such as the European Human Brain Project * or the US BRAIN Initiative . The increasing number of open-source tools that can be freely accessed and customized to enrich multiscale brain models just confirms how broad and multidisciplinary is the community e ort[ -].

  Domenico et al. ( ), 'Structural reducibility of multilayer networks', [ ]: De Domenico et al. ( ), 'Mapping Multiplex Hubs in Human Functional Brain Networks',

  . a)[ -]. This representation is particularly useful for functional brain networks with a broad frequency content, such as in those obtained from electrophysiology, EEG or MEG signals. Although less frequent than multiplexes, this type of representation has a great potential for characterizing whole brain cross-frequency coupling, which has been recently shown to be crucial for many cognitive and pathological mental states [ ].

Figure . :

 . Figure . : Structural reducibility of multifrequency brain networks. Panel a) For each combination of layers a quality function measures the amount of new information added with respect to an equivalent single-layer model. Panel b) Median values of quality function obtained from fMRI multifrequency brain networks in healthy subjects. Shaded areas indicate the standard deviation around each value. Pictures and captions adapted from [ ].

Figure . :

 . Figure . : Emergent properties in full multilayer brain networks Panel a) Intralayer and interlayer links in the multifrequency MEG network. -link between regions at the same frequency; -link of the same area between di erent frequency bands; -link between di erent nodes at di erent frequency bands. Panel b) Algebraic connectivity 2 as a function of the total interlayer connectivity ( ). The vertical solid line corresponds to the actual value of interlayer connectivity, i.e., without modifying their weights. Pictures and captions adapted from [ ].

Figure . :

 . Figure . : Multiplex motif analysis of multimodal brain networks. Panel a) Structural-functional -layer brain network. Interlayer links between replica nodes are omitted for the sake of visibility. Five nontrivial multiplex motifs of two nodes are possible based on the type of connectivity in the DTI structural layer and in the fMRI functional layer. The Z-scores show the motifs that are overrepresented and underrepresented as compared to equivalent random networks. Panel b) Patterns of multiplex triangles comprising directed structural tuples (solid connections) closed by a functional link (dashed connections). The overall motif counts normalized by equivalent random multiplexes are illustrated as a function of basal activation parameters P and Q of the Wilson-Cowan model. Pictures and captions are adapted from [ ], (panel a), and [ ] (panel b).

  Jirsa et al. ( ), 'Cross-frequency coupling in real and virtual brain networks', [ ]: Garcés et al. ( ), 'Multimodal description of whole brain connectivity: A comparison of resting state MEG, fMRI, and DWI', to oscillatory dynamics, as well as in other regions of the parameter space (Fig. .

Figure . :

 . Figure . : Multiplex core-periphery structure of the human connectome. Scatter plot of multiplex coreness against single-layer corenesses obtained from structural (DTI) and functional (fMRI) layers. Labels indicate brain areas whose multiplex coreness cannot be predicted by looking at the coreness values in the respective structural and functional layer. Picture and caption adapted from [ ].

Figure . :

 . Figure . : Temporal network flexibility correlates with brain perfomance. Panel a) An overview of network switching (or flexibility) in a temporal network. Red and blue colors identify the nodes belonging to two di erent communities according to the multilayer network modularity metric. Panel b) Brain maps of switching rate and dynamic fMRI connectivity. Values were normalized into z-scores to ensure both connectivity dynamics and switching values were scaled equally. Pictures and captions adapted from [ ].

Figure . :

 . Figure . : Stabilization of critical dynamics in multilayer glia-neuronal networks. Panel a) Left side: Glia cells redistribute metabolic resources from the bloodstream to neural synapses. Right side: Associated two-layer network model. Black arrows indicate neural synaptic interactions. Arrow thickness indicates synaptic strength which evolves according to spike timedependent plasticity (STDP). Red arrows which terminate on black arrows represent the resource supply to the corresponding synapse. Panel b) Stability analysis of the two-layer STDP model. The largest eigenvalue of the neuronal network layer and the total resource R of all glia and synapses are illustrated as a function of time. The data plotted in black correspond to a 'baseline' condition. For the data plotted in red (labelled 'instability'), the initial evolution is the same as for the baseline data up until the di usion of resources between the glial cells is turned o (vertical arrow). Pictures and captions adapted from [ ].

Figure . :

 . Figure . : Temporal network flexibility predicts future learning rate. Significant predictive Spearman correlations between flexibility in session and learning in session (black curve,≈ 0.001) and between flexibility in session and learning in session ( ≈ 0.009). Each point corresponds to a subject. Note that relationships between learning and fMRI network flexibility in the same experimental sessions ( and ) were not significant; > 0.13 was obtained using permutation tests. Picture and caption adapted from [ ].

  hubs in the human brain', [ ]: Uhlhaas ( ), 'Dysconnectivity, large-scale networks and neuronal dynamics in schizophrenia', [ ]: Uhlhaas et al. (

Figure . :

 . Figure . : Temporal network flexibility as a clinical marker of shizophrenia genetic risk. Panel a) Significant increases in the mean dynamic reconfiguration of modular fMRI brain networks in una ected first-grade relatives (gray bar, REL) and patients with schizophrenia (black bar, SZ) in comparison to matched healthy controls (white bar, HC) [F( , ) . , P . ]. Bars indicate mean values, and whiskers represent standard error means (SEMs). Panel b) Significant increases in the mean dynamic reconfiguration of modular brain networks in healthy controls after application of dextrometorphan (DXM) [dark gray bars; repeated measures ANOVA placebo (PLA) versus DXM: F( , ) . , P . ] relative to PLA (light gray bars). Pictures and captions adapted from [ ].

  Multiplex networks have been used to provide a unified description of AD brain reorganization across multiple MEG frequency bands (Fig. . a). [ ] used di erent multiplex nodal metrics (e.g., overlapping clustering, local e ciency and betweenness centrality) and consistently showed that physiological multilayer hub regions, including posterior parts of the DMN, were severely impacted by AD (Fig. . b).

  core-periphery structure of multimodal brain networks in Alzheimer's disease', [ ]: Guillon et al. ( ), 'Disrupted core-periphery structure of multimodal brain networks in Alzheimer's disease', [ ]: Chételat et al. ( ), 'Relationship between atrophy and -amyloid deposition in Alzheimer disease', [ ]: Yu et al. ( ), 'Selective impairment of hippocampus and posterior hub areas in Alzheimer's disease', [ ]: Canal-Garcia et al. ( ), 'Multiplex Connectome Changes across the Alzheimer's Disease Spectrum Using Gray Matter and Amyloid Data', [ ]: Guillon et al. ( ), 'Disrupted core-periphery structure of multimodal brain networks in Alzheimer's disease', [ ]: Kubicki et al. ( ), 'The Early Indicators of Functional Decrease in Mild Cognitive Impairment', [ ]: Albers et al. ( ), 'At the interface of sensory and motor dysfunctions and Alzheimer's disease',

Figure . :

 . Figure . : Multifrequency and temporal reorganization of brain networks in Alzheimer's disease. Panel a) Multiplex brain networks are constructed by layering di erent frequency-specific networks, while temporal networks were constructed by concatenating time-specific networks within frequency bands. Panel b) Top: hub disruption of MEG multifrequency networks in patients with Alzheimer's disease. Each point corresponds to a di erent brain area; k slope of the regressing line. Bottom: Brain regions with significant between-group di erence in overlapping weigthed degree. PCUN.R right precuneus; HIP.L left hippocampus; IPL.R right inferior parietal, but supramarginal and angular gyri; SPG.R right superior parietal gyrus; MOG.L left middle occipital gyrus; SOG.L left superior occipital gyrus; IOG.L left inferior occipital gyrus. Panel c) Scatter plot showing the Mahalanobis distance of each subject from AD or control group when combining of multiplex clustering coe cient (MCC) and participation coe (MPC) extracted from time-varying networks (gray line indicates equal distance). Pictures and captions adapted from [ ] (panel a, b) and from [ ] (panel b).

Figure

  Figure . : Multimodal brain networks reveal disrupted core-periphery structure in Alzheimer's disease. Panel a) Multimodal brain networks (multiplex) are constructed by layering DTI, fMRI, and several frequency-based MEG brain connectivity. Panel b) Spearman correlation (R . , p . ) between the coreness disruption index ( ) and the memory impairment of AD patients as measured by the free recall (FR) test. Panel c) Boxplots show the values of coreness disruption index ( ) obtained by progressively removing the links preferentially connected to the multiplex periphery of the HC group. The blue (x-axis HC) and red (x-axis AD) boxplots illustrate respectively the values for the HC and AD groups. Pictures and captions adapted from [ ] (Panel a) and [ ] (Panel b,c).

Figure . :

 . Figure . : Topological duality in multilayer networks.The unitary element of a multilayer network is the duplet ( , ) corresponding to the node and the layer . Because of the intrinsic node-layer isomorphism, the system can be equivalently represented by nodes plunged into layers or by layers plunged into nodes. Depending on the representation side, two complementary descriptions of the same system can be obtained. In the primal nodewise, connectivity is integrated across the layers and the topological properties of the nodes are provided. In the dual layerwise, connectivity is integrated across the nodes and the topological properties of the layers are obtained instead. In the networks, solid lines intralayer links, grey lines interlayer links, dashed lines replica links.

  [ ]: De Domenico et al. ( ), 'Structural reducibility of multilayer networks', . . Primal and dual descriptions are independent at the 'second order'

Figure . :

 . Figure . : Nodewise and layerwise distributions are non-constrained from each other. The MC distributions are normalized in order to have the same mean. a) Nodewise (blue) and layerwise (orange) MC distribution of the EuAir multiplex network [ ]. The normalized nodewise standard deviation is X = 0.73. The normalized layerwise standard deviation is Y = 0.45.b) Eu-Air multiplex network is rewired in a way that layerwise degree distribution stays fixed and the nodewise one is modified using the stochastic rewiring algorithm ( = 1, = 1, see Sec. . . ). The normalized standard deviations become X = 0.10, Y = 0.45.c) EuAir multiplex network is rewired in a way that nodewise degree distribution stays fixed and the layerwise one is modified using the stochastic rewiring algorithm ( = 0.66, = 1, see Sec.

  multiplex and multilayer networks (see definitions in Sec. .

Figure . :

 . Figure . : E ect of the density on the matching between analytical and numerical distances. Curves are the analytical distances as obtained with Eqs. . , . for the layer multidegre centrality (MC) and Eqs. . , . for the node multidegree centrality.Dots are the numerical distances obtained with the stochastic rewiring algorithm (see Appendix ). The numerical expected MC is obtained by rewiring of the initial networks. Error bars represents times the standard error of the mean. a) Nodewise distance in function of r for random multilayer networks ( = 20, = 20) at di erent densities ( = 1). b) Layerwise distance in function of r for random multilayer networks ( = 20, = 20) at di erent densities ( = 1). c) Nodewise distance in function of r for random multilayer networks ( = 20, = 20) at di erent densities ( = 1). d) Layerwise distance in function of r for random multilayer networks.

Figure . :

 . Figure . : Linear relation between layerwise and nodewisee multidegree centrality distances in random multilayer networks. Lower slopes (higher X ) are obtained for > . Higher slopes (higher Y ) are obtained for > ). Solid lines correspond to Eqs. . , . in the case of random multilayer networks with = = 200, = 0.0005, by varying the entire range of possible rewiring (color line). Scattered points correspond to synthetic random networks simulated with the same parameters (one synthetic random network di erent for each point).

Figure . :

 . Figure . : Linear relation between layerwise and nodewisee multidegree centrality distances in random multiplex networks. Lower slopes (higher X ) are obtained for > . Higher slopes (higher Y ) are obtained for > ). Solid lines correspond to Eqs. . , . in the case of random networks with = = 200, = 0.01 by varying the entire range of possible rewiring (color line). Scattered points correspond to synthetic random mutiplex networks simulated with the same parameters (one synthetic random network di erent for each point).

  Domenico et al. ( ), 'Structural reducibility of multilayer networks',

Figure . :

 . Figure . : Dual characterization of real-world multiplex networks.Scatter plot of the nodewise (x-axis) and layerwise (y-axis) distances of real networks' MCs from uniformly rewired counterparts. To avoid network-size and density biases, all values are further divided by the distances obtained from equivalent random networks. The logarithm of the aforementioned quantities are shown in the plot. The optimal partition through k-mean clustering separates the network into two clusters (green and purple). Zones indicate in which cluster the networks belong to.

  ]. These clustering results are confirmed by a hierarchical clustering algorithm (Fig..

Figure . :

 . Figure . : Optimal partition of the dual characterization of real-world multiplex networks (k-means clustering). a) Elbow method. Sum of the squared error in function of the number of clusters. The elbow point corresponds to the partition of space into two clusters which is considered as the optimum trade-o between error and the number of clusters (with knee-point detected using python package kneed). b) The silhouette coe cient in function of the number of clusters. The bigger is the silhouette coe cient, the more cohesive and separated the clusters. Its values ranges from -to . Here the maximum is obtained for two clusters, = 0.70.

Figure . :

 . Figure . : Clusters of real-world multiplex networks using hierarchical clustering. Colors of each network indicate in which cluster they belong to (green or purple).

Figure . :

 . Figure . : Projection plots of real-world multiplex networks. . In the nodewise description, layers are represented by lines and nodes are represented by markers (blue). In the layerwise description, layers are represented by markers (orange) and nodes by lines. Nodewise and layerwise plots are put on the same scale. a) Projection plots of the C.Elegans multiplex connectome ( nodes, layers). b) Projection plots of the Arxiv multiplex network ( nodes, layers). c)Projection plots of the Chemnitz multiplex network in the German Transport dataset ( nodes, layers). d) Projection plots of the PierreAuger multiplex network ( nodes, layers). Inset in the nodewise projection plot represents a zoom where markers are concentrated.

  izing interactions in online social networks during exceptional events', [ ]: De Domenico et al. ( ), 'Unraveling the Origin of Social Bursts in Collective Attention',

Figure . :

 . Figure . : Number of nodes in function of the number of layers for real-world multiplex networks. Size of the markers is proportional to the logarithm of the density of each multiplex network. Purple indicates networks with a strong spatial connotation and green, networks with no or weak spatial correlation.

Figure . :Figure

 . Figure . : Distribution of the total multistrength centrality of AD patients compared to the group-averaged HC. The total multistrength centralities of AD patients (orange bars) are compared to the one of the group-averaged HC(blue line). For every tested layer configuration, the ratio between the total multistrength centrality of the group-averaged AD patient over the one of the group-averaged HC is constant equal to 1.06. All AD patients except outliers have a total multistrength centrality that is very close to the one of the group-averaged HC. a) layers configuration. b) layers configuration. c) layers configuration. d) layers configuration. The layers (frequencies) are removed evenly on the range to Hz.

Figure . :

 . Figure . : Distribution of connectivity changes in AD multifrequency brain networks across di erent number of layers. Connectivity changes are computed as the absolute value of the di erence between the links' weights of each AD network and those from the group-averaged HC brain network. Colored histograms correspond to di erences between replica (green), intralayer (blue) and interlayer links (yellow). The number of nodes (brain areas) stays the same, i.e., N . a) layers configuration. b) layers configuration. c) layers configuration. d) layers configuration. The layers (frequencies) are removed evenly on the range to Hz.

Figure :

 : Figure : Alzheimer's disease multilayer brain networks' global disruption in region-wise and frequency-wise dimensions. a) Di erence between the average node multistrength centrality distribution of AD patients and HC in function of the average node multistrength centrality of HC. Regression line coe cient, i.e. , nodewise disruption index is Δ X = -0.44( 2 = 0.23). b) Di erence between the average layer multistrength centrality distribution of AD patients and HC in function of the average layer multistrength centrality of HC. Regression line coe cient, i.e., layerwise disruption index is Δ Y = -0.38( 2 = 0.78). c) Cognitive decline of AD patients (y-axis, MMSE) as a function of the individual nodewise disruption indices Δ X (x-axis). Spearman correlation test returns = -0.12, = 0.60 between Δ X and MMSE scores. d) Cognitive decline of AD patients (y-axis, MMSE) as a function of the individual layerwise disruption indices Δ Y (x-axis). Spearman correlation test returns = 0.63, = 1.9 -3 between Δ Y and MMSE scores.

Figure . :

 . Figure . : Alzheimer's disease multilayer brain networks' local disruption in region-wise and frequency-wise dimensions. Multilayer brain networks are inferred from source-reconstructed MEG signals using cross-frequency coupling. In the primal dimension, nodes correspond to brain regions ( = 70) and layers to di erent frequency bins ( = 77). a) Statistical di erence (Wilcoxon test, Z-score) between node multistrength centralities of AD patients and healthy controls (HC) across di erent brain regions. b) Statistical di erence (Wilcoxon test, Z-score) between layer multistrength centralities of AD patients and healthy controls (HC) across di erent brain frequencies. c) Cognitive decline of AD patients (y-axis, MMSE) as a function of the multistrength centrality decrement in the left anterior caudal cingulate cortex (x-axis, X ). Spearman correlation test returns = 0.54, = 8 -3 between X and MMSE scores. d) Cognitive decline of AD patients (y-axis, MMSE) as a function of the multistrength centrality decrement at Hz (x-axis, Y ). Spearman correlation test returns = 0.61, = 2 -3 between Y and MMSE scores. Regressing curves resulting from a linear fit are shown for illustrative purposes.

  and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates', [ ]: Daducci et al. ( ), 'The Connectome Mapper: An Open-Source Processing Pipeline to Map Connectomes with MRI', preprocessed epochs per participant. The source reconstruction of the MEG activity was performed in two steps [ , ]. First, to solve the direct problem, we used the individual segmented T -weighted images [ , ] from which we imported cortical surfaces in the Brainstorm software [ ]. The software was used to model approximately , equivalent current dipoles. Then, to solve the inverse problem, we applied the weighted minimum norm estimate approach with overlapping spheres [ ]. . . Functional connectivity Nodes of the MEG-based multifrequency brain networks are Regions of Interests (ROIs), N , defined by the Lausanne cortical atlas parcellation (subcortical ROIs not included) [ ]. Layers are frequency bins ranging from Hz to

  present the stochastic rewiring algorithm in Python. The code with mock data is available on GitHub (link for the online version). import numpy as np import pickle from random import choices def to_ supra_ index(N,uplet): """Convert tensorial indices (layer,layer,node,node) into flattened matrix indices (row,col) Input: N -number of nodes uplet -tensor indices (layer,layer,node,node) Output: row, col -matrix indices (row,col)""" k, l, i, j uplet[ ],uplet[ ],uplet[ ],uplet[ ] row k*N i col l*N j return row,col def to_ tensor_ index(N,uplet): """Convert flattened matrix indices (row,col) into tensorial indices (layer,layer,node,node) Input: N -number of nodes uplet -matrix indices (row,col) Output: row, col -tensor indices (layer,layer,node,node)""" row,col uplet[ ],uplet[ ] k row // N l col // N i row N j col N return k,l,i,jdef build_ instructions(uplet,com): """Set which type of rewiring are allowed or not with a boolean instruction ex: com["be_ interlayer] False forbids the rewiring algorithm to rewire links in an interlayer space Input: uplet -proposed tensor indices (layer,layer,node,node) rewire in com -dictionary with two keys ["be _ interlayer"] and ["replica"] Output: boolean -False if the proposed rewiring is forbidden """ a,b,c,d uplet[ ], uplet[ ], uplet[ ], uplet[ ] if com["be_ interlayer"]False: CARE: present a matrix that has no interlayer links otherwise it runs during infinity if a! b: if we are in an interlayer space for the new_ uplet... return False ... then forbids this change by False elif com["replica"] False: if (a! b and c d): we selected a replica as new_ uplet return False the change is forbidden return True def select_ intra(A,N,uplet,rng,com,trial): """Return the hole where 'uplet' is rewired according to the p_ node event, which keeps the (inter)layer indices fixed Input: A -supra-adjacency matrix N -number of nodes uplet -tensor indices (layer,layer,node,node) of the link to rewire rng -random number generator com -dictionary with two keys ["be interlayer"] and ["replica"] trial -number of proposed new rewired position Output: [ , , , ] -Position where the link "uplet" is rewired, it links node a in layer k to node b in layer l """ k, l, i, j uplet[ ], uplet[ ], uplet[ ], uplet[ ] the coordinates of the selected uplet rand_ N list(rng.integers(N,size [ , 2])) for a,b in rand_ N: avoid_ self_ loops not((a b) and (k l)) condition to avoid rewiring in self-loops conditions build_ instructions([ , , , ],com) and avoid_ self_ loops and avoid_ rewire_ inter and avoid_ double_ replica row,col to_ supra_ index(N,[ , , , ]) if A[ , ]

  ["replica"] trial -number of proposed new rewired position Output: [ , , , ] -Position where the link 'uplet' is rewired, it links node c in layer a to node d in layer b """ k, l, i, juplet[ ], uplet[ ], uplet[ ], uplet[ ] the coordinates of the selected uplet rand_ M list(rng.integers(M, size [trial, ])) rand_ N list(rng.integers(N, size [trial, ])) rand_ NM [[rand_ 1[ ][0],rand_ 1[ ][1],rand_ 1[ ][0],rand_-1[ ][1]] for i in range(len(rand_ M ))]for a,b,c,d in rand_ NM: teleport_ condition ((a ! k or b ! l) and (c! i or d! j)) or ((a ! l or b ! k) and (c! j or d! i)) minimal condition for the teleport event teleport_ condition_ ((a ! k or b ! l) and (c! j or d! i)) or ((a ! l or b ! k) and (c! i or d! j)) minimal condition for the teleport event if teleport_ condition and teleport_ condition_ : avoid_ self_ loops not ((d c) and (a b)) condition to avoid rewiring in self-loops avoid_ rewire_ inter not ((c i and d j) or (c j and d i)) condition to avoid making p_ layer events row, col to_ supra_ index(N, [ , , , ]) conditions instructions([ , , , ], com) and avoid_ self_ loops and avoid_ rewire_ inter if A

  link were rewired already random_ index rng.integers( ,len(row), size ) classic_ uplet [ [ _ ][0], [ _ ][0]] always the first one of the random list uplet to_ tensor_ index(N,classic_ uplet) new_ uplet None security

  uplet[ ], uplet[ ], uplet[ ], uplet[ ] u, v, w, x new_ uplet[ ], new_ uplet[ ], new_ uplet[ ], new_ uplet[ ] Update the multidegree sequences with the rewiring k,l,i,j -u,v,w,x l_ Ktranspose[k] -l_ Ktranspose[l] -l_ K[i] -l_ K[j] -l_ Ktranspose[u] l_ Ktranspose[v] l_ K[x] l_ K[w] classic_ new_uplet to_ supra_ index(N,new_ uplet) Update the supra-adjacency matrix with the rewiring k,l,i,j -u,v,w,x s_ B[classic_ new_ uplet[ ],classic_ new_ uplet[ ]] indicates to the function that this links was rewired already s_ B[classic_ new_ uplet[ ], classic_ new_ uplet[ ]] s_ B[classic_ uplet[ ], classic_ uplet[ ]] s_ B[classic_ uplet[ ], classic_ uplet[ ]] return s_ B,l_ K,l_ Ktranspose,number_ of_ unrewire_ edge if _ _ name_ _ "_ _ main_ _ ": Parameters print("Begin") name_ path 'mock_ full_ multilayer_ N_ _ M_ .dat' save_ path "//" with open(name_ path,'rb') as f: loading the data parameter pickle.load(f) s_ A parameter['sparse'] supra-adjacency amtrix of the network to be rewired multi_ degree parameter['l_ K'] node multidegree sequence multi_ degree_ DS parameter['l_ Ktranspose'] layer multidegree sequence R proportion of edges to rewire number_ of_ trial_ per_ event for each trial draw are tested

  Network Science . Classical networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Basics and concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . Measures on networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Null models in network science . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Network comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Multilayer networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Types of multilayer networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Multilayer network measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Null models of multilayer networks . . . . . . . . . . . . . . . . . . . . . . . . . . . .

-Renyi model. The simplest

	[ ]: Barabàsi, Albert-László (	), Net-
	work science,
	[ ]: Maslov et al. (	), 'Specificity
	and Stability in Topology of Protein Net-
	works',
	: It exists also the soft configuration
	model that generates an ensemble of	model to generate random network
	networks with fixed expected degree is the Erdos-Renyi model [ ]. In its more widely used variation [ ], sequence.
	the	, model, it is defined by two parameters: the number of nodes

  , these results demonstrated the clinical value of considering multiplex/multilayer network approaches to derive more reliable neuromarkers of consciousness disorders.Alzheimer's disease (AD) is a neurodegenerative disorder and the most common form of dementia. Clinically, it is characterized by mild memory impairments that gradually evolve up to severe cognitive impairments and eventually to death. In 2016, people a ected by AD and other dementias were around 44 million worldwide and this incidence is likely to augment because of longer life expectancy [ ]. At cellular level, AD is characterized by the progressive accumulation of -tangles and -amyloid plaques that cause neurons and synapses to die, thus leading to brain atrophy and disordered dysconnection patterns.

	. Multilayer network characterization of
	Alzheimer's disease
	[ ]: Nichols et al. (	), 'Global,
	regional, and national burden of
	Alzheimer's disease and other demen-
	tias,	-	: a systematic analysis
	for the Global Burden of Disease Study
		',
	[ ]: Tijms et al. (	), 'Alzheimer's dis-
	ease: connecting findings from graph
	theoretical studies of brain networks',
	[	]: Gaubert et al. (	), 'EEG evi-
	dence of compensatory mechanisms
	in preclinical Alzheimer's disease',
	[ ]: Yu et al. (	), 'Selective impair-
	ment of hippocampus and posterior
	hub areas in Alzheimer's disease',
	[ ]: Guillon et al. (	), 'Loss of brain
	inter-frequency hubs in Alzheimer's dis-
	ease',
	[	]: Echegoyen et al. (	), 'From sin-
	gle layer to multilayer networks in mild
	cognitive impairment and Alzheimer's
	disease',
	[	]: Cai et al. (	), 'Functional Inte-
	gration and Segregation in Multiplex
	Brain Networks for Alzheimer's Dis-
	ease',

liminary

  . )

	Figure . : Invisibility zones for mul-
	tidegree centrality distances. Node-
	wise distances ( X ) increase linearly
	with	(x-axis) and	(white diago-
	nals) but they cannot see link displace-
	ments that keep nodes unchanged
	(	→ 1). Layerwise distances ( Y )
	increase linearly with	(y-axis)
	and	(white diagonals) but they are
	blind to link displacements that keep
	layers unchanged (	→ 1).

Table . :

 . Average

	nodewise and layer-wise standard deviations of the multi-degree centrality distribution for real-	Classes German Transport	< X > .	< Y > .
	world multiplex networks.	Arxiv	.
		PierreAuger	.
		Uganda villages	.	.
		Twitter events	.
		HumanMicrobiome	.	.
		C.Elegans	.	.
		FAO trade
		EuAir	.	.
		Genetic	.

  The study involved Alzheimer's disease (AD) patients ( women) and healthy agematched control (HC) subjects ( women). All participants underwent the Mini-Mental State Examination (MMSE) for global cognition. Inclusion criteria for all participants were: age between and . absence of general evolutive pathology. no previous history of psychiatric diseases. French as a mother tongue. clinical diagnosis of Alzheimer's disease. Mini-Mental State Examination (MMSE) score greater or equal to .

  before iplet is declared None event_ proba [ , , ] [ _ , _ , _ ] giving the probability of each related rewiring event rng np.random.default_ rng() dic_ instructions dic_ instructions["be interlayer"] True if False, we dont rewire into interlayer space, False is needed for multiplex networks) dic_ instructions["replica"] True if False, we can't rewire nodes in replica Running rewiring algorithm s_ B,l_ K,l_ Ktranspose,number_ unrewire stochastic_ rewiring(s_-A,multi_ degree,multi_ degree_ DS, event_ proba,R,dic_ instructions,number_-of_ trial_ per_ event number_ of_ trial_ per_ event) Rosa Crum et al. 'Population-Based Norms for the Mini-Mental State Examination by Age and Educational Level'. In: JAMA . (May ). _eprint: https://jamanetwork.com/journals/jama/articlepdf/ Baillet et al. 'Evaluation of inverse methods and head models for EEG source localization using a human skull phantom'. In: Physics in Medicine and Biology . (Jan. Hsuan Lin et al. 'Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates'. en. In: NeuroImage . Lucas Lacasa et al. 'From time series to complex networks: The visibility graph'. en. In: Proceedings of the National Academy of Sciences .

	[ ] /jam
		_ _	.pdf, pp.	-		.	: .	/jama.	.	(cited on page ).
	[ ] S ), pp. -.	:
		.	/	-	/ / /	. (Visited on / /	) (cited on page	).
	[ ] B. He. 'Brain electric source imaging: scalp Laplacian mapping and cortical imaging'. eng. In:
		Critical Reviews in Biomedical Engineering . -(	), pp.	-	(cited on page	).
	[ ] Bruce Fischl et al. 'Whole Brain Segmentation'. en. In: Neuron . (Jan.	), pp.	-.	:
		.	/S	-	( )	-X. (Visited on / /	) (cited on page	).
	[	] Bruce Fischl et al. 'Sequence-independent segmentation of magnetic resonance images'. en. In:
		NeuroImage (Jan.	), S -S .	: .	/j.neuroimage.	. .	. (Visited on
		/ /	) (cited on page	).
	Saving data [ ] François Tadel et al. 'Brainstorm: A User-Friendly Application for MEG/EEG Analysis'. en. In: Com-
		dico putational Intelligence and Neuroscience	(	), pp. -.	:	.	/	/	.
		(Visited on / /	dico[ ) (cited on page	] s_ B ).
	dico[ _ ] l_ K dico[ _ dico[ [ ] Fa-(May ] l_ Ktranspose ] event_ proba dico[ ] len(l_ K) . /j.neuroimage. . . . (Visited on / / ) (cited on page ), pp.	-. ).	:
						dico[ dico[	] len(l_ Ktranspose) ] s_ A.nnz	:	.	/journal.pone.
			(cited on page	dico[ ] R ).
						dico[	_	] event_ proba
		Physical Review	dico[ . (Aug. ), pp.	_ _ -.	:	.	_	] number_ unrewire /PhysRev. . . (Visited on
		/ /	) (cited on page	with open(save_ path "EuAir_ idx_ " str(idx), 'wb') as f: ).
							pickle.dump(dico, f)
		(Apr.	), pp.	-.		: .	/RevModPhys. .	. (Visited on / /	) (cited on
		page	).
	[ ] Franz Kaiser et al. 'Dual communities in spatial networks'. In: Nature Communications . (Dec.
		), p.	.	: .			/s	-	-	-(cited on page	).
	[ ] Dmitri Krioukov and Massimo Ostilli. 'Duality between equilibrium and growing networks'. en. In:
		Physical Review E . (Aug.	), p.	.	: .	/PhysRevE. .	. (Visited on
		/ /	) (cited on page	).
	[ ] Charles Murphy et al. Duality between predictability and reconstructability in complex systems.
		arXiv:	.	[cond-mat, physics:physics]. Feb.	.	: http://arxiv.org/abs/	.
		(visited on / /		) (cited on page	).
	[ ] (Apr.	), pp.	-	.	:	.	/pnas.
			. (Visited on / /	) (cited on page	).
	[	] Andriana S. L. O. Campanharo et al. 'Duality between Time Series and Networks'. In: PLOS ONE .
		(Aug.	), pp. -.	: .	/journal.pone.	(cited on page	).

[ ] Alessandro Daducci et al. 'The Connectome Mapper: An Open-Source Processing Pipeline to Map Connectomes with MRI'. In: PLOS ONE . (Dec. ), pp. -. [ ] H. A. Kramers and G. H. Wannier. 'Statistics of the Two-Dimensional Ferromagnet. Part I'. en. In: [ ] Robert Savit. 'Duality in field theory and statistical systems'. en. In: Reviews of Modern Physics .

Acknowledgements

I would like to express my deepest appreciation to my supervisor, Dr.Fabrizio De Vico Fallani that gave me the opportunity to prepare my PhD thesis at his side. By developing my autonomy, by supporting me, by helping me framing my ideas with texts, presentations and of course pictures, he taught me the job of research scientist.

. Stochastic rewiring model

In this section, we first present the model that allows us to change the local connectivity of a multilayer system. Then, we derive anaytically the expected multidegree centrality (MC) sequence in the nodewise and layerwise description according to the parameters of the model. Next, we propose a MC-based Euclidean distance to summarize the nodewise and layerwise changes in the general case. In the case of random multilayer and random multiplex networks, we provide a closed-form of these distances establishing a formal duality between nodewise and layerwise description. Finally, we elaborate on properties of these distances according to the parameters.

. . Notations

Latin letters encode for node indices and greek letters for layer indices in the nodewise description.

is the number of nodes in the nodewise description (the classical one).

: number of layers in the nodewise description. X : node MC sequence. Y : layer MC sequence. : total number of unique links in the system. : total possible number of unique links in the system. = ( -1) 2

for multiplex networks with no self-loops. ( ): weight of link between node i in layer and node j in layer after rewiring. For binary networks it equals or .

(0) is the related weights before any rewiring.

Π

, ( ): probability to select the link connecting a in layer to node b in layer at the -th rewiring.

, → , ( ): transition probability to move the link connecting a in layer to node b in layer to link connecting i in layer and j in layer at the -th rewiring .

( = ): probability to select a given rewiring rule. It can take the value = { , , } associated with the probabilities , , respectively. Note that + + = 1.

.

. Stochastic rewiring model

The model selects links uniformly at random in a given multilayer network. For each selected link a new position is randomly drawn based on the relative importance of the rewiring parameters , , and , which determine in a probabilistic fashion the type of link reassignment (Fig. . ): . New node indices are drawn uniformaly at random from all possible node indices. Layer indices stay fixed.

. New layer indices are drawn uniformaly at random from all possible layer indices. Node indices stay fixed.

the manuscript (e.g. , ( ) = ( )):

Note that we forbid the link to be rewired in the exact same place as it is. Once selected, the link connecting (i, ) to (j, ) necessarily goes somewhere in the multilayer network i.e, ( , )≠( , )

, → , = 1. The transition probability is conditionned by the selected rewiring rule whose probability is ( = ). Using ( = ) = 1 the previous formula read:

Each rule forces the link to be rewired in a particular way, e.g, for = , the link connecting (i, ) in (j, ) is necessarily rewired in a position that conserves the layer indices (see below). The transition probability conditionned to the rewiring rules read: is the total possible number of unique links in the (inter)layer , , is the total possible number of unique links for the pair i,j. Since we deal with multilayer networks composed of replica nodes, the possible number of unique links for a (inter)layer and for a pair of nodes does not depend on the node and layer indices. Consequently, we note = and = . For convenience we also write -Δ = --+ 1. By expliciting Eq. . , using Eqs. . , . ,we obtain:

( . )

: Spearman correlation between the AD patients' frequency multistrength centrality and the MMSE clinical score. Only the most significant correlations are reported here ( < 0.05).

In bold, we show the frequencies whose significance is lower than the mean false discovery rate (mFDR) for independent tests ( < 0.0253). Nine frequencies in the alpha band are close to the statistical significance after correction with the Benjamini-Hochberg method.

Here are references in citation order.

[