Efficient exploration of multiplex networks

Efficient techniques to navigate networks with local information are fundamental to sample large-scale online social systems and to retrieve resources in peer-to-peer systems. Biased random walks, i.e. walks whose motion is biased on properties of neighbouring nodes, have been largely exploited to design smart local strategies to explore a network, for instance by constructing maximally mixing trajectories or by allowing an almost uniform sampling of the nodes. Here we introduce and study biased random walks on multiplex networks, graphs where the nodes are related through different types of links organised in distinct and interacting layers, and we provide analytical solutions for their long-time properties, including the stationary occupation probability distribution and the entropy rate. We focus on degree-biased random walks and distinguish between two classes of walks, namely those whose transition probability depends on a number of parameters which is extensive in the number of layers, and those whose motion depends on intrinsically multiplex properties of the neighbouring nodes. We analyse the effect of the structure of the multiplex network on the steady-state behaviour of the walkers, and we find that heterogeneous degree distributions as well as the presence of inter-layer degree correlations and edge overlap determine the extent to which a multiplex can be efficiently explored by a biased walk. Finally we show that, in real-world multiplex transportation networks, the trade-off between efficient navigation and resilience to link failure has resulted into systems whose diffusion properties are qualitatively different from those of appropriately randomised multiplex graphs. This fact suggests that multiplexity is an important ingredient to include in the modelling of real-world systems.

The network paradigm has proven to be a successful framework to study the intricate patterns of relations among the constituents of real-world complex systems, from the Internet to the human brain [1,2], and has revealed that the dynamical behaviours observed in such systems, such as information spreading, diffusion, opinion formation and synchronisation, are quite often affected -and to some extent determined-by the structure of the underlying interaction network [3][4][5][6].However, the recent availability of massive data sets of social, technological and biological systems has suggested that the classical complex network approach might fall somehow short in modelling systems whose elementary units can interact through more than one type of connections.This is typical of many real-world systems, such as social networks, where people are connected through a variety of social relationships, e.g.kinship, friendship, collaboration, competition, or transportation systems, which often exploit different communication channels [7][8][9][10].Such systems can be treated in terms of multi-layer or multiplex networks [11][12][13][14], where each layer describes a particular type of interaction among the nodes of the graph.Some recent works have confirmed that multi-layer networks are characterised by new levels of complexity [15] and that the interaction of multiple network layers can produce new interesting dynamical behaviours [16][17][18][19][20]22].
In the realm of dynamical processes on networks [23] the simplicity and -still-the richness of random walks has attracted much attention in recent years [24,25].Ran-dom walks are the most simple way to explore a network using only local information, and the steady-state properties of a walk, including characteristic times, limiting occupation probability, and coverage, have tight relationships with the structure of the graph upon which the walk takes place [26,27].For this reason, random walks have also been successfully used as probes of network properties, with applications ranging from community detection [28][29][30] to taxonomy of real-world networks [31].Moreover, specific flavours of random walks are widely used for the exploration of online social networks, information networks and the like.
A class of random walkers of particular interest is that of walkers whose motion is biased on the structural properties of the network [32].In its simplest possible version, the considered biased random walks are Markov processes whose transition probability is a parametric function of the topological properties of the destination node.In this way, by tuning the parameters of the biasing function one can force the walk to preferentially visit, or avoid, nodes exhibiting high or low values of given topological descriptors, such as the degree, clustering or betweenness.In particular, degree-biased random walks have been used to define new centrality measures [33,34], identify communities [35], and provide optimal exploration of a network using only local information [36].It has also been found that the dynamics of degree-biased random walks is strongly affected by the presence of degree-degree correlations in the structure of the network [37][38][39], so that an appropriate choice of the structural bias can be used to perform efficient sampling of unknown networks.
In this Article we study several ways in which random walks can be extended to multi-layer networks, and we show how to devise appropriate ways to bias the walkers on the topological properties of the nodes at each layer in order to perform an efficient exploration of such systems.We notice that random walks have already been applied to multi-layer networks, e.g. to quantify the impact of failures in interconnected systems [40].However, we will focus here on biased random walks and will investigate how the biasing function affects the dispersiveness of the walk and the steady-state occupation probability distribution.The aim is to find walks which visit far away regions of a multiplex network within a relatively small number of steps, a property related to the dispersiveness of the walk, and, at the same time, guarantee that the probability for a walker to visit any node in the system is as close as possible to uniform, thus allowing to sample unknown graphs in an efficient way.
The presence of many interdependent layers allows to construct several classes of biased random walks, and in particular what we call extensive walks and intensive walks, where the difference between the two classes is in the dependence of the parameters of the biasing function on the number of layers of the system.In the former case, the biasing function depends on the structural properties of the destination node at all the layers of the system (thus, the number of parameters is extensive in the number of layers), while in the latter case the bias depends on intrinsically multiplex properties of the destination node, which do not depend explicitly on the number of layers of the network.
For both classes of biasing functions, we provide analytical closed forms for the long-time properties of the random walks, in terms of stationary probability distribution and entropy rate [41], and we study the effect of different structural properties, including the number of layers, the presence and sign of inter-layer degree correlations, the redundancy of edges across layers, the density of the multiplex and the heterogeneity of the degree distributions, on the steady-state behaviour of these walks.We find that all these properties have a remarkable effect on the maximal dispersiveness and on the steady-state occupation probability of biased random walks.
Finally, we study the diffusion properties of several real-world multiplex networks, namely the six continental airline transportation networks, and we show that in those cases the pressure to provide robust route alternatives has somehow hindered the overall diffusion properties of those systems.

