Effect of shortest path multiplicity on congestion of multiplex networks

Shortest paths in graphs are defined as paths between two nodes such that the sum of the weights of the links composing the path is minimized [1]. The computation of shortest paths is of utmost importance to determine the efficiency of the network to exchange information [2], to compute the load of nodes (defined as the number of shortest path traversing it) [3–6], to predict and alleviate congestion [7–10] or to classify networks [11], to mention some. The main applications of the finding of shortest paths are routing strategies [12, 13], analysis of road networks [14], epidemic spreading [15], or analysis of brain activity [16], among others.

The development of efficient numerical algorithms for both the exact and approximate calculation of shortest paths is a problem by itself that still attracts the attention of computer scientists [17]. The exact numerical calculation of all the shortest paths from a single node is commonly addressed using the well-known Dijkstra's algorithm whose running time is $O(M+N\mathrm{log}N)$, where N is the number of nodes and M is the number of edges of the graph. Given that, the numerical calculation of shortest paths in large networks is still computationally expensive.

In particular, the computation of shortest paths is essential for the determination of a fundamental measure of centrality: edge and node betweenness. The node (or link) betweenness is normally defined as the fraction of shortest paths between node pairs that pass through the node (or link) of interest. The ranking according to betweenness [18] is informative about the bottlenecks of a certain network structure, and can be used to determine the most critical node or link with respect to any process traversing the network using shortest paths. In particular, the node with the largest betweenness in a network defines the onset of congestion [7, 19].

The extension of centrality descriptors to multiplex networks [20–22] prompts for the quantification of the distribution of shortest paths in this new scenario [23, 24]. Multiplex networks are defined as a set of interconnected layers of networks, where every node has its own replica through all layers, and is connected with them [25–28]. Shortest paths in multiplex networks are essentially governed by the shortest paths at each layer, but also by those emerging from partial paths at different layers and the switch between layers. Here, we present a method to analytically determine the distribution of shortest paths in multiplex networks, with special emphasis on the multiplicity of shortest paths occurring in the mentioned set up. This computation allows us to determine semi-analytically the feasible region of congestion induced by the multiplex structure [29].

The paper is structured as follows. First, we present the definition of shortest paths in multiplex networks. Next, we show how to compute the distribution of shortest paths in two-layer multiplex networks. Once we provide a way to compute the joint distribution of shortest paths and their multiplicities in sigle layer Erdős–Rényi networks, we extend the method to multiplex network, and evaluate the accuracy of our analytical predictions. Finally, we apply our results to the analysis of congestion in multiplex networks, being able to predict the region of parameters in which the multiplex structure may induce congestion.

Multiplex networks may, at first sight, seem equivalent to single layer networks, just with the difference of having two kinds of edges: intralayer links between nodes in the same layer, and interlayer links connecting the replicas of each node in the different layers (see figure 1). However, the fact that the replicas of each node refer exactly to the same entity, has important consequences on the structural and dynamical properties of multiplex networks.

Zoom In Zoom Out Reset image size

Suppose, for example, we have the transportation network of a city, formed by metro and bus connections. Nodes represent the locations of the metro stations and bus stops, and interlayer links allow the transfer from metro to bus and vice versa. There exists a cost associated with each link, which accounts for the time needed to travel from one location to another (intralayer links), and the time spent to make a transfer (interlayer links). A path in this network must take into account both types of cost, thus we could use any of the standard methods to find the shortest paths in graphs. Now, consider that we want to go from a source location A to a destination location B taking the shortest route, and that we may start the trip either with bus (A1) or metro (A2).The algorithm provides us with four different types of paths: two in which the origin and destination are both in the same layer, $A1\to B1$ and $A2\to B2;$ and two where the endpoints belong to different layers, $A1\to B2$ and $A2\to B1$. Note that neither of the paths are necessary contained within a single layer, e.g. a shortest path between A1 and B1 could change layer twice. Since we are only interested in going from A to B, and it is irrelevant if we start or end the trip in the bus or metro, some of these shortest paths should be discarded.For example, if the length of the shortest paths $A1\to B1$ is larger than that for $A1\to B2$, we should discard ending our traversal at layer 1. These means that, in two-layer multiplex networks, shortest paths between two nodes may exist in just one layer, in the other layer, in a path with transfers (what we call multiplex paths), or in a combination of them (see figure 1). Consequently, any property depending on paths in multiplex networks must consider the concept of source and destination as being intrinsically different to that of nodes, such as in the calculation of centrality measures [22–24], the analysis of interdependence [30], the analysis of random walks [21], or the study of congestion [29]. The former definition of shortest paths in multiplex networks imposes the necessity of correctly computing their distribution and multiplicity.

