Lévy random walks on multiplex networks

Random walks constitute a fundamental mechanism for many dynamics taking place on complex networks. Besides, as a more realistic description of our society, multiplex networks have been receiving a growing interest, as well as the dynamical processes that occur on top of them. Here, inspired by one specific model of random walks that seems to be ubiquitous across many scientific fields, the Lévy flight, we study a new navigation strategy on top of multiplex networks. Capitalizing on spectral graph and stochastic matrix theories, we derive analytical expressions for the mean first passage time and the average time to reach a node on these networks. Moreover, we also explore the efficiency of Lévy random walks, which we found to be very different as compared to the single layer scenario, accounting for the structure and dynamics inherent to the multiplex network. Finally, by comparing with some other important random walk processes defined on multiplex networks, we find that in some region of the parameters, a Lévy random walk is the most efficient strategy. Our results give us a deeper understanding of Lévy random walks and show the importance of considering the topological structure of multiplex networks when trying to find efficient navigation strategies.


Introduction
2][3][4] Many diverse dynamical processes have been explored on top of networks, including diffusion processes, [5][6][7] synchronization, 8,9 percolation, 10,11 to cite just a few. 123][24] A diversity of random walk processes can be defined and studied, however, most of them rely on the classical random walk process whose dynamics occurs according to the topology of the network. 21n the later scenario, the random walker can only hop to one of the nearest neighbours of the node where it is at any given time, with some -generally the same-probability.Another common random walk process, named Lévy flight, represents the best strategy for randomly searching a target in an unknown environment.This latter kind of random walk dynamics has been widely observed in many animal species. 25,26 n its simplest schematization, this stochastic process could drive a walker over very long distances in a single step event that is called 'flight'. 27The length of the jump, l, obeys a power-law probability distribution in the form of P(l) ∼ l −α , 25 which makes it possible for the random walker to hop from one node to any other node.
On the other hand, multiplex networks, [28][29][30] i.e., networks composed by many different layers of interactions, are gaining much attention recently.The social and technological revolution brought by the Internet and mobile connections, chats, on-line social networks, and a plethora of other human-to-human machine mediated channels of communications have revealed the need to consider that networks might be made up by many different layers of interactions.The same occurs in other fields, like in contemporary biology, where the needs to integrate multiple sets of omic data naturally leads to a multiplex network as a schematization of the system under study.Also in the traditional field of transportation networks, the concept of multiplex networks has a natural translation in different modes of transportations connecting the same physical location in a city, a country, or on the globe.Finally, in the area of engineering and critical infrastructures, it applies to the interdependence of different lifelines. 317][38] For example, it has been shown that a diffusion process can have an enhanced-diffusive behaviour 7 on a multiplex network, which means that the time scale associated to it is shorter than that occurring on a single layer network.
All the already existing studies of random walks on multiplex networks adopt a nearest-neighbour navigation strategy. 39,40 e aim of this paper is to generalize Lévy flights random walks to multiplex networks, which means that a random walker has a certain probability to move to any other node without the need of a direct connection as far as the network is concerned.At each step, the random walker has three options: the first one is to stay at the same node, the second one is to jump to other nodes on the same layer and the last one is to switch to one of its counterparts on other layers, as illustrated in Fig. 1.According to the definition of the dynamics, we obtain the expression for the stationary distribution and the random walk centrality. 21esides, with the help of stochastic matrix theory, 17,41 we derive the exact expression for the mean first passage time (MFPT).The MFPT is used to describe the expected time needed for a random walker starting from a source point to reach a given target point. 42Finally, we also compare the results for Lévy flights with other random walks dynamics obtained also on multiplex networks, 39 finding that, under certain conditions, the Lévy random walk is the most efficient from a global viewpoint.