GENERAL FEATURES OF BIASED RANDOM WALKS
Let us consider a M -layer multiplex network of N nodes, i.e. a multi-layer graph in which each node can interact with the other ones by means of M different kinds of relationships.A multiplex is fully described by the Mdimensional array of the adjacency matrices of its layers A = {A [1] , A [2] , . . ., A [M ] }, where and a [α] ij = 1 if node i and node j are connected at layer α.
In the following we assume the layers to be unweighted, but all the results can be easily extended to to the case of weighted multiplexes.
In general, a random walker on a multiplex is not constrained on a single layer and can exploit all the connections pointing out of the current node, at all layers.A synthetic -yet incomplete-description of the topology of a multiplex is provided by the overlapping adjacency ma- ij account for the total number of connections between two nodes across all layers [12].In particular, we consider the class of Markovian random walks defined by the transition probabilities: This set up is very general and allows for a variety of different motion rules.In fact, f j can be either a function of some topological multiplex properties of the arrival node j, or an informative combination of some structural features of the destination node, measured at all or at a fraction of the layers.Notice that the unbiased random walk on the multiplex is obtained by setting f j = 1, ∀j ∈ V .
In this case a walker jumps out of node i by traversing one of the edges incident on i chosen with uniform probability and independently on the layer to which it belongs.It is worth noting that the use of the overlapping adjacency matrix {o ij } does not automatically make the walk in Eq. ( 1) equivalent to a random walk on the aggregated graph obtained by flattening all the layers in a single network.In general, if the biasing function f j depends, either explicitly or implicitly, on the structural properties of node j in the multiplex network, the walk in Eq. ( 1) cannot be directly mapped on an equivalent walk on the aggregated graph.Stationary probability distribution.-Starting from the one-step transition probability given in Eq. 1 we derive closed forms for several asymptotic properties of the walk.Following an approach similar to that used in Ref. [32], we now show that for any choice of the biasing function f j the stationary probability distribution p * = {p * i } of biased walks on multiplex networks can be analytically derived, under the hypotheses that i) the topological overlapping matrix O is primitive and that ii) f j is a time-invariant function of any property of the destination node j.We start by considering the probability p i→j (t) that a walker starting at node i will be found on node j after exactly t time steps: and the dual probability p j→i (t): Comparing Eq. ( 2) with Eq. ( 3) and considering that the multiplex is undirected (i.e., o ij = o ji ), we obtain where c i = j o ij f j .If the matrix O is primitive, then a stationary probability distribution exists and lim t→∞ p i→j (t) = p * j , leading to the expression: By imposing the normalisation condition j p * j = 1 we finally get: We notice that Eq. ( 6) is quite general, since it does not explicitly depend on the form of the biasing function or on the actual structure of each layer or of the topological overlapping matrix O.
In many real-world application scenarios, e.g. in crawling the structure of online social networks, it is important to guarantee that for long enough times the walk will end up visiting all the nodes of the graph with the same probability.It is easy to prove that an unbiased random walk is not a good choice in this case, since its steady-state occupation probability distribution is proportional to the degree sequence, hence an appropriate bias should be used to avoid to visit hubs more frequently than poorly-connected nodes.In practice, it is not always possible to find a walk which produces exactly the same stationary occupation probability distribution for all the nodes, i.e. p * i = p = 1/N, ∀i.However, one could instead require that the resulting stationary probability distribution, although not equal for all nodes, has the minimum possible variance.In particular, in the following we will focus on the normalised standard deviation of the stationary probability distribution: where µ(p * ) and σ(p * ) are the average and the standard deviation of p * , respectively.We will look for suitable combinations of the parameters of the walk that produce the smallest possible value of η(p * ), corresponding to the maximum uniformity of the accessibility of the nodes attainable on a certain multiplex network.Entropy rate.-One classical measure to quantify the mixedness or dispersiveness of a walk on a graph is the entropy rate h = lim t→∞ S t /t [41], where S t is the Shannon entropy of the set of all the trajectories of length t generated from the walk rule, and h is the minimum amount of information necessary to describe the process [41].In particular, h = 0 only if the walk generates exactly one possible trajectory, while h is maximum when all the trajectories are equiprobable.Intuitively, walks with a high mixedness can explore remote regions of a graph within a relatively small number of time-steps.This property is again desirable for the efficient exploration of unknown networks, where only local information is available.In particular, it is interesting to find a biasing function which guarantees that the walk does not remain trapped for too long in any region of the graph, and this is usually obtained by maximising the dispersiveness of the walk.
It is possible to show that the entropy rate for a Markov process can be expressed as which means that h depends only on the walk rule π ij and on the stationary probability distribution [32].By substituting the analytical expression for p * given in Eq. ( 6) into Eq.( 8) we get: (9) This expression has a natural upper bound, which reflects the case of random walks where all trajectories of the same length have equal probability.It is interesting to notice that, as shown by Burda et al. in Ref. [42], the maximal value of entropy rate attainable by any walk on a given single-layer graph depends on the structure of the graph, and in particular for an undirected graph it is equal to ln λ max , where λ max is the maximum eigenvalue of the adjacency matrix of the graph.
This result can be extended to the case of walks on multiplex networks as follows.The total number of trajectories of length t generated by a walk defined as in Eq. ( 1) is equal to N t = i,j (O t ) ij , where O t is the t-th power of the overlapping adjacency matrix.In the limit of large t, we have where λ max is now the maximum eigenvalue of the overlapping adjacency matrix O (this result is a direct consequence of the application of the power method).In general, the maximal value of the entropy rate attainable with a particular motion rule will be smaller than or at most equal to h max .Since obtaining high mixedness is a desirable property of a walk in many real-world applications, such as when searching for a given resource on a graph, in the following we will look for combinations of the parameters of different motion rules which can produce high values of h, to better approximate the corresponding value of h max allowed by the structure of the network.
Heterogeneous mean-field.-In the particular case in which the bias function f i depends only on the (vectorial) degree ij is the degree of node i at layer α, the expression for the stationary probability distribution can be considerably simplified.Let us consider a heterogeneous mean-field, in which all the nodes belonging to the same degree class k are structurally indistinguishable.Under these assumption, and since f i depends only on the degree, then for all the nodes i having the same degree and similarly: where C is an appropriate normalisation constant to ensure that k p * k = 1.Eq. ( 11) means that all the nodes in the same degree class will have the same steady-state probability of being visited by the walk.Notice that o kk is the expected number of edges connecting two nodes whose multiplex degree is respectively equal to k and to k .If we assume that there are no edge correlations, i.e. that the probability of having a since the expected number o kk of edges between a node with degree k and a node with degree k is actually equal to the sum of the expected number of edges connecting these two nodes at each of the M layers (we indicate by k [α] the degree at layer α of a node whose vectorial degree is equal to k ).If we additionally assume that there are no intra-layer correlations, then: where P (k [α] ) is the degree distribution at layer α.In the end we find: This expression for p * k is quite general, and in particular it is valid even in the presence of inter-layer degreecorrelations [43].Since the heterogeneous mean-field discards intra-layer and edge correlations, which usually contribute to hinder the dispersiveness of a walk, Eq. ( 14) can be readily plugged into the expression of the entropy rate in Eq. ( 8) to obtain an estimate of the maximum value of h attainable with a given biasing function on a multiplex network with an assigned multiplex degree sequence {k i }.