The problem we are going to address is how the distribution of shortest paths changes when we take two different networks (the layers), and connect them to build a multiplex network. The prescription we propose can be easily extended to more layers, however its combinatorics makes the mathematical analysis more involved. In particular, we are interested in two important parameters [29]: λ, the fraction of all shortest paths fully contained in one layer (i.e. paths which do not make use of the multiplex option of changing layer using interlayer links); and μ_α, the fraction of non-multiplex shortest paths using only layer α. Thus, ${\sum }_{\alpha }{\mu }_{\alpha }=1$, and 1 − λ amounts for the fraction of multiplex shortest paths. It is important to remark that shortest paths are usually degenerated, i.e. there are many shortest paths with the same length, with some of them being multiplex paths; we will refer to this degeneration as the multiplicity of the shortest paths.

Previous works have provided methods to estimate the complementary cumulative probability F(d) that two randomly selected nodes are separated by a distance larger than d, for certain classes of monoplex (single layer) random networks [31–35]. Thus, the probability distribution f(d) of two nodes being exactly at distance d can be expressed in terms of the complementary cumulative probability as

$$\begin{eqnarray}&&f(d)=F(d-1)-F(d).\end{eqnarray} \tag{ 1 }$$

We denote by f₁(d), f₂(d) and f_M(d) the respective distributions of shortest paths fully contained in layer 1, in layer 2, or using the full multiplex structure. Assuming an unweighted and undirected multiplex network with two layers and no degree correlations, we may use this relation to obtain the fraction of multiplex shortest paths:

$$\begin{eqnarray}\begin{array}{rcl}1-\lambda & = & \displaystyle \sum _{d=3}^{\infty }{f}_{M}(d)\left[{F}_{1}(d){F}_{2}(d)+{\theta }_{M1}(d){f}_{1}(d){F}_{2}(d)+{\theta }_{M2}(d){F}_{1}(d){f}_{2}(d)\right.\\ & & +\left.{\theta }_{M12}(d){f}_{1}(d){f}_{2}(d)\right].\end{array}\end{eqnarray} \tag{ 2 }$$

The sum starts with d = 3, since the minimum length of a multiplex path is 3, i.e. at least one hop in the first layer, another in the second layer, and one to perform the change of layer. Equation (2) expresses the four possible ways in which shortest paths can appear in the multiplex: only multiplex paths; multiplex paths mixed with paths in layer 1; multiplex paths mixed with paths in layer 2; and multiplex paths mixed with paths in layers 1 and 2.  Consider for example the case of multiplex shortest paths mixed with shortest paths in layer 1. This correspond to the paths that have exactly distance d in layer 1, f₁(d), and in the multiplex, f_M(d), and distance larger than d in layer 2, F₂(d). For the terms with mixed contributions, the θ(d) factors capture the fraction of shortest paths that correspond to multiplex paths. This means we cannot calculate λ just with the knowledge of the shortest paths distributions f(d) and F(d), we also need to estimate the multiplicity of the paths of each type.

Let us denote by P₁(d, ϕ), P₂(d, ϕ) and P_M(d, ϕ) the probabilities that randomly chosen source and destination (discarding layer, i.e. locations in our example) are exactly at distance d with multiplicity ϕ, for paths in first layer, second layer, and in the multiplex, respectively. The average multiplicities $\langle \phi {\rangle }_{1}(d)$, $\langle \phi {\rangle }_{2}(d)$ and $\langle \phi {\rangle }_{M}(d)$ for each kind of path can be expressed as:

$$\begin{eqnarray}&&\langle \phi {\rangle }_{\alpha }(d)=\displaystyle \frac{{\displaystyle \sum }_{\phi =1}^{\infty }\phi \,{P}_{\alpha }(d,\phi )}{{\displaystyle \sum }_{\phi =1}^{\infty }{P}_{\alpha }(d,\phi )},\quad \alpha \in \{1,2\},\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&\langle \phi {\rangle }_{M}(d)=\displaystyle \frac{{\displaystyle \sum }_{\phi =1}^{\infty }\phi \,{P}_{M}(d,\phi )}{{\displaystyle \sum }_{\phi =1}^{\infty }{P}_{M}(d,\phi )}.\end{eqnarray} \tag{ 3b }$$

Therefore,

$$\begin{eqnarray}&&{\theta }_{M\alpha }(d)=\displaystyle \frac{\langle \phi {\rangle }_{M}(d)}{\langle \phi {\rangle }_{M}(d)+\langle \phi {\rangle }_{\alpha }(d)},\quad \alpha \in \{1,2\},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{\theta }_{M12}(d)=\displaystyle \frac{\langle \phi {\rangle }_{M}(d)}{\langle \phi {\rangle }_{M}(d)+\langle \phi {\rangle }_{1}(d)+\langle \phi {\rangle }_{2}(d)}.\end{eqnarray} \tag{ 4b }$$

Additionally, the distributions f(d) can be recovered as marginals of P(d, ϕ):

$$\begin{eqnarray}&&{f}_{\alpha }(d)=\displaystyle \sum _{\phi =1}^{\infty }{P}_{\alpha }(d,\phi ),\quad \alpha \in \{1,2\},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{f}_{M}(d)=\displaystyle \sum _{\phi =1}^{\infty }{P}_{M}(d,\phi ).\end{eqnarray} \tag{ 5b }$$

Similarly to (2), we may use the distributions F(d) and f(d), and the multiplicity factors θ(d), to calculate how the non-multiplex shortest paths are distributed among the layers. The expression for the first layer reads:

$$\begin{eqnarray}&&{\mu }_{1}=\displaystyle \sum _{d=1}^{\infty }{f}_{1}(d)\left[{F}_{2}(d)+{\theta }_{1}(d){f}_{2}(d)\right],\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{\theta }_{1}(d)=\displaystyle \frac{\langle \phi {\rangle }_{1}(d)}{\langle \phi {\rangle }_{1}(d)+\langle \phi {\rangle }_{2}(d)}.\end{eqnarray} \tag{ 7 }$$

Equivalent expressions hold for μ₂, and μ₁ + μ₂ = 1.

Summarizing, we have set all the necessary ingredients to analyze the use of the multiplex structure and of the layers in terms of the joint distributions P₁(d, ϕ), P₂(d, ϕ) and P_M(d, ϕ). In the next sections we develop expressions for each of these probabilities.

The analytical calculation of the joint probability that two nodes are exactly at distance d with multiplicity ϕ, P(d, ϕ), for any type of monoplex networks, is quite involved, thus we are going to restrict our analysis to Erdős–Rényi networks. However, it would be possible to try with other types of uncorrelated random networks making use of hidden variables as in [31].

An Erdős–Rényi network is characterized by two parameters, the number of nodes N, and the probability p that an edge exists between any pair of nodes. Thus, we will use the notation P_{p, N}(d, ϕ) to refer to the joint distance and multiplicity distribution for these networks. Clearly, at d = 1,

$$\begin{eqnarray}&&{P}_{p,N}(1,1)=p,\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&{P}_{p,N}(1,\phi )=0,\quad \forall \phi \gt 1,\end{eqnarray} \tag{ 8b }$$

since we do not allow multiple edges between pairs of nodes. The case d = 2 can be obtained as the probability of not having a shortest path of length 1 times the probability of having exactly ϕ different shortest paths of length 2. Since paths of length 2 have just one intermediate node, we have to choose ϕ nodes among the N − 2 available to build the ϕ shortest paths. Formally

$$\begin{eqnarray}&&{P}_{p,N}(2,\phi )=(1-p)\displaystyle \left(\genfrac{}{}{0em}{}{N-2}{\phi }\right){p}^{2\phi }{\left(1-{p}^{2}\right)}^{(N-2)-\phi }.\end{eqnarray} \tag{ 9 }$$