Results
In this work we consider undirected connected node-aligned multiplex networks. 29A node-aligned multiplex network is made up of L layers with N nodes i = {1, 2, . . ., N} on each layer.An adjacency matrix A α = {a α i j } N×N , with α = {1, 2, ..., L}, is associated to each layer α.Besides, a coupling matrix C = {c αβ i j } NL×NL describes the coupling between nodes in different layers; since each node is coupled only with its counterparts in different layers, then, only the elements of the type c αβ ii are different from zero.
The whole multiplex network can be described by the supra-adjacency matrix A = a i j NL×NL = A α +C. 29Additionally, we consider another set of matrices associated with the multiplex network, that is, we consider a distance matrix D α = {d α } N×N associated to each layer α, where the element d α i j encodes the length (number of steps) of the shortest path connecting node i to node j in layer α. 25 We indicate the probability to find a random walker in node j of layer β at time t starting from node i of layer α at t = 0 by p αβ i j (t).The discrete-time master equation is given: where w αβ i j is the transition probability of moving from node i of layer α to node j of layer β .To account for the inter-layer connections, we introduce D αβ ii to quantify the "cost" to switch from layer α to layer β at node i, while D αα ii quantifies the "cost" of staying in the same node and in the same layer.We can now define the transition probabilities w αβ i j to be where ii is the strength of node i with respect to its connections in the multiplex network, which takes into account the probability of staying at this node and of switching to another layer.As in the case of traditional single layer networks, 25 the transition probabilities w αβ i j define a dynamical process where a random walker can visit not only the nearest neighbours of a node, but also nodes without direct connections with it on the same layer, while the random walker can switch layer only staying at the same node.Since (d α i j ) −θ has an exponential decay according to the shortest path between the source node and the target node, the farther they are, the smaller the probability to hop to one from the other.The parameter θ , which varies in the range [0, ∞), controls the decay of this probability.
In the following we will derive the mean first passage time (MFPT), that is a characteristic quantity related to a random walk. 17By iterating Eq.( 1),we get an explicit expression for p αβ i j (t): Comparing according to the definition in Eq.( 2), we get For the stationary solution, which corresponds to the infinite time limit, we can get lim t→∞ p αβ i j (t) = p β j . 25Hence, Eq.( 4) implies that p β j s α i = p α i s β j and the probability p α i reduces to where s = ∑ α ∑ i s α i characterises the strength of the whole multiplex network.The expression of the stationary distribution p α i shows that the larger the strength of node i, the more often it will be visited, which is valid for any undirected network. 43he average of the MFPT over the stationary distribution is (see Methods for details of the derivation) whew λ k are the eigenvalues of the matrix W = w Besides, the random walk centrality of node i, as introduced in, 21 is given by (see Methods) Hence, we have derived the exact expression of the transition probability p αβ i j and the MFPT T .In addition, in order to analyse the navigation of Lévy walks, we average all the τ α i 's over the whole network, which means τ = 1 NL ∑ i,α τ α i .Note that with respect to the local index τ i , which represents the average time needed to reach node i from a randomly chosen node, τ gives the average number of steps needed to reach any node independently of the initial condition. 25ext, we proceed to characterise Lévy random walks on multiplex networks drawing on the exact analytic results given above.For the sake of simplicity, following Ref, 39 we assign the same value D X to all the D αβ ii , i.e., switching layers has the same cost at any node.In Fig. 2 we show τ vs θ for different topologies and different values of the cost D X .It is worth noting that, while for large values of the parameter D X the behaviour of the time τ when varying θ is qualitatively similar to the classical case of single layer networks, for small and intermediate values of D X it deviates significantly from the classical case.In particular, the relationship between τ and θ appears to be of three different kinds depending on D X : when D X is sufficiently small (D X = 0.1 in panel (a)) τ decreases quickly for small value of θ , while it remains more or less constant for large values of θ , ; when D X is sufficiently large (D X = 10 in panel (c)) τ increases monotonically with θ , as in classical single layer networks, with the speed of the increase being much smaller when θ is small.Furthermore, for intermediate values of D X ( D X = 1 in panel (b)) τ shows a clear minimum for a given value of θ .This phenomenology, that is, the fact that the efficiency depends on the coupling, constitutes the central finding of our study.
In single layer networks, setting θ = 0 is always the best strategy to navigate a network, as the global time τ is minimum.However, in multiplex networks the value of θ -it's optimal-that minimizes τ depends on the value of the coupling D X .Interestingly enough, for low values of D X , the limiting case θ → ∞, which corresponds to the normal random walks on networks, can be more efficient than θ = 0. Other scenarios worth inspecting are given by the topologies of the networks that made up each layer of the multiplex.In particular, a multiplex network can be made up of different combinations of homogeneous (Erdos-Renyi (ER)) and heterogeneous (Scale-Free (SF)) networks.We have also explored these scenarios numerically for different regimes of the coupling parameter.For D X = 1, different structures lead to different relationships between τ and θ .When θ is small, a SF-SF multiplex network has a much smaller τ than an ER-ER or an ER-SF multiplex network; however, if D X is bigger than 1, the difference is evident only when θ 0. Altogether, the previous results show that whether the optimal value of the Lévy walk index θ is constant across different multiplex topologies depends on the value of the coupling strength D X .