CLASSES OF BIASED RANDOM WALKS
The introduction of a biasing function in the motion rule is mainly motivated by the necessity to obtain an exploration of the graph which is more efficient, i.e., faster with respect to the time needed to visit all the nodes, or more homogeneous, i.e., avoiding heterogeneities in the stationary distribution probability, in order to explore with the same probability each node of the graph.In single layer networks these two aims are in general antithetical.For instance, a biasing function which maximises the mixing of the walk (corresponding to higher values of entropy rate) usually produces a quite heterogeneous stationary occupation probability, mainly due to the fact that a better mixing is obtained by exploiting the central role played by hubs.High values of h are usually achieved in a single-layer uncorrelated graph by a degreebiased walk π ji ∼ k b j with b = 1, and in general with a bias b > 0 in graphs with non-trivial degree-degree correlations [32].On the other hand, a uniform stationary occupation probability is obtained by using π ji ∼ k b j with b = −1 in uncorrelated graphs, and in general by a value of b < 0 for graphs with degree-degree correlations, which corresponds to forcing the walkers to preferentially move towards poorly connected nodes [38].
The richness of multi-layer networks allows the exploration of more complex biasing functions and, as we will show in the following, usually produces quite interesting dynamics.The reason of such richness is that the multiplex degree of a node i is a vectorial rather than a scalar quantity, a fact that allows to construct several degreebased biasing functions.In the following we present two particular classes of such biasing functions, which we call extensive and intensive biases, respectively.
Extensive bias functions.-We call extensive those walks whose motion rule depends on a function of the degrees of the destination node at each of the M layers.A first example is that of additive degree-biased walks, defined by transition probabilities of the form: where b α ∈ R is the bias exponent associated to layer α.
Another example is that of multiplicative degree-biased walks, whose transition probability takes the form: We named these walks "extensive" since the number of free parameters in the motion rule, namely the exponents b α , increases with the number of layers M .This peculiar property of extensive walks allows for a fine-grained setting of the bias in order to avoid nodes whose replicas on each of the M layers belong to a specific degree class.For instance, in the case of a two-layer multiplex, if we set b 1 > 0 and b 2 < 0 then the walkers will preferentially move towards node having, at the same time, high degree on layer 1 and low degree on layer 2. It might sometimes be desirable for a walker to have such sophisticated motion rules.An example is that of multiplex collaboration networks, in which nodes are scientists and layers represent co-authorship patterns in different fields.In that case, we might use an appropriately biased multiplex random walk which prefers to move towards nodes having a higher degree in a particular field, whose stationary probability distribution will represent a measure of the relative importance of each author in that field.
However, having a number of parameters which scales with the number of layers is not always a desirable property, especially if one wants to tune these parameters in order to obtain a walk with certain dynamical properties (e.g., either in terms of stationary probability or in terms of entropy rate).This problem is efficiently solved by intensive bias functions.
Intensive bias functions.-We call intensive those multiplex walks whose motion bias depends on one or more intrinsically multiplex properties of the destination node.In the following we will focus on the intensive walk whose transition probability reads: where is the overlapping degree of node j and P j is the multiplex participation coefficient of j, and is defined as [12]: We notice that by considering o • and P • we are effectively using information about the distribution of the edges of the destination node across the layers.In particular, for fixed number of layers M , o i is proportional to the average of the distribution defined by i , . . ., k [M ] i }, while P i gives information about the homogeneity of k i , with i ∀α (i.e., if node i has roughly the same degree at all layers) and P i ∼ 0 if almost all the edges of node i lie on just one layer.
We notice that when b o > 0 the walkers will preferentially move towards hubs, while for b o < 0 they tend to visit the poorly connected nodes more often.Similarly, for positive values of b p the walkers will preferentially move towards truly multiplex nodes, i.e. nodes whose distribution of edges across the M layers is more homogeneous, while for b p < 0 the walkers prefer to move towards focused nodes, i.e. those having the majority of their connections in just one or a few of the M layers [12].In general, by tuning the two parameters b o and b p we can obtain a rich variety of different walks.For instance, for b o > 0 and b p > 0, the walkers will be attracted by truly multiplex hubs (i.e., nodes with many links, almost equally distributed across the layers).Conversely, when b o > 0 and b p < 0 focused hubs are visited often and multiplex poorly connected nodes are strongly avoided, and so forth.The unbiased multiplex walk is recovered for b p = b o = 0.
The most interesting characteristic of the intensive walk defined by Eq. ( 17) is that the number of free parameters is fixed and does not scale with the number of layers, as instead happens for extensive walks.We will show in the following that intensive walks usually perform at least as well as extensive walks, e.g. with respect to the maximisation of entropy rate or to the minimisation of heterogeneity in the stationary occupation probability distribution.
It is worth noting that in the case of a duplex, i.e. when M = 2, even if the number of biasing parameters in intensive and extensive walks is the same, their effect on the motion of the walkers is different.Differently from b 1 and b 2 , intensive biases do not allow to bias the walkers towards nodes with given properties in a particular layer but always consider intrinsically multiplex features, such as their total number of connections and their heterogeneity.
In order to explore the differences in the dynamical properties (i.e., the entropy rate h and the normalised standard deviation of the stationary occupation probability distribution η(p * )) of biased multiplex walks, in the top panel of Fig. 1 we report the values of h obtained by additive, multiplicative and intensive random walks as a function of the two bias exponents in a two-layer multiplex network whose layers have the same average degree k and power-law degree distributions P (k) ∼ k −γ with γ = 2.5, with no inter-layer correlations and no edge overlap [45].We notice that also in this simple case the three walks have remarkably different behaviours.In particular, the additive walk exhibits a relative small sensitivity to the values of the biasing exponents, which results in smaller variations of h.In fact, there is a large region of b 1 (i.e.0 < b 1 < 2) within which the entropy rate is almost constant and not very different from the Conversely, the picture is much richer and less trivial in the case of multiplicative and intensive walks, for which the maximum of h is obtained for a relatively small range of parameters, usually corresponding to positive exponents.We obtain slightly different results when we consider two layers with different power-law degree distributions P (k [1] ) ∼ k [1] γ1 and P (k [2] ) ∼ k [2] γ2 , namely with exponents γ 1 = 2.2 and γ 2 = 2.7 respectively.In this case, the symmetry in the additive and multiplicative phase diagrams is broken, and the maximum values of h are found by biasing the walk towards nodes with high degree on both layers, with a higher biasing exponent on the degree of the second layer, which has a more homogeneous degree distribution.Also the phase diagram for the intensive walk is modified, with the line of maximum values becoming thinner.Similar considerations hold for the phase diagram of η(p * ), reported in Fig. 2. In this case, the minimum variance (yielding a more homogeneous exploration of nodes) is obtained for negative values of the two bias exponents.Moreover, the phase diagram exhibits quite small varia-tions in the case of additive walk, while we observe more heterogeneity in the case of multiplicative and intensive walks.Again, the symmetry of the phase diagrams of the extensive walks is broken when pairs of layers with different power-law exponents γ 1 , γ 2 are considered, with the region b 2 > b 1 showing greater variations than for b 2 < b 1 .Qualitatively similar differences can be obtained with asymmetric layers with respect to other statistical properties, such as density.
All the results for synthetic networks, both in the current and following sections, have been obtained for layers with N = 10 4 nodes and averaged over 1000 realisations.