Note that the existence of paths of length 2 has probability p², thus p^2ϕ expresses the probability of having ϕ paths of length 2, and ${\left(1-{p}^{2}\right)}^{(N-2)-\phi }$ the probability that the rest of the (N − 2) − ϕ nodes are not used to build shortest paths of length 2. Using the same procedure, the generalization to larger distances becomes

$$\begin{eqnarray}&&{P}_{p,N}(d,\phi )=\left[(1-p)\displaystyle \prod _{k=2}^{d-1}{\left(1-{p}^{k}\right)}^{\eta (k)}\right]B(\phi ;\eta (d),{p}^{d}),\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&\eta (k)=\displaystyle \prod _{i=0}^{k-2}(N-2-i)\end{eqnarray} \tag{ 11 }$$

stands for the amount of possible paths of length $k\geqslant 2$, and

$$\begin{eqnarray}&&B(\phi ;\eta (d),{p}^{d})=\displaystyle \left(\genfrac{}{}{0em}{}{\eta (d)}{\phi }\right){p}^{d\phi }{\left(1-{p}^{d}\right)}^{\eta (d)-\phi }\end{eqnarray} \tag{ 12 }$$

is the probability density function of the binomial distribution with probability p^d. The term within square brackets in (10) computes the probability of not having any path at distance lower than d.  Although (9) is exact, (10) is an approximation because it does not consider that some shortest paths may share some (but not all) of their links.

We now address the computation of the joint distribution P_M(d, ϕ) for a multiplex network composed of two Erdős–Rényi layers, with N nodes each, edge probabilities p₁ and p₂, and without interlayer degree correlations; we will refer to it as ${P}_{{p}_{1},{p}_{2},N}(d,\phi )$, to emphasize the dependence on these structural parameters. Multiplex shortest paths are characterized by the number of interlayer links they contain, i.e. how many times the path changes from one layer to the other, see figure 2. We restrict our analysis to shortest paths that change layer only once, since the number of shortest paths with multiple jumps between layers is usually very small, thus making their contribution negligible, and the generalization of the mathematical formulation for two or more jumps is not difficult but generates very long expressions.

Zoom In Zoom Out Reset image size

For each multiplex shortest path which has just one change of layer, as in figure 2(a), once we have selected the origin and destination, we may choose the intermediate switch node among the remaining N − 2. Since we are considering a multiplicity ϕ of multiplex shortest paths between these origin and destination nodes, the paths may be distributed in many different ways. For example, if ϕ = 10, we could select 7 shortest paths to use one intermediate switch node, and the remaining 3 to choose another one. The set of all structurally different distributions of one-jump multiplex shortest paths is given by the partition set of the given integer multiplicity ϕ. The finding and counting of the number of partitions of an integer number constitutes a classical problem in number theory [36], e.g. they can be enumerated using Young diagrams [37]. We may express the set of available partitions as:

$$\begin{eqnarray}&&{ \mathcal F }(\phi )=\{{\boldsymbol{\Phi }}\in \bigcup _{r=1}^{\phi }\{1,\ldots ,\phi \}{}^{r}:\displaystyle \sum _{i=1}^{| {\boldsymbol{\Phi }}| }{\phi }_{i}=\phi ,\ {\phi }_{1}\geqslant \,\cdots \,\geqslant {\phi }_{| {\boldsymbol{\Phi }}| },\ | {\boldsymbol{\Phi }}| \leqslant N-2\}.\end{eqnarray} \tag{ 13 }$$

For example, ϕ = 4 can be partitioned in five different ways: ${ \mathcal F }(4)=\{(4)$, (3, 1), (2, 2), (2, 1, 1), (1, 1, 1, 1)}. If N = 5, we would have only 3 possible intermediate nodes, and partition (1, 1, 1, 1) should be discarded, hence the $| {\boldsymbol{\Phi }}| \leqslant N-2$ condition added to (13).

Once a partition Φ is selected, we need to consider the different ways in which we can choose the intermediate nodes, and the different ways in which the multiplicities ϕ_i are assigned to the selected set of intermediate nodes. Gathering together all these contributions, we may write the joint distribution of multiplex shortest paths as:

$$\begin{eqnarray}&&{P}_{{p}_{1},{p}_{2},N}(d,\phi )=\displaystyle \sum _{{\boldsymbol{\Phi }}\in { \mathcal F }(\phi )}\displaystyle \left(\genfrac{}{}{0em}{}{N-2}{| {\boldsymbol{\Phi }}| }\right)C({\boldsymbol{\Phi }})\left[{R}_{{p}_{1},{p}_{2},N}{\left(d\right)}^{N-2-| {\boldsymbol{\Phi }}| }\displaystyle \prod _{i=1}^{| {\boldsymbol{\Phi }}| }{Q}_{{p}_{1},{p}_{2},N}(d,{{\rm{\Phi }}}_{i})\right],\end{eqnarray} \tag{ 14 }$$

where the combinatorial number accounts for the selection of intermediate switch nodes, $C({\boldsymbol{\Phi }})$ for the assignment of individual multiplicities, ${R}_{{p}_{1},{p}_{2},N}(d)$ stands for the probability of not having multiplex paths of length lower or equal than d, and ${Q}_{{p}_{1},{p}_{2},N}(d,{{\rm{\Phi }}}_{i})$ stands for the joint probability of having one-switch multiplex shortest paths with known intermediate jump node. Factor C(Φ) is just the multinomial coefficient corresponding to the frequencies of the components of Φ. For example, if Φ = (14, 3, 3), there are three components, the 14 with frequency 1 and the 3 with frequency 2, leading to C(Φ) = 3!/(1! 2!) = 3, which corresponds to the 3 ways in which we can sort the components of Φ, i.e. (14, 3, 3), (3, 14, 3) and (3, 3, 14).

Probability ${R}_{{p}_{1},{p}_{2},N}(d)$ admits a simple expression:

$$\begin{eqnarray}&&{R}_{{p}_{1},{p}_{2},N}(d)=1-\displaystyle \sum _{s=1}^{d}\displaystyle \sum _{\phi =1}^{\infty }{Q}_{{p}_{1},{p}_{2},N}(s,\phi ),\end{eqnarray} \tag{ 15 }$$

However, ${Q}_{{p}_{1},{p}_{2},N}(d,\phi )$ still requires further decompositions. Although we know the origen, destination and intermediate nodes, the total length d of the paths, and the total multiplicity ϕ, it remains to establish: the layer α at which the path starts; the length ℓ covered in layer 1 (the rest of the path in layer 2 will have length d − ℓ −1); the distribution of the multiplicities between the two layers. Moreover, when ϕ > 1, we may have a combination of paths starting at different layers, with different lengths in each layer, and with different distribution of the multiplicities per layer. Denoting by ψ_α,_{ℓ, β} the multiplicity of paths in layer β, when they start in layer α and have length ℓ in layer 1, the set of possible decompositions of the multiplicity can be written as:

$$\begin{eqnarray}&&{ \mathcal G }(d,\phi )=\{{\boldsymbol{\Psi }}\in \{0,\,\ldots ,\,\phi \}{}^{2\times (d-2)\times 2}:\displaystyle \sum _{\alpha =1}^{2}\displaystyle \sum _{{\ell }=1}^{d-2}{\psi }_{\alpha ,{\ell },1}{\psi }_{\alpha ,{\ell },2}=\phi \}.\end{eqnarray} \tag{ 16 }$$

Note that a paths arriving to the switching node in one layer, and b paths departing from it in the other layer, generate ab different paths in the multiplex, hence the products inside the sums in (16). An example of a decomposition Ψ of ${ \mathcal G }(5,10)$ could be: 2 paths of length 2 starting in layer 1, followed by 3 paths of length 2 in layer 2 (amounting 6 paths), plus 1 path of length 1 starting in layer 2, followed by 4 paths of length 3 in layer 1 (amounting 4 paths), which correspond to ψ_1,2,1 = 2, ψ_1,2,2 = 3, ψ_2,3,2 = 1 and ψ_2,3,1 = 4 (the rest of the components are zero). In practical terms, the calculation of ${ \mathcal G }(d,\phi )$ involves: finding the partitions of ϕ; factorizing the parts as products of two integers, by calculating all their divisors; choosing initial layer; and dividing the total length as the sum of the lengths in each layer plus one.

Using the decomposition set in (16), we may finally calculate the remaining probabilities in (14) and (15):