4/12
Next, in order to provide more numerical evidences of what we have found analytically, we present the results obtained when the fraction of covered nodes is taken into consideration.Figure 3 shows this magnitude as a function of time for different multiplex networks.In the first case (panel a), the network is made up of two ER networks with the same structure.For a small value of D X (upmost left figure), it can be seen that the bigger the value of θ is, the higher the efficiency of the Lévy random walk is.This also confirms the results in Fig. 2, since a larger value of θ leads to a smaller value of τ.With respect to other values of D X , such as D X = 1 (middle figures in all panels), D X = 10 (upmost right figures in all panels), the results for τ in Fig. 2 are also confirmed.In the case of other kinds of arrangements for the networks at each layer, we obtain the same results, as can be seen from panels b and c in Fig. 3. Furthermore, the comparison of the results obtained for different combinations shows that the topologies of the networks in each layer do not play a significant role.The previous results indicate that the coupling strength between layers is a crucial factor determining the structure and the dynamical behavior of the system. 45In addition, as described above, being used to characterize the cost for a walker to switch between layers, the value of D X also has distinct effects on Lévy random walks on multiplex networks.In order to further  103 ) in each layer.From left to right, the structure of network in each layer is ER-ER, ER-SF and SF-SF, respectively.The rectangles highlights the areas that show the largest differences due to the multiplex structure of the system.explore the details of these effects, as a function of θ , we show in Fig. 4 the time needed to cover the 50% of all nodes as a function of bothD X and θ .As shown in the figure, an interesting phenomenology appears.Firstly, the highest values of the time needed to cover half of the network locate at the up-right and down-left corners, where the values of D X and θ are the biggest and the smallest.Moreover, in a significant range of values of θ , increasing D X does not change greatly the time needed to cover 50% of the nodes.However, if θ is large enough, the increase of D X have a larger impact.Note that these results also confirm the analytical findings about τ, as can be seen in Fig. 2. Another result worth highlighting that connects our results for τ with the structural properties of the multiplex involves the second largest eigenvalue λ 2 .As shown for the smallest eigenvalue of the supra-Laplacian, 45 there is a transition point that separates two different regimes in interdependent networks: in one regime, all the layers are structurally decoupled and in the other regime, the system behaves as a single layer.The same result holds for the second smallest eigenvalue of the generalized supra-Laplacian of a multiplex network when increasing the coupling strength D X .Specifically, the generalized S F -S F Figure 6.Time needed to cover 50% of all nodes on three types of multiplex networks, as a function of the inter-layer weight D X .We compare the results obtained for the Lévy random walk (RWL) studied here with three other scenarios for the walks, as discussed in the text.The values of θ considered are 1, 5, 10, respectively.Each layer has 10 3 nodes and all the simulations were averaged over 100 realizations.

supra-Laplacian is
and I is the N × N identity matrix.
Figure5 shows the dependency of λ 2 with D X for different values of the Lévy flight parameter Θ.Also in this case (see 39 ), X , regardless of the network structure as showed it can be seen in the different panels of Fig. 5. Finally, for the sake of comparison with results obtained for other random walk dynamics, we compare their efficiency 39 with that of the Lévy flight.In the first case (RWC), the random walker in node i can move to any one of its neighbors j on the same layer with the transition probability w αα i j = , where k i is the degree of node i.Secondly, we also consider the case (RWD) in which the random walker is allowed to jump to any other node with probability w αα i j = s α i s max , where s max = max {i,α} s α i .Lastly, a third scenario (RWP) considers that it is possible for a random walker to switch layers and jump to another neighborhood at the same step.
In Fig. 6, we show the time τ needed to cover 50% of all the nodes as a function of the value of D X with respect to different topological structures.For the Lévy random walk, we study three different cases where the index θ equals 1, 5, 10, respectively.Comparing the Lévy case with the three others mentioned above, one can get further insights on the different strategies for navigation, that is to say, there is no strategy that is always the most efficient for any network and an arbitrary coupling strength.For instance, taking the Lévy random walk as an example: when θ is small (θ = 1), the time τ appears to be the smallest in 7/12 the range 1 < D X < 10, but as θ if further increased, the time needed to cover the 50% of the nodes of the network is almost the same as compared to that needed for a classical random walk, which in its turn is not the most efficient.This is easy to understand because in a Lévy random walk, when θ → ∞, the transition probability w i j = d −θ i j s i → 0 if the shortest path d i j is larger than 1.Therefore, the fist of the cases to which we compare -i.e., the classical random walk-is a special case of a Lévy walk when θ → ∞. 25

Discussion
In summary, we have studied Lévy random walks on multiplex networks.With the help of stochastic matrix theory, we have calculated analytically the expression of the stationary distribution and MFPT from any node to any other node.Besides, we have also obtained an exact expression of the average time τ needed to reach a node regardless of the source node.This dynamics on multiplex networks shows a strong dependence on the inter-layer weight D X and the Lévy index θ .Our main result is that when D X is small enough, contrary to the case of a Lévy random walk on single layer networks, the bigger the index θ is, the more efficient the Lévy random walk is.In order words, in that region of parameter values, although it is not very likely for any given walker to jump directly to other nodes far away, the total average time τ needed to visit any node independently of the initial condition is smaller.Interestingly, if the value of θ is not too large, for instance for θ < 4, D X does not have a significative impact on τ.The present results add to previous works that explored other kinds of random walk processes on multiplex networks, and allow to have a more complete picture that highlights the importance of considering the interconnected nature of many systems if we aim at finding the best navigation strategies and develop searching and navigability algorithms for such interdependent networked systems.