HOW THE STRUCTURE OF A MULTIPLEX AFFECTS THE WALK
In this section we illustrate how the structure of the multiplex network affects the maximal entropy rate and the minimum heterogeneity of the stationary occupation probability distribution achievable in the system.
We focus on five structural aspects, namely i) the presence and sign of inter-layer degree-degree correlations, ii) the existence of edge overlap across layers, iii) the number M of layers of the multiplex, iv) the power-law exponent γ of the degree distribution of the layers, and v) FIG. 2: Heat-maps of the normalised standard deviation of the stationary occupation probability distribution η(p * ) of different multiplex biased random walks.Legend as in Fig. 1.In extensive walks, the minimum of η is always attained for negative values of the two exponents, while in intensive walks the minimum of η is obtained for bo < 0 and bP 0, meaning that walkers tend to preferentially move towards nodes with small degrees on both layers.their density, measured through the average degree k .Since our focus is on the construction of efficient walks (in terms of maximal dispersiveness and of homogeneity of the stationary occupation probability) the parameters of interest in all the cases are the overall maximum value of entropy rate, denoted by h max , and the minimum value of the normalised standard deviation, denoted by η min , obtained by extensive and intensive biased random walks as a function of the biasing parameters.
Effect of inter-layer degree correlations.-In a recent work [43] the authors have shown that real-world multiplex networks are usually characterised by non-trivial inter-layer degree correlation patterns.In the same paper the authors propose several methods to quantify the presence of inter-layer correlations between a pair of layers, including the rank correlation among the two degree sequences, as measured by the Spearman's coefficient ρ.