$$\begin{eqnarray}&&{Q}_{{p}_{1},{p}_{2},N}(d,\phi )=\displaystyle \sum _{{\boldsymbol{\Psi }}\in { \mathcal G }(d,\phi )}\mathop{\displaystyle \prod _{\alpha =1}^{2}\displaystyle \prod _{{\ell }=1}^{d-2}}\limits_{{\psi }_{\alpha ,{\ell },1}{\psi }_{\alpha ,{\ell },2}\ne 0}{P}_{{p}_{1},N}({\ell },{\psi }_{\alpha ,{\ell },1}){P}_{{p}_{2},N}(d-{\ell }-1,{\psi }_{\alpha ,{\ell },2}).\end{eqnarray} \tag{ 17 }$$

Note the presence of the joint probabilities (10) for single layer Erdős–Rényi networks, which account for the subpaths of the multiplex shortest paths fully contained in each of the layers.

To check the validity of our calculations of the distribution of paths in multiplex networks, we have generated a large set of two-layer multiplex networks with uncorrelated Erdős–Rényi layers, and compared the predicted values with the experimental ones. In particular, we have analyzed multiplex networks with 500 nodes in each layer, for a total of 50² different configurations. These configurations correspond to 50 logarithmically spaced values of the average degree $\langle k\rangle$ of each layer, with values ranging between 5 and 35. For each configuration, we generate 100 different multiplex networks, calculate their experimental values of λ and μ₁, and take averages. These averages are then compared with the predicted values using (2) and (6), respectively. The analytical values only consider multiplex shortest paths with just one change of layer (as discussed above), and contributions from multiplicities larger than 100 have been discarded. The results are presented in figure 3, which shows an excellent agreement between theory and experiments.

Zoom In Zoom Out Reset image size

When the average degree of an Erdős–Rényi network is large, its average shortest path length is small [31], thus it is difficult to find multiplex shortest paths since the overhead of changing layer has to be compensated with very short paths in the layers. As a consequence, the fraction λ of shortest paths fully included in one of the layers tends to 1, as shown in figure 3(a). If we reduce these average degrees, the shortest paths of the layers become larger, and the opportunities to find multiplex shortest paths increases, yielding lower values of λ. In the extreme cases in which the average degrees are very small, multiplex shortest paths become more common, and it is even possible to find shortest paths with two or more changes of layer; this is the reason for the small deviations between the experimental and predicted values of λ for small $\langle {k}_{1}\rangle +\langle {k}_{2}\rangle$.

Figure 3(b) shows also that non-multiplex shortest paths tend to be concentrated in the layer with largest average degree, i.e. with smaller average shortest path length. Thus, parameter μ₁ approaches value 1 when $\langle {k}_{1}\rangle \gg \langle {k}_{2}\rangle$, and 0 when $\langle {k}_{1}\rangle \ll \langle {k}_{2}\rangle$.

Shortest paths play an important role in congestion phenomena in complex networks. When elements (packages, vehicles, etc.) travel a network using shortest paths, the load of the nodes is directly related to the number of shortest paths that make use of them. If the capacity of some nodes to process these elements is lower than their incoming rate, congestion emerges [7]. In the general case of multiplex networks [29], it was shown that congestion appears when the injection rate per node, ρ, reaches a critical value ρ_c:

$$\begin{eqnarray}&&{\rho }_{c}=\tau {L}^{-1}\displaystyle \frac{N-1}{{{ \mathcal B }}^{* }},\end{eqnarray} \tag{ 18 }$$

where τ is the processing capacity of the nodes, and ${{ \mathcal B }}^{* }$ is the maximum betweenness of all nodes in all layers of the multiplex network. Additionaly, as it is shown in [29], the maximum betweenness found in the multiplex can be approximated in terms of the betweenness of the most efficient layer, ℓ, as ${{ \mathcal B }}^{* }\approx \lambda {\mu }_{{\ell }}{{ \mathcal B }}_{{\ell }}^{* }$. Specifically, ℓ is the layer for which the onset of congestion is larger than for the rest of the layers, when layers are considered as independent networks. Thus, the onset of congestion ρ_c becomes