Methods
In the following, by using the formalism of generating functions, 44 we will get the analytical result for MFPT.The first passage probability q αβ i j (t) from node i of layer α to node j of layer β after t steps satisfies the relation Let T αβ i j denote the MFPT from node i of layer α to node j of layer β , then here, as proposed in, 17 we introduce the following generating functions: where |x| < 1, inserting Eq.8 into 11, we get Since , the problem of solving for the MFPT is reduced to calculate the derivative of Qαβ i j (x) and evaluate it at x = 1.
We will address this point making use of the stochastic matrix theory. 41For the sake of simplicity, we use the matrix W = w αβ i j NL×NL and the matrix S to describe the transition probabilities and node 8/12 strengths, respectively.It is clear that the matrix W is a stochastic matrix, since for any node i its elements satisfy that ∑ NL j=1 w i j = 1.W is an antisymmetric matrix.Because of that, we introduce the matrix which is symmetric and similar to W. Thus, they have the same eigenvalues.Since Γ can be diagonalized and the eigenvalues are all real, we define λ 1 , λ 2 , • • • λ NL as its eigenvalues.These eigenvalues satisfy that 1 denote the corresponding normalized, real-valued, and mutually orthogonal eigenvectors.As a result, the matrix Γ can be written as which, together with 13, leads to Then considering the master equation Eq.( 1), we can get P(t) = P(0)W t = S − 1 2 Γ t S 1 2 , whose element denoted by p αβ i j (t) represents the transition probability from node i of layer α to node j of layer β in t steps.Note that the elements of the matrix P(0) fulfill the following relations Then, inserting Eq.15 into the expression of P(t), one has where φ i = (φ i1 , φ i2 , • • • φ iNL ) T and they satisfy φ T i φ j = 1 if i = j, else φ i φ j = 0, which means According to the definition given above, the MFPT T αβ i j can be calculated by differentiating Qαβ i j (x) In addition, using Eq.20, we have where the time T is the average of the MFPT over the stationary distribution, obviously it does not depend on i, and it is known as the Kemeny's constant. 25Besides, as introduced in, 21 we can calculate the quantity C α i = (τ α i ) −1 , that is the random walk centrality of node i, where τ α i is defined as {p αα ii (t) − p α i }/p α i .Combining Eq. 17, τ α i is given by (23)

Figure 1 .
Figure 1.Illustration of the Lévy flight navigation strategy on a multiplex network.In the toy model, we consider a three-layer multiplex network and show two different paths that can lead the walker to the yellow node (one involves a Lévy flight and the other implies that the walker follows the topological path of the graph).The right panel summarizes the different elementary steps that a walker can adopt in our model as indicated.

Figure 2 .
Figure 2. The quantity τ vs the Lévy flight index θ for different two-layer multiplex networks.Each layer is an ER network or a SF network with 1000 nodes as indicated in the legends.The values of the coupling strength D X between the two layers are: (a) D X =0.1,(b) D X =1, (c) D X =10.

Figure 3 .
Figure 3. Number of visited nodes versus time for Lévy random walks on multiplex networks.The structures of the multiplex networks considered are: ER-ER (top panels a), ER-SF (middle panels b) and SF-SF (bottom panels c).In each configuration, the synthetic networks of each layer are composed of 10 3 nodes.From left to right, the values of D X are 0.1, 1, 10.The Lévy index θ used are indicated in the legend.

Figure 4 .
Figure 4.The effects of inter-layer weight D X and Lévy index θ on the efficiency of Lévy random walks.The color-coded map describes the time needed to cover 50% of all the nodes (103 ) in each layer.From left to right, the structure of network in each layer is ER-ER, ER-SF and SF-SF, respectively.The rectangles highlights the areas that show the largest differences due to the multiplex structure of the system.

Figure 5 .
Figure 5.The second smallest eigenvalue of the generalized supra-laplacian matrix as a function of the inter-layer weight D X for three different multiplex topologies, which from left to right are ER-ER, ER-SF, SF-SF, respectively.Each panel describes Lévy random walk with different Lévy index θ (1, 5, 10).The solid line corresponds to D −1 X .

φ ki φ k j = 0 ( 18 )
We have now an expression for p αβ i j (t), plugging it into Eq.11, it is easy to obtain Pαβ i j (x) = Eq.18 and the relation φ 1i φ 1 j s