If we call R
[α] i the rank of node i due to its degree on layer α, the Spearman rank correlation coefficient between layer α and layer β reads: where R [α] and R [β] are the average ranks of nodes respectively at layer α and layer β.The coefficient ρ takes values in [−1, 1], so that ρ = 1 if the two degree sequences are perfectly correlated (meaning that a hub at layer α is also a hub at layer β), while ρ = −1 when the two degree sequence are perfectly anti-correlated, i.e. when a hub on layer α is always a poorly connected node on the other layer, and viceversa.Intermediate positive (resp.negative) values of ρ indicate weaker positive (negative) inter-layer correlations, while ρ 0 when the two degree sequences are uncorrelated.
In Fig. 3(a) we report the plot of h max and η min for extensive and intensive walks on two-layer multiplex networks with same average degree and power-law degree distributions P (k) ∼ k −γ with γ = 2.5, for different levels of inter-layer degree correlations.As made evident by the figure, intensive walks usually perform at least as well as extensive walks with respect to both maximisation of entropy and minimisation of the heterogeneity of the stationary occupation probability distribution.This suggests that, aside from the actual differences in the phase space, intensive walks are able to span the same range of values of entropy and η(p * ) by using only two parameters, irrespective of the actual numbers of layers in the multiplex.
Effect of edge overlap.-We now investigate the impact of the presence of edge overlap on the long-term dynamics of extensive and intensive walks.We recall here the definition of overlap for an edge (i, j), which is the fraction of layers in which the edge (i, j) exists [12,44], i.e.: The edge overlap of a multi-layer network is defined as the average of ω ij over all the node pairs for which o ij = 0 (i.e., for all pairs of nodes which are connected by at least one edge): where K is the number of pairs of nodes which are connected in at least one of the M layers.Notice that the average edge overlap ω is equal to 1 only if all the M layers are identical, while ω = 1/M when every edge is present on exactly one of the M layers.We started from two-layer multiplex networks obtained by coupling identical layers (thus having edge overlap equal to 1) with power-law degree distributions P (k) ∼ k −γ with γ = 2.5, and then we obtained multiplex networks with prescribed values of edge overlap by rewiring a certain percentage of the edges of one of the two layers in order to maintain the degree sequence unaltered.Notice that by construction the resulting multiplex networks have maximally positive inter-layer degree correlations (i.e., ρ = 1).As shown in Fig. 3 as ω increases, meaning that higher values of edge overlap correspond to a more heterogeneous stationary state probability distribution.Conversely, h max decreases with ω, in accordance with the fact that higher edge overlap tends to hinder the dispersiveness of the walk, since a smaller number of distinct walks can originate from each node.Summing up, multiplex networks having smaller values of edge overlap are overall preferable in order to maximise the dispersiveness of the walk and to obtain a more homogeneous stationary occupation probability.
In other words, a small edge overlap guarantees a more effective exploration of a multiplex network and, at the same time, a more homogeneous distribution of the probability of visiting each node.
Effect of the number of layers.-It is also interesting to study how the dynamical properties of intensive walks change when the number of layers M is progressively increased.To this aim, we constructed multiplex networks with different number of layers, with no inter-layer degree correlations and negligible edge overlap, where all the layers had power-law degree distributions P (k) ∼ k −γ with γ = 2.5.As shown in Fig. 3(c), h max is an increasing function of M , while η min decreases as the number of layers grows.In general, the addition of layers in absence of inter-layer correlation flattens the structural differences among the nodes of the multiplex, and provides better dispersiveness and less heterogeneity in the occupation probability distribution.
Effect of the heterogeneity of the degree distribution.
-We investigate here how the heterogeneity of the the degree distribution of each layer affects h max and η min .
To this aim, we considered pairs of uncorrelated layers with the same power-law degree distribution P (k) ∼ k −γ for different values of γ, maintaining fixed the average degree of the networks k .The plots in Fig. 4(a) confirm that both h max and η min grow as γ increases, i.e. as the degree distribution of the layers becomes more homogeneous.We notice though that the variation in η min appears to be relatively smaller, especially for multiplicative and intensive walks.This result can be explained by considering that dispersiveness is favoured by more homogeneous degree distributions.Layers with different power-law exponents γ 1 and γ 2 have been considered in the previous section.
Effect of layer density.-Finally, we focus on the effect of layer density, measured through the average degree of the layers k .Once again we report here the case of uncorrelated layers with power-law exponent γ = 2.5, but similar results have been obtained for other values of γ.As shown in Fig. 4(b), both h max and η min increase as a function of k .Layers with different average degrees k [1] and k [2] break the symmetry of the phase diagrams for h and η qualitatively in a similar way as pairing layers with different power-law exponents.
Summing up, the analysis of the impact of structural properties on the values of h max and η min attainable on a multiplex network confirms that positive inter-layer degree correlations, small edge overlap, large number of layers, and more homogeneous layers all concur towards allowing biased random walks with nearly-optimal dispersiveness and closely-to-homogeneous steady-state visiting probability.In other words, a multiplex network with a large number of layers and small edge overlap, where nodes have roughly the same number of links at all layers, can be explored ways more efficiently than a similar multiplex network where nodes have disassortative inter-layer correlations and edges are redundant across layers.
In the following section we show that the multiplex airline transportation networks of all the six continents have evolved towards a structure which provides a good trade-off between efficient exploration and robustness.

APPLICATIONS TO REAL-WORLD AIRLINE TRANSPORTATION NETWORKS
As an application, we study here the dynamical properties of multiplex biased walks on a set of real-world systems, namely the six continental airline transportation networks.In such systems nodes represent airports, edges indicate the existence of a route between two airports and each layer is associated to an airline company, i.e. all the edges in a layer represent the routes operated by the corresponding airline.Such networks have been introduced and extensively studied in Ref. [43].As shown in Table I, all such multiplex networks consist of a relatively high number of layers.For this reason, we will use intensive walks to compute the maximal entropy rate h max and the minimum value of the standard deviation of the stationary distribution η min .In Table I, we also report for each multiplex the average number of layers M × ω where each edge exists, the theoretical upper value of entropy rate ln λ, and the values of h max and η min obtained by optimising intensive walks.We notice that the efficiency of a transportation system is usually measured in terms of the accessibility of the locations it serves.In particular, in an ideal transportation system it should be easy to travel between any pair of far-apart regions of the network, mostly irrespective of where exactly those locations are located.Now, discarding the cost associated to the distance between the nodes of an airline transportation network, high accessibility can be obtained by guaranteeing that a traveller can reach remote locations in the system without large effort, in terms of number of interchanges, and that all locations can be visited with comparable effort.We have seen that in the language of random walks these two criteria correspond, respectively, to the maximisation of dispersiveness and to the minimisation of the standard deviation of the visiting probability.
Hence, we can ask whether the six continental air transportation systems can guarantee a good level of navigability, i.e. an optimal trade-off between dispersiveness and homogeneity of the visiting probability.We reckon that a more informative analysis of the efficiency of these systems would require more detailed information about the actual patterns of trips travelled by passengers, the cost associated to each route, the presence of non-Markovian effects (people often come back to their original place at the end of a trip), the non-stationarity of the system due to seasonality, etc..However, we argue that biased random walks can still provide useful, yet coarse-grained, information about the overall navigability of those systems.Since we cannot modify the degree distributions of each of the layers, or the patterns of inter-layer correlations, or the actual number of layers in each continental air transportation system, we focus here in particular on the effects of edge overlap.
In the previous section we showed that networks with high edge overlap ω achieve lower maximal values of dispersiveness of the walk and larger heterogeneity of the equilibrium occupation probability distribution.When two nodes are connected by more than one edge, indeed, from a dynamical point of the view some connections are wasted, since redundant links do not allow for new paths in the network.However, their redundancy might often be important for a transportation system, since it makes specific connections more robust to single link failures.It is not unrealistic to assume that multi-layer transportation systems from the real-world have developed by satisfying a trade-off between the necessity to provide, at the same time, high diffusivity together with reasonable levels of robustness.
Because of the large heterogeneity in the size and number of layers of the six continental transportation systems, it is necessary to introduce some kind of normalisation which allow to compare the results observed in different systems.In order to test the effect of edge overlap on the diffusion properties of real-world systems, for each of the six multiplex networks we computed the zscore of the average edge overlap: where ω and σ(ω) represent respectively the average value and the standard deviation of the overlap com-puted on an ensemble of suitably randomised multiplex networks.In particular, for each continental airline system we sampled 1000 multi-layer graphs from the configuration model which maintains fixed the degree sequence of all the layers and rearranges the links on each layer, pairing edge stubs at random.We computed also z(h max ) and z(η min ), i.e. the z-scores of the maximal entropy rate and minimum variance over those 1000 multiplex graphs.The results reported in Fig. 5 confirm that also in realworld systems h max is negatively correlated with edge overlap, in agreement with the results obtained on synthetic networks.Similarly, η max is positively correlated with ω.Notice that we have z(ω) > 0 in all the six continents, meaning that the edge overlap of the real-world systems is always higher than that of the null-model, in agreement with the observation that real-world transportation networks tend to guarantee a certain level of robustness to failures.However, the quest for robustness has a cost in terms of dispersiveness and accessibility.In fact, h max is consistently smaller than the value observed in the randomised systems (z(h max ) < 0) for all continents, and similarly the steady-state probability distribution is consistently larger than that observed in the null model (z(η min ) > 0) It is quite interesting to note that the two multiplex networks with smallest overlap and overall better diffusion properties are the continental networks of Oceania and Europe, which span the least geographical space.We can speculate that in such systems some nodes representing cities in different countries are connected comparably well by different modes of transport, such as trains and bus, suited for relatively short distances and not included in our analysis.This might potentially explain the relative low number of redundant edges in those two airline transportation systems.Conversely, the necessity to provide route redundancy has somehow forced the airline transportation networks of Asia, South America and North America towards slightly less efficient configurations.