$$\begin{eqnarray}&&{\rho }_{c}\approx \displaystyle \frac{1}{L\lambda {\mu }_{{\ell }}}{\rho }_{c}^{({\ell })},\end{eqnarray} \tag{ 19 }$$

where ${\rho }_{c}^{({\ell })}$ is the critical injection rate of layer ℓ, given by (see [7])

$$\begin{eqnarray}&&{\rho }_{c}^{({\ell })}=\tau \displaystyle \frac{N-1}{{{ \mathcal B }}_{{\ell }}^{* }}.\end{eqnarray} \tag{ 20 }$$

For example, if layer one is the most efficient, ℓ = 1, and we have two layers, L = 2, the first layer of the multiplex accepts much more load than the second layer before congestion appears, thus ${\rho }_{c}^{(1)}\gt {\rho }_{c}^{(2)}$. When we connect these two layers to form a multiplex, there exists a migration of shortest paths from the second to the first layer, which increases its load, eventually becoming responsible for the onset of congestion of the multiplex. This increases the congestion of the most efficient layer, as shown in (19). An important consequence of this redistribution of shortest paths in the multiplex, and the appearance of multiplex shortest paths, is that the multiplex may attain congestion with lower load than for any of its separated layers, the so-called congestion induced by the multiplex structure [29]. In our example, this happens when ${\rho }_{c}\lt {\rho }_{c}^{(2)}$. By combining this inequality with (19) and (20), we obtain the condition for having congestion induced by the multiplex:

$$\begin{eqnarray}&&L\lambda {\mu }_{1}\leqslant \displaystyle \frac{{{ \mathcal B }}_{2}^{* }}{{{ \mathcal B }}_{1}^{* }}.\end{eqnarray} \tag{ 21 }$$

We have shown above how to analytically estimate the values of λ and μ₁, by using (2) and (6), respectively. However, the estimation of ${{ \mathcal B }}_{\alpha }^{* }$ in terms of the structural properties of the network is not straightforward.We have experimentally observed that the relationship between the maximum betweenness (excluding endpoints) and $\langle {k}_{\alpha }\rangle$ follows a power law with exponent −1, with a constant proportional to N − 1 and to the average multiplicity found in Erdős–Rényi networks of the same size, $\langle \phi {\rangle }^{* }$. Thus, we may express the maximum betweenness as

$$\begin{eqnarray}&&{{ \mathcal B }}_{\alpha }^{* }\approx (N-1)\langle \phi {\rangle }^{* }\langle {k}_{\alpha }{\rangle }^{-1}+2(N-1),\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}&&\langle \phi {\rangle }^{* }=\mathop{\max }\limits_{2\leqslant \langle k\rangle \leqslant N}\displaystyle \frac{1}{{d}_{\langle k\rangle }}\displaystyle \sum _{d=1}^{{d}_{\langle k\rangle }}\langle \phi {\rangle }_{\langle k\rangle }(d).\end{eqnarray} \tag{ 23 }$$

The term ${d}_{\langle k\rangle }$ corresponds to the maximum expected distance in an Erdős–Rényi network with average degree $\langle k\rangle$, and $\langle \phi {\rangle }_{\langle k\rangle }(d)$ to the average multiplicity in (3a) for that network.

We show in figure 4 the estimated region of the parameters space in which congestion is predicted to be induced by the multiplex structure, for multiplex networks composed of Erdős–Rényi layers, in good agreement with the experimental results.

Zoom In Zoom Out Reset image size

We have presented a method to compute the multiplicity and distribution of shortest paths in multiplex networks. Using the method, we have analytically determined the distribution and multiplicity of shortest paths in duplex of Erdős–Rényi networks. This results are essential to determine the onset of congestion on multiplex structures, as for example in many multimodal transportation networks. We can use the analytical findings to determine the area of the phase diagram where the multiplex can induce congestion. This work is relevant for any dynamical process that uses shortest path as a routing strategy in multilayer networks.

This work has been supported by Ministerio de Economía y Competitividad (grant FIS2015-71582-C2-1), Generalitat de Catalunya (grant 2017SGR-896), and Universitat Rovira i Virgili (grant 2017PFR-URV-B2-41). AA acknowledges partial financial support from ICREA Academia and the James S McDonnell Foundation.