Efficient exploration of multiplex networks

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 M-dimensional array of the adjacency matrices of its layers ${ \mathcal A }=\{{A}^{[1]},{A}^{[2]},\ldots ,{A}^{[M]}\}$, where ${A}^{[\alpha ]}=\{{a}_{{ij}}^{[\alpha ]}\}\in {{\mathbb{R}}}^{N\times N}$ and ${a}_{{ij}}^{[\alpha ]}=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 matrix ${ \mathcal O }={o}_{{ij}}$, whose entries ${o}_{{ij}}={\sum }_{\alpha }{a}_{{ij}}^{[\alpha ]}$ 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:

$$\begin{eqnarray}&&{\pi }_{{ji}}=\displaystyle \frac{{o}_{{ij}}{f}_{j}}{{\displaystyle \sum }_{j}\;{o}_{{ij}}{f}_{j}}.\end{eqnarray} \tag{ 1 }$$

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,\;\forall j\in 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 equation (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 equation (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 equation (1) we derive closed forms for several asymptotic properties of the walk. Following an approach similar to that used in [32], we now show that for any choice of the biasing function f_j the stationary probability distribution ${{\boldsymbol{p}}}^{*}=\{{p}_{i}^{*}\}$ of biased walks on multiplex networks can be analytically derived, under the hypotheses that (i) the topological overlapping matrix ${ \mathcal 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\to j}(t)$ that a walker starting at node i will be found on node j after exactly t time steps:

$$\begin{eqnarray}&&{p}_{i\to j}(t)=\displaystyle \sum _{{j}_{1},{j}_{2},\ldots ,{j}_{t-1}}\;{\pi }_{{j}_{1},i}\times {\pi }_{{j}_{2},{j}_{1}}\times \;\ldots \;\times {\pi }_{j,{j}_{t-1}},\end{eqnarray} \tag{ 2 }$$

and the dual probability ${p}_{j\to i}(t)$:

$$\begin{eqnarray}&&{p}_{j\to i}(t)=\displaystyle \sum _{{j}_{1},{j}_{2},\ldots ,{j}_{t-1}}\;{\pi }_{{j}_{1},j}\times {\pi }_{{j}_{2},{j}_{1}}\times \ldots \times {\pi }_{i,{j}_{t-1}}.\end{eqnarray} \tag{ 3 }$$

Comparing equation (2) with equation (3) and considering that the multiplex is undirected (i.e., ${o}_{{ij}}={o}_{{ji}}$), we obtain

$$\begin{eqnarray}&&{c}_{i}{f}_{i}{p}_{i\to j}(t)={c}_{j}{f}_{j}{p}_{j\to i}(t),\quad \forall i,j\in V,\end{eqnarray} \tag{ 4 }$$

where ${c}_{i}={\sum }_{j}\;{o}_{{ij}}{f}_{j}$. If the matrix ${ \mathcal O }$ is primitive, then a stationary probability distribution exists and ${\mathrm{lim}}_{t\to \infty }\;{p}_{i\to j}(t)={p}_{j}^{*}$, leading to the expression:

$$\begin{eqnarray}&&{c}_{i}{f}_{i}{p}_{j}^{*}\;=\;{c}_{j}{f}_{j}{p}_{i}^{*}.\end{eqnarray} \tag{ 5 }$$

By imposing the normalisation condition ${\sum }_{j}\;{p}_{j}^{*}\;=\;1$ we finally get:

$$\begin{eqnarray}&&{p}_{i}^{*}\;=\;\displaystyle \frac{{c}_{i}{f}_{i}}{{\displaystyle \sum }_{{\ell }}\;{c}_{{\ell }}{f}_{{\ell }}}.\end{eqnarray} \tag{ 6 }$$

We notice that equation (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 ${ \mathcal 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}^{*}\quad =\quad \bar{p}=1/N,\;\forall 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:

$$\begin{eqnarray}&&\eta ({p}^{*})=\displaystyle \frac{\sigma ({p}^{*})}{\mu ({p}^{*})},\end{eqnarray} \tag{ 7 }$$

where μ(${p}^{*}$) and σ(${p}^{*}$) are the average and the standard deviation of ${{\boldsymbol{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={\mathrm{lim}}_{t\to \infty }\;{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 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 of a Markov process can be expressed as

$$\begin{eqnarray}&&h=-\displaystyle \sum _{i,j}\;{\pi }_{{ji}}{p}_{i}^{*}\mathrm{ln}({\pi }_{{ji}}),\end{eqnarray} \tag{ 8 }$$

which means that h depends only on the walk rule π_ij and on the stationary probability distribution [32]. By substituting the analytical expression for ${{\boldsymbol{p}}}^{*}$ given in equation (6) into equation (8) we get:

$$\begin{eqnarray}&&h=-\displaystyle \frac{1}{{\displaystyle \sum }_{i}{c}_{i}{f}_{i}}\;\left[\displaystyle \sum _{i}\;{f}_{i}\displaystyle \sum _{j}\;{o}_{{ij}}{f}_{j}\mathrm{ln}({o}_{{ij}}{f}_{j})-\displaystyle \sum _{i}\;{f}_{i}{c}_{i}\mathrm{ln}({c}_{i})\right].\end{eqnarray} \tag{ 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 [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 $\mathrm{ln}{\lambda }_{{\rm{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 equation (1) is equal to ${N}_{t}={\sum }_{i,j}\;{({{ \mathcal O }}^{t})}_{{ij}}$, where ${{ \mathcal O }}^{t}$ is the t-th power of the overlapping adjacency matrix. In the limit of large t, we have

$$\begin{eqnarray}&&{\tilde{h}}_{{\rm{max}}}=\underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{\mathrm{ln}{N}_{t}}{t}=\mathrm{ln}{\lambda }_{{\rm{max}}},\end{eqnarray} \tag{ 10 }$$

where λ_max is now the maximum eigenvalue of the overlapping adjacency matrix ${ \mathcal 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 ${\tilde{h}}_{{\rm{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 ${\tilde{h}}_{{\rm{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 ${{\boldsymbol{k}}}_{i}=\{{k}_{i}^{[1]},{k}_{i}^{[2]},\ldots ,{k}_{i}^{[M]}\}$ of node i, where by definition ${k}_{i}^{[\alpha ]}={\sum }_{j}\;{a}_{{ij}}^{[\alpha ]}$ 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 ${\boldsymbol{k}}$ are structurally indistinguishable. Under this assumption, and since f_i depends only on the degree, then for all the nodes i having the same degree ${{\boldsymbol{k}}}_{i}={\boldsymbol{k}}$ we have ${f}_{i}={f}_{{{\boldsymbol{k}}}_{i}}={f}_{{\boldsymbol{k}}}$, but also ${c}_{i}={c}_{{{\boldsymbol{k}}}_{i}}={c}_{{\boldsymbol{k}}}$, and similarly:

$$\begin{eqnarray}&&{p}_{{\boldsymbol{k}}}^{*}=\displaystyle \frac{1}{C}{f}_{{\boldsymbol{k}}}{c}_{{\boldsymbol{k}}}=\displaystyle \frac{1}{C}{f}_{{\boldsymbol{k}}}\displaystyle \sum _{{{\boldsymbol{k}}}^{\prime }}\;{o}_{{\boldsymbol{k}}{{\boldsymbol{k}}}^{\prime }}{f}_{{{\boldsymbol{k}}}^{\prime }}\end{eqnarray} \tag{ 11 }$$

where C is an appropriate normalisation constant to ensure that ${\sum }_{{\boldsymbol{k}}}\;{p}_{{\boldsymbol{k}}}^{*}\quad =\quad 1$. Equation (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}_{{\boldsymbol{k}}{{\boldsymbol{k}}}^{\prime }}$ is the expected number of edges connecting two nodes whose multiplex degree is respectively equal to ${\boldsymbol{k}}$ and to ${{\boldsymbol{k}}}^{\prime }$. If we assume that there are no edge correlations, i.e. that the probability of having ${a}_{{ij}}^{[\alpha ]}=1$ does not depend on the probability of having ${a}_{{ij}}^{[\beta ]}=1$ for all the possible $\beta \ne \alpha$, then we can write:

$$\begin{eqnarray}&&{p}_{{\boldsymbol{k}}}^{*}=\displaystyle \frac{1}{C}{f}_{{\boldsymbol{k}}}\displaystyle \sum _{{{\boldsymbol{k}}}^{\prime }}\;{f}_{{{\boldsymbol{k}}}^{\prime }}\displaystyle \sum _{\alpha =1}^{M}\;{k}^{[\alpha ]}P({{k}^{\prime }}^{[\alpha ]}| {k}^{[\alpha ]})\end{eqnarray} \tag{ 12 }$$

since the expected number ${o}_{{\boldsymbol{k}}{{\boldsymbol{k}}}^{\prime }}$ of edges between a node with degree ${\boldsymbol{k}}$ and a node with degree ${{\boldsymbol{k}}}^{\prime }$ 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}^{\prime }}^{[\alpha ]}$ the degree at layer α of a node whose vectorial degree is equal to ${{\boldsymbol{k}}}^{\prime })$. If we additionally assume that there are no intra-layer correlations, then:

$$\begin{eqnarray}&&P({{k}^{\prime }}^{[\alpha ]}| {k}^{[\alpha ]})={q}_{k{^{\prime} }^{[\alpha ]}}=\displaystyle \frac{{{k}^{\prime }}^{[\alpha ]}P({{k}^{\prime }}^{[\alpha ]})}{\langle {{k}^{\prime }}^{[\alpha ]}\rangle },\end{eqnarray} \tag{ 13 }$$

where $P({{k}^{\prime }}^{[\alpha ]})$ is the degree distribution at layer α. In the end we find:

$$\begin{eqnarray}&&{p}_{{\boldsymbol{k}}}^{*}=\displaystyle \frac{1}{C}{f}_{{\boldsymbol{k}}}\displaystyle \sum _{{{\boldsymbol{k}}}^{\prime }}\;{f}_{{{\boldsymbol{k}}}^{\prime }}\displaystyle \sum _{\alpha =1}^{M}\;\displaystyle \frac{{k}^{[\alpha ]}{{k}^{\prime }}^{[\alpha ]}P({{k}^{\prime }}^{[\alpha ]})}{\langle {{k}^{\prime }}^{[\alpha ]}\rangle }.\end{eqnarray} \tag{ 14 }$$

This expression for ${p}_{{\boldsymbol{k}}}^{*}$ is quite general, and in particular it is valid even in the presence of inter-layer degree-correlations [43]. Since the heterogeneous mean-field discards intra-layer and edge correlations, which usually contribute to hinder the dispersiveness of a walk, equation (14) can be readily plugged into the expression of the entropy rate in equation (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 $\{{{\boldsymbol{k}}}_{i}\}$.

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 degree-biased walk ${\pi }_{{ji}}\sim {k}_{j}^{b}$ 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 ${\pi }_{{ji}}\sim {k}_{j}^{b}$ 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 degree-based 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:

$$\begin{eqnarray}&&{\pi }_{{ji}}\propto \displaystyle \sum _{\alpha =1}^{M}{({k}_{j}^{[\alpha ]})}^{{b}_{\alpha }},\end{eqnarray} \tag{ 15 }$$

where ${b}_{\alpha }\in {\mathbb{R}}$ is the bias exponent associated to layer α. Another example is that of multiplicative degree-biased walks, whose transition probability takes the form:

$$\begin{eqnarray}&&{\pi }_{{ji}}\propto \displaystyle \prod _{\alpha =1}^{M}{({k}_{j}^{[\alpha ]})}^{{b}_{\alpha }}.\end{eqnarray} \tag{ 16 }$$

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₁ > 0 and b₂ < 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 desirable, 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:

$$\begin{eqnarray}&&{\pi }_{{ji}}=\displaystyle \frac{{o}_{{ij}}({o}_{j}^{{b}_{o}}{{ \mathcal P }}_{j}^{{b}_{p}})}{{\displaystyle \sum }_{{\ell }}\;{o}_{i{\ell }}({o}_{{\ell }}^{{b}_{o}}{{ \mathcal P }}_{{\ell }}^{{b}_{p}})},\end{eqnarray} \tag{ 17 }$$

where ${o}_{j}={\sum }_{\alpha }\;{k}_{j}^{[\alpha ]}$ is the overlapping degree of node j and ${{ \mathcal P }}_{j}$ is the multiplex participation coefficient of j, and is defined as [12]:

$$\begin{eqnarray}&&{{ \mathcal P }}_{i}=\displaystyle \frac{M}{M-1}\left[1-\displaystyle \sum _{\alpha =1}^{M}\;{\left(\displaystyle \frac{{k}_{i}^{[\alpha ]}}{{o}_{i}}\right)}^{2}\right].\end{eqnarray} \tag{ 18 }$$

We notice that by considering o_• and ${{ \mathcal P }}_{\bullet }$ 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 ${{\boldsymbol{k}}}_{i}=\{{k}_{i}^{[1]},{k}_{i}^{[2]},\ldots ,{k}_{i}^{[M]}\}$, while ${{ \mathcal P }}_{i}$ gives information about the homogeneity of ${{\boldsymbol{k}}}_{i}$, with ${{ \mathcal P }}_{i}=1$ if ${k}_{i}^{[\alpha ]}=\tfrac{1}{M}{\sum }_{\beta }{k}_{i}^{[\beta ]}\;\forall \alpha$ (i.e., if node i has the same degree at all layers) and ${{ \mathcal 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}\gt 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 equation (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₁ and b₂, 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 $\eta ({p}^{*})$) of biased multiplex walks, in the top panels of figure 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 $\langle k\rangle$ and power-law degree distributions P(k) ∼ k^−γ with γ = 2.5, with no inter-layer correlations and no edge overlap³ . 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₁ (i.e. 0 < b₁ < 2) within which the entropy rate is almost constant and not very different from the absolute maximum for a relatively large range of values of the other exponent b₂, i.e. −5 < b₂ < 2 (the same reasoning is valid for 0 < b₂ < 2 and −5 < b₁ < 2, due to the symmetry of the additive bias function).

Zoom In Zoom Out Reset image size

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]})\sim {({k}^{[1]})}^{{\gamma }_{1}}$ and $P({k}^{[2]})\sim {({k}^{[2]})}^{{\gamma }_{2}}$, namely with exponents γ₁ = 2.2 and γ₂ = 2.7 respectively (bottom panels of figure 1). 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 figure 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 variations 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 γ₁, γ₂ are considered, with the region b₂ > b₁ showing greater variations than for b₂ < b₁. Qualitatively similar differences can be obtained with asymmetric layers with respect to other statistical properties, such as density.

Zoom In Zoom Out Reset image size

All the results for synthetic networks, both in the current and following sections, have been obtained for layers with N = 10⁴ nodes and averaged over 1000 realisations.

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) their density, measured through the average degree $\langle k\rangle$. 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}^{[\alpha ]}$ the rank of node i due to its degree on layer α, the Spearman rank correlation coefficient between layer α and layer β reads:

$$\begin{eqnarray}&&{\rho }_{\alpha ,\beta }=\displaystyle \frac{{\displaystyle \sum }_{i}({R}_{i}^{[\alpha ]}-\bar{{R}^{[\alpha ]}})({R}_{i}^{[\beta ]}-\bar{{R}^{[\beta ]}})}{\sqrt{{\displaystyle \sum }_{i}{({R}_{i}^{[\alpha ]}-\bar{{R}^{[\alpha ]}})}^{2}{\displaystyle \sum }_{j}{({R}_{j}^{[\beta ]}-\bar{{R}^{[\beta ]}})}^{2}}},\end{eqnarray} \tag{ 19 }$$

where $\bar{{R}^{[\alpha ]}}$ and $\bar{{R}^{[\beta ]}}$ 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 sequences are perfectly anti-correlated, i.e. when a hub on layer α is always a poorly connected node on the other layer, and vice versa. Intermediate positive (respectively negative) values of ρ indicate weaker positive (negative) inter-layer correlations, while ρ ≃ 0 when the two degree sequences are uncorrelated.

In figure 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\omega }_{{ij}}=\displaystyle \frac{1}{M}\displaystyle \sum _{\alpha =1}^{M}\;{a}_{{ij}}^{[\alpha ]}.\end{eqnarray} \tag{ 20 }$$

The edge overlap of a multi-layer network is defined as the average of ω_ij over all the node pairs for which ${o}_{{ij}}\ne 0$ (i.e., for all pairs of nodes which are connected by at least one edge):

$$\begin{eqnarray}&&\omega =\displaystyle \frac{1}{{\displaystyle \sum }_{i,j}\;(1-{\delta }_{{o}_{{ij}},0})}\displaystyle \sum _{{ij}}\;{\omega }_{{ij}}=\displaystyle \frac{1}{2K}\displaystyle \sum _{{ij}}\;{\omega }_{{ij}},\end{eqnarray} \tag{ 21 }$$

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 figure 3(b), η_min becomes higher 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 trajectories 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)\sim {k}^{-\gamma }$ with γ = 2.5. As shown in figure 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 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 $\langle k\rangle$. The plots in figure 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 small, 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 γ₁ and γ₂ have been considered in the previous section.

Zoom In Zoom Out Reset image size

Effect of layer density. Finally, we focus on the effect of layer density, measured through the average degree of the layers $\langle k\rangle$. 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 figure 4(b), both h_max and η_min increase as a function of $\langle k\rangle$. Layers with different average degrees $\langle {k}^{[1]}\rangle$ and $\langle {k}^{[2]}\rangle$ 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.

As an application, we study here the dynamical properties of multiplex biased walks on a set of real-world systems, namely the six continental air transportation networks. In such systems nodes represent airports, edges indicate the existence of a direct 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. These networks have been introduced and extensively studied in [43]. As shown in table 1, 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 1, we also report for each multiplex the average number of layers M × ω where each edge exists, the theoretical upper value of entropy rate $\mathrm{ln}{\lambda }_{\mathrm{max}}$, 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 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 allows the comparison of 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 z-score of the average edge overlap:

$$\begin{eqnarray}&&z(\omega )=\displaystyle \frac{\omega -\langle \omega \rangle }{\sigma (\omega )},\end{eqnarray} \tag{ 22 }$$

where $\langle \omega \rangle$ and σ(ω) represent respectively the average value and the standard deviation of the overlap computed 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 figure 5 confirm that also in real-world 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).

Zoom In Zoom Out Reset image size

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 air transportation systems. Conversely, the necessity to provide route redundancy has somehow forced the air transportation networks of Asia, South America and North America towards slightly less efficient configurations.

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 walks (where the node properties, as the degree, are considered separately at each of the layers with different biasing parameters) and intensive walks (biased on a function of two intrinsically multiplex properties, 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 we provided simplified heterogeneous mean-field expressions for the case in which the multiplex has no correlations.

We thoroughly investigated how structural properties of a 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 walks. 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.

FB, VN, and VL acknowledge the support of the EU Project LASAGNE, Contract no. 318132 (STREP). VL also acknowledges support from EPSRC project GALE Grant EP/K020633/1. This research utilized Queen Mary's MidPlus computational facilities, supported by QMUL Research-IT and funded by EPSRC grant EP/K000128/1.