CONCLUSIONS
In our work we have explored how to extend biased random walks to the case of multiplex networks, showing that the richness of multi-layer systems allows to define several different classes of walks.In particular we studied the general features of the so-called extensive walkers (where the node properties, as the degree, are considered separately at each of the layers with different biasing parameters) and intensive walkers (biased on of the product two intrinsically multiplex, namely the overlapping degree and the participation coefficient) finding closed forms for the stationary occupation probability of these walks and for the entropy rate, and provided simplified heterogeneous mean-field expressions for the case in which the multiplex has no correlations.
We thoroughly investigated how structural properties of the multiplex, such as its number of layers, the presence of edge overlap and/or inter-layer degree correlations, the density of the layers and the heterogeneity of their degree distribution affect the dynamics of the random walkers.We found that number of layers, edge overlap and inter-layer degree correlations have a substantial impact on the diffusion properties of the walks.Also, we found that intensive random walkers perform at least as well as extensive random walkers in all the considered scenarios, with the advantage that the number of bias parameters does not scale with the number of layers.
Finally, the study of the diffusion properties of six real-world multiplex networks, namely the continental airline transportation networks of Africa, Asia, Europe, Oceania, North and South America, has shed some new light on the interplay between efficiency and robustness in multi-layer transportation systems.In particular, we found that the emerging necessity to provide some resilience to single link failures, which corresponds to the introduction of some level of edge overlap, has shaped these systems in such a way that their navigability, in terms of entropy rate and heterogeneity of the node occupation probability, has somehow been sacrificed in favour of robustness.The results of the present work represent a valuable theoretical contribution to the development of efficient strategies to explore, search or navigate multiplex networks, and confirm the importance of appropriately taking into account the multiplexity of interactions when modelling intrinsically multi-dimensional systems.
not depend on the probability of having a[β] ij = 1 for all the possible β = α, then we can write:

FIG. 1 :
FIG.1: Heat-maps of the value of entropy rate h of different multiplex biased walks as a function of the parameters of the biasing function.The panels correspond, respectively, to additive [right, (a) and (d)], multiplicative [middle, (b) and (e)] and intensive walks [left, (c) and (f)] on uncorrelated duplex networks (in the top panels the two layers have the same power-law degree distribution P (k) ∼ k −γ with γ = 2.5, while in the bottom panels the two layers have power-law degree distributions with different exponents, namely γ1 = 2.2 and γ2 = 2.7.In general, the maximum of h is obtained for positive values of the two bias parameters, corresponding to extensive walks which move preferentially towards nodes having high degrees on both layers, and to intensive walks whose motion rule is biased towards truly multiplex nodes.

FIG. 3 :
FIG.3: Values of hmax (top panels) and ηmin (bottom panels) as a function of the the inter-layer degree correlation coefficient ρ (a), the average edge overlap ω (b) and the number of layers M (c), respectively for additive (triangles), multiplicative (squares) and intensive (circles) walks.For the entropy rate, we also show the value of hmax = ln λ corresponding to the maximum entropy random walk (solid line).(a) For all walks, hmax is an increasing function of the inter-layer degree correlation coefficient ρ, and provides a very good approximation of the maximum theoretical entropy rate hmax.Notice that intensive walks perform at least as well as the extensive ones.(b) As the overlap increases, the estimates of hmax obtained by the biased walks become less precise, while ηmin increases as a function of ω.(c) hmax increases and ηmin decreases as a function of M .In this case we only performed simulations for intensive walks.

FIG. 4 :
FIG.4: Values of hmax (top panels) and ηmin (bottom panels) as a function of the exponent γ of the the power-law distribution of each layer (a) and of the average degree k (b), respectively for additive (triangles), multiplicative (squares) and intensive (circles) walks.For the entropy rate, we also show the value of hmax = ln λ corresponding to the maximum entropy random walk (solid line).As shown, for all walks hmax appears to increase as a function of both γ and k .Smaller variations are also found in the values of ηmin.

FIG. 5 :
FIG.5: z-score of the average edge overlap ω versus the z-scores of the maximal entropy rate hmax (a) and the minimum standard deviation of the stationary distribution ηmin (b) obtained through intensive walks.In agreement with findings for synthetic networks, also in real-world systems z(hmax) is negatively correlated with z(ω) -Pearson's correlation coefficient r = −0.70-, whereas z(ηmin) is positively correlated -r = 0.30, which increases to r = 0.67 excluding the outliner North America -.

TABLE I :
Structural properties of the six continental airline transportation systems.For each multiplex, we report the number of nodes N , the number of layers M , the average number of layers in which an edge exists M ×ω, the theoretical upper value of entropy rate ln λmax and the extremal values hmax and ηmin obtained by optimising intensive walks.