Unique superdiffusion induced by directionality in multiplex networks

Multilayer networks provide a proper mathematical representation of complex systems with different types of interactions, e.g. social interactions among individuals across different social platforms such as Facebook and Twitter, and enable the understanding of dynamics acting on or evolving within those systems [1, 2]. A number of dynamical processes are studied on the framework of multilayer networks: for example, innovation/information diffusion [3], the control of formation and task allocation of swarms [4–6], synchronization of oscillators [7], the understanding of functional brain connectivity [8], and many others. More specifically, recent studies based on a multilayer network representation have revealed unforeseen dynamical regimes and behaviors, such as new spreading regimes when more than one diseases in spreading in a given population [9–12], and enhanced stability of synchronized states [13]. Moreover, conversely from building up a multilayer structure, gradually dismantling a multilayer network triggers an abrupt transition [14], suggesting a non-additive effect of a multilayer structure from the integration of its layers. In terms of diffusion dynamics, multilayer network studies revealed that for any multilayer configuration, the diffusive behavior of the overall system (multilayer network) could be tuned to be faster than that of the slowest layer for significant values of interlayer diffusivity (in comparison to intralayer diffusivity). More surprisingly, but only for certain multilayer network configurations, the overall system can exhibit a superdiffusive behavior for large values of interlayer diffusivity, where the diffusive dynamics on the multilayer network is faster than in any of the individual layers, when those are considered in independently [15–17].

Most of the studies made in the field of multilayer networks adopt undirected structures, i.e. the interactions represented by both interlayer and intralayer links are assumed to be undirected. However, most real world systems are inherently directed. Examples are the World Wide Web, in which hyperlinks run in one direction from one Webpage to another; food webs, in which energy flows from prey to predator; phone call networks; metabolic networks; citation networks; social networks of followers and followee; gene regulation networks, to name a few. Additionally, some dynamical processes are governed by physical laws and inherently evolve in a directional manner despite of the underlying undirected connectivity topology, like current/information flow from high voltage to low voltage, virus transmissions from infected to susceptible individuals. Moreover, multi-leveled unidirectional interactions are commonly found in interacting systems. For example, a gene regulation network aims to properly characterize both protein-DNA and protein-protein directed interactions between genes and proteins [18]. Therein, multilayer networks with directionality provide a realistic representation for such systems consisting of multilevel and directed interactions.

Multilayer networks where at least one of the layers is directed give rise to new theoretical challenges in analyzing their dynamical behaviors. As diffusion dynamics on undirected multilayer networks is fully characterized by the spectrum of network matrices, like supra-Laplacian [15], network directionality leads to nonsymmetrical graph matrices, which mostly features complex and non-orthogonal spectrum. Consequently, spectrum related analysis on undirected multilayer networks might be proven as invalid or inaccurate. Few recent studies on directed multilayer networks reveal new and unexpected physical behaviors [16, 19], as well as increased sensitivity to structural changes [20].

The contribution of this paper is to provide a theoretical framework for diffusive dynamics on multilayer networks with directed layers. We show that although certain undirected multilayer networks can result in enhanced diffusive (faster than the slowest layer) or superdiffusive (faster than any of the layers) regimes for large values of the interlayer coupling, the diffusive rate of the overall system never exceeds the diffusion on the equivalent aggregated system (the direct sum of each layer). We also show that certain directed multilayer networks can also show enhanced diffusive and superdiffusive regimes for large values of interlayer coupling. Remarkably, directed multilayer networks, differently from their undirected counterparts, can exhibit a unique superdiffusive behavior for intermediate values of coupling, where diffusion dynamics are even faster than the corresponding aggregated system. This unique superdiffusion regime only emerges when certain structural conditions are satisfied. We analytically and numerically uncover conditions driving the unique superdiffusion in directed multilayer networks. In addition, we propose a model to configure directed multilayer networks that achieves the unique superdiffusion with tunable magnitude.

The paper is organized as follows. Diffusive dynamics on multilayer networks with direction is described in the Model section 2. The section 3 of Results presents the conditions for diffusion regimes below the diffusivity in the corresponding aggregated system in section 3.1 and for the superdiffusion regime above the diffusivity in the aggregated system in section 3.2, followed by a model for the construction of a multilayer network with superdiffusion in section 3.3. Section 4 concludes the work.

Multilayer networks consist of nodes interacting both within the same layer and across different layers via intralinks and interlinks which encodes multiple types of interactions. In this study, we focus on a particular type of multilayer networks called multiplex networks, characterized by layers consisting of the same set of nodes, but possibly different connectivity (layer topology); and layers interacting with each other only via counterpart node. The topological structure of the multiplex is described by the so-called supra-adjacency matrix $A=\left({a}_{ij}\right)$, which is a block diagonal matrix. Each block is an N × N matrix, where N is the number of nodes per layer. Each diagonal block corresponds to the adjacency matrix of a layer (intra-layer connectivity, and therefore an entry a_ij = 1 in such a block corresponds to an interlayer link), while the off-diagonal blocks, given the definition of a multiplex network, are just N × N identity matrices, representing the inter-layer links between replica nodes across layers. Note that a directed multiplex, is a multiplex where at least one of the diagonal blocks is not symmetric.

Let Q be the corresponding Laplacian matrix, defined as Q = D − A where $D=\mathrm{diag}\left({d}_{i}\right)$ and ${d}_{i}={\sum }_{i=1}^{N}{a}_{ij}$ denoting the out-going degree of a node i. We consider without loss of generality a multiplex network consisting of two directed layers and interconnected via links of weight p ⩾ 0, the corresponding Laplacian matrix Q can be written in a block form as

$$\begin{equation}Q=\left[\begin{matrix}\hfill {Q}_{1}+pI\hfill & \hfill -pI\hfill \\ \hfill -pI\hfill & \hfill {Q}_{2}+pI\hfill \end{matrix}\right]\end{equation} \tag{ 1 }$$

where Q_i = D_i − A_i , i = 1, 2, denotes the Laplacian matrix for the directed graph in each layer and pI encodes the weighted interconnections across layers. The above defined Laplacian matrix for multilayer networks with directed layers is non-symmetrical, which might lead to complex rather than real spectrum. As the Laplacian satisfies a zero row sum, i.e. Qu = 0 with u denoting the all one vector, value 0 is an eigenvalue with right eigenvector u and left eigenvector denoted as y. For a nonsymmetric matrix, left and right eigenvectors are in general not the same.

To analyze the effect of network directionality on dynamical processes, we consider the standard diffusive dynamics on multilayer networks with directed layers and determine the convergence rate to a fixed point solution. Let x_i (t) denotes the state of a node i at time t. Changes of the state x_i (t) with respect to time is governed by ${\dot {x}}_{i}\left(t\right)=-{d}_{i}{x}_{i}\left(t\right)+{\sum }_{i}{a}_{ij}{x}_{i}\left(t\right)$. The diffusive dynamics on a multilayer network can be written in a matrix form as

$$\begin{equation}\frac{\mathrm{d}x{\left(t\right)}^{\mathrm{T}}}{\mathrm{d}t}=-x{\left(t\right)}^{\mathrm{T}}Q\end{equation} \tag{ 2 }$$

where ()^T denotes transposition, and x(t) is the vector encoding the state of all nodes at time t. The solution of the above governing equation follows x(t)^T = x(0)^Te^−Qt . The fixed point x^*T at which ${\dot {x}}^{{\ast}\mathrm{T}}=0$ is calculated employing the conservation law x(0)^T u = x(∞)^T u, where all one vector u is the right eigenvector corresponding to eigenvalue 0 of the Laplacian matrix Q. Let y denote the left eigenvector of Q corresponding to eigenvalue 0. Combining the conservation law x(0)^T u = x^*T u and y^T Q = 0 yields the solution of fixed point

$$\begin{equation}{x}^{{\ast}}=\left(\frac{x{\left(0\right)}^{\mathrm{T}}u}{{y}^{\mathrm{T}}u}\right)y.\end{equation} \tag{ 3 }$$

We show in appendix A a convergence to the fixed point in equation (3) is guaranteed for a connected directed graph. The convergence rate to the fixed point x* is characterized by the second smallest real part [16], denoted as Re λ₂(Q), of the eigenspectrum of the Laplacian matrix.

2.1. General analysis for diffusive behavior

To analyze the convergence behavior of diffusive dynamics on multiplex networks with directed layers, we study the spectrum of the governing Laplacian matrix and compare with the integrated system. The characteristic polynomial for the Laplacian matrix Q in equation (1) can be calculated by, applying the Schur complement theorem on the block Laplacian matrix,

$$\begin{equation}\enspace \mathrm{det}\left(\lambda I-Q\right)=\mathrm{det}\left(\left(\lambda I-\frac{{Q}_{1}+{Q}_{2}}{2}\right)\left(\left(\lambda -2p\right)I-\frac{{Q}_{1}+{Q}_{2}}{2}\right)-{\Delta}Q\right)\end{equation} \tag{ 4 }$$

where $4{\Delta}Q={\left({Q}_{1}-{Q}_{2}\right)}^{2}-2\left({Q}_{1}{Q}_{2}-{Q}_{2}{Q}_{1}\right)$. The eigen-pair of eigenvalue 2p and eigenvector $\left(u,-u\right)$ holds for directed multilayer networks, which is previously found in undirected multilayered networks [12, 21]. The joint effect of the Laplacian matrices Q₁ and Q₂ for each layer, encoded in the matrix ΔQ, determines the deviation of the convergence rate of the whole system from that of the integrated multilayer.

Depending on relations between the diffusive rate of the system as a whole and that of the integrated system, we refer two complementary regimes, superdiffusive and non-superdiffusive, as follows: (i) if Re λ₂(Q) of the overall system is greater than λ₂((Q₁ + Q₂)/2) of the integrated system, we refer to as superdiffusive, (ii) otherwise as non-superdiffusive.

$$\begin{equation}\text{Re}\enspace {\lambda }_{2}\left(Q\right)\begin{cases}{\lambda }_{2}\left(\left({Q}_{1}+{Q}_{2}\right)/2\right)\hfill & \text{Superdiffusive}\hfill \\ {\leqslant}{\lambda }_{2}\left(\left({Q}_{1}+{Q}_{2}\right)/2\right)\hfill & \text{Non}-\text{superdiffusive}\hfill \end{cases}\end{equation} \tag{ 5 }$$

We mention that the superdiffusive regime defined in this work slightly differs from the definition in literature [15, 16], where superdiffusion implies a faster diffusion on the overall system than the diffusion on the fastest layer. Here, we mainly focus on diffusive behavior comparing the overall system and the correspondent integrated system, given the later serves as the asymptotic diffusion as intercoupling strength goes to infinity and, most importantly, serves as a theoretical upper limit for diffusion on undirected multiplex networks [as shown in the following equation (8)].

Under the special case of Q₁ = Q₂, solving equation (4) yields $\text{Re}\enspace {\lambda }_{2}\left(Q\right)=\mathrm{min}\left({\lambda }_{2}\left(\left({Q}_{1}+{Q}_{2}\right)/2\right),2p\right)$ and therefore only the non-superdiffusive regime is exhibited. General cases of Q₁ ≠ Q₂ result in possibilities of both non-superdiffusive and superdiffusive regimes. Figure 1 exemplifies three typologies including (a) coupled directed and undirected layers, (b) coupled directed layers with the opposite direction, and (c) coupled directed layers with the same direction. Figure 1(a) and (c) show only non-superdiffusive regimes and figure 1(b) exhibits superdiffusion. Figures 2(a) and (c) show the non-superdiffusive and superdiffusive behavior on real-world multiplex networks, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the rich dynamical behavior of diffusion in directed multiplex, it is of paramount importance to identify the key structural properties underpinning each of the observed regimes. In this study, we mainly focus on identifying the underlying conditions driving the non-superdiffusive and superdiffusive regimes employing the principal submatrix approach and the Cauchy interlacing theorem [22] (or the Poincare separation theorem). In addition, we employ the eigenvalue property of a normal matrix, which satisfies XX^T = X^T X. For a normal matrix, the real parts of eigenvalues are the eigenvalues of the Hermitian part Re X = (X + X^T)/2 of the matrix X.

The Laplacian matrix is similar to the following matrix $\tilde {Q}$ in a transformed basis

$$\begin{equation}\tilde {Q}=\frac{1}{\sqrt{2}}\left[\begin{matrix}\hfill I\hfill & \hfill I\hfill \\ \hfill -I\hfill & \hfill I\hfill \end{matrix}\right]Q\left[\begin{matrix}\hfill I\hfill & \hfill -I\hfill \\ \hfill I\hfill & \hfill I\hfill \end{matrix}\right]\frac{1}{\sqrt{2}}\end{equation} \tag{ 6 }$$

which can be simplified as

$$\begin{equation}\tilde {Q}=\left[\begin{matrix}\hfill \frac{{Q}_{1}+{Q}_{2}}{2}\hfill & \hfill \frac{{Q}_{2}-{Q}_{1}}{2}\hfill \\ \hfill \frac{{Q}_{2}-{Q}_{1}}{2}\hfill & \hfill \frac{{Q}_{1}+{Q}_{2}}{2}+2pI\hfill \end{matrix}\right].\end{equation} \tag{ 7 }$$

The integrated Laplacian $\frac{{Q}_{1}+{Q}_{2}}{2}$ appears as a principal submatrix of the transformed Laplacian $\tilde {Q}$ of the whole system. In addition, eigenvalues of the matrix $\tilde {Q}$ are the same with eigenvalues of Q due to the similarity relation between Q and $\tilde {Q}$. Applying the Cauchy interlacing theorem on the transformed matrix $\tilde {Q}$ and the integrated system $\frac{{Q}_{1}+{Q}_{2}}{2}$ as a principal submatrix, undirected multilayer networks always satisfy the following interlacing relation

$$\begin{equation}{\lambda }_{2}\left(Q\right){\leqslant}{\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right).\end{equation} \tag{ 8 }$$

Hence, there never occurs a superdiffusion in undirected multilayer networks, regardless of the strength of intercoupling. Equation (8) serves as a theoretical limit and implies that although multilayer structure of undirected layers enhances the less diffusive layer, the diffusive rate is upper bounded by the integrated system.

3.1. Conditions for the non-superdiffusive regime

For multilayer networks with directed layers, the conditions for non-superdiffusive and superdiffusive regimes are translated into conditions when the Cauchy interlacing theorem between a non-symmetric matrix and its principle submatrix is satisfied. Starting from the special case of coupled identical layers, it follows

$$\begin{equation}\text{Re}\enspace {\lambda }_{2}\left(Q\right)={\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)\end{equation} \tag{ 9 }$$

which implies a non-superdiffusive regime for directed multilayer networks with identical directed layers for $p{\geqslant}{\lambda }_{2}\left(\left({Q}_{1}+{Q}_{2}\right)/2\right)/2$. When the transformed Laplacian $\tilde {Q}$ and all the principal matrices including the integrated matrix (Q₂ + Q₁)/2 are normal, the interlacing is also satisfied. Because for a principally normal matrix, all eigenvalues of matrices Q and (Q₂ + Q₁)/2 lie on the same complex line with an angle θ in the complex plane [23, 24]. For a principally normal matrix, the partial order of generalized interlacing on a complex plain is obtained:

$$\begin{equation}\vert {\lambda }_{2}\left(Q\right)\vert {\leqslant}\left\vert {\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)\right\vert .\end{equation} \tag{ 10 }$$

As the real part of eigenvalue reads $\vert {\lambda }_{2}\left(Q\right)\vert {\text{e}}^{\text{i}\theta }$ and θ is the same for all eigenvalues, it implies an interlacing between real parts of eigenvalues and accordingly a non-superdiffusive regime for multilayer networks with principally normal Laplacian matrices.

When the requirement of normality of all principal submatrices is relaxed to the normality of a single principal submatrix, generalized interlacing of lexicographic order can be applied to $\tilde {Q}$ and (Q₂ + Q₁)/2. For a normal Laplacian matrix and a normal principal submatrix $\frac{{Q}_{1}+{Q}_{2}}{2}$, we have that

$$\begin{equation}{\lambda }_{2}\left(Q\right){{\leqslant}}_{\theta }{\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)\end{equation} \tag{ 11 }$$

where ⩽_θ represents the lexicographic order [25] characterized by the positive cone H := {a + bi: a > 0 or a = 0 and b > 0}. Equation (11) means that ${\text{e}}^{-\text{i}\theta }\enspace {\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)-{\text{e}}^{-\text{i}\theta }\enspace {\lambda }_{2}\left(Q\right)$ lies in a positive cone, which implies ${\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right){\geqslant}\text{Re}\enspace {\lambda }_{2}\left(Q\right)$. Therefore, directed multilayer networks with normal Laplacian matrix and normal principal submatrix $\frac{{Q}_{1}+{Q}_{2}}{2}$ exhibit only the non-superdiffusive regime.

3.2. Conditions for the superdiffusive regime

Diffusion on undirected multilayer never exceeds the diffusion of integrated part. However, certain multilayer networks with directed layers break the constraint and exhibit a unique superdiffusion. In this subsection, we analyze conditions for the occurrence of such superdiffusion. For the generalized interlacing theorem to hold, it requires both the whole system and the integrated system to be normal or simultaneously close to normal. When the condition is relaxed such that the integrated system is normal or close to normal, while the whole system deviates from normality, it might occur with superdiffusion. We provide analytical evidences by firstly relating the real part of eigenvalues and eigenvalues of the Hermitian part of the original matrix, which follows

$$\begin{equation}\sum\limits _{j=N-k}^{N}\text{Re}\enspace {\lambda }_{j}=\sum\limits _{j=N-k}^{N}{x}_{j}^{{\ast}}\text{Re}\enspace \tilde {Q}{x}_{j}{\leqslant}\sum\limits _{j=N-k}^{N}{\lambda }_{j}\left(\text{Re}\enspace \tilde {Q}\right)\end{equation} \tag{ 12 }$$

where $\text{Re}\enspace \tilde {Q}$ denotes the Hermitian part of matrix $\tilde {Q}$. Because ${\sum }_{j=1}^{N}\text{Re}\enspace {\lambda }_{j}={\sum }_{j=1}^{N}{\lambda }_{j}\left(\text{Re}\enspace \tilde {Q}\right)$ and ${\lambda }_{1}\left(\text{Re}\enspace \tilde {Q}\right)=0$ for a normal subgraph, we have that $\text{Re}\enspace {\lambda }_{2}\left(Q\right){\geqslant}{\lambda }_{2}\left(\text{Re}\enspace \tilde {Q}\right)$. The equality, ${\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)={\lambda }_{2}\left(\text{Re}\enspace \frac{{Q}_{1}+{Q}_{2}}{2}\right)$, is achieved if the integrated system is normal.

Secondly, we analyze the superdiffusion condition by relating eigenvalues of the Hermitian part of the integrated system and Hermitian part of the whole system. On one hand, observing that the Hermitian part of the integrated system coincides with the principal submatrix of the Hermitian part of the whole system, we have that ${\lambda }_{2}\left(\text{Re}\enspace \tilde {Q}\right){\leqslant}{\lambda }_{2}\left(\text{Re}\enspace \frac{{Q}_{1}+{Q}_{2}}{2}\right)$ due to the Cauchy interlacing theorem for Hermitian matrices. On the other hand, ${\lambda }_{2}\left(\text{Re}\enspace \tilde {Q}\right){\geqslant}{\lambda }_{2}\left(\text{Re}\enspace \frac{{Q}_{1}+{Q}_{2}}{2}\right)+{\lambda }_{2}\left(\text{Re}\enspace \frac{{Q}_{2}-{Q}_{1}}{2}\right)$ due to Ky Fan majorization [26] theorem $\lambda \left(A+B\right)\prec \lambda \left(A\right)+\lambda \left(B\right)$. Therefore, if the difference ${\Vert}\text{Re}\left({Q}_{2}-{Q}_{1}\right){{\Vert}}_{2}$ is bounded, it leads to ${\lambda }_{2}\left(\text{Re}\enspace \tilde {Q}\right)\approx {\lambda }_{2}\left(\text{Re}\enspace \frac{{Q}_{1}+{Q}_{2}}{2}\right)$.

A superdiffusion is thus established

$$\begin{equation}\text{Re}\enspace {\lambda }_{2}\left(Q\right){\gtrapprox}{\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)\end{equation} \tag{ 13 }$$

if (i) the integrated system is or close to normal and (ii) the structure difference of coupled layers Re Q₂ − Re Q₁ is negligibly small.

The equality of ${\lambda }_{2}\left(\text{Re}\enspace \tilde {Q}\right)={\lambda }_{2}\left(\text{Re}\enspace \frac{{Q}_{1}+{Q}_{2}}{2}\right)$ is reached for a multilayer network consisting of a directed graph coupled with its reversed graph, ${Q}_{2}={Q}_{1}^{\mathrm{T}}$. In this case, it follows Re Q₂ = Re Q₁ and there exists a superdiffusion of $\text{Re}\enspace {\lambda }_{2}\left(Q\right){ >}{\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)$.

It was reported that the superdiffusion occurs when a fully connected network coupled with a network with a certain level of directionality, quantified by a metric called network directionality index [16]. Here, we prove in appendix B an open problem of the normative property of the network directionality index (NDI). In addition, we show that the condition of sufficient directionality as quantified by NDI is in line with the derived normality conditions in this study (as shown in figure 3). The derived condition for superdiffusion in this work, i.e. a close to normal and a deviation from normal (upper left panel in figure 3) for the integrated and overall multiplex, respectively, is equivalent to the condition of a sufficient level of NDI. The nonnormality condition for the non-superdiffusive regime (lower right panel in figure 3) also agrees with the NDI condition. However, the nornormality condition is widely applicable for multiplex networks with any number of directed layers, compared to a single directed layer dealt by NDI condition.

Zoom In Zoom Out Reset image size

Though the results are demonstrated for two-layer multiplex networks, numerical results on multiplex networks with a larger number of layers show analogous superdiffusive and non-superdiffusive behaviors, also depending on the nonnormality and NDI conditions of the structure of each layer. However, in those cases it is far more challenging to provide a complete mathematical proof due to the difficulty in isolating the Laplacian of each layer from the block Laplacian.

3.3. Formulation of multilayer networks with superdiffusion

To interpret the normality conditions for superdiffusion, we propose a model to construct multiplayer networks exhibiting superdiffusion with a tunable level of magnitude. When the coupled layers have identical Hermitian part of Laplacian, the minimization of the normality level of the integrated structure is translated into

$$\begin{equation}\enspace \mathrm{min}\enspace {\Vert}{Q}_{1}{Q}_{2}^{\mathrm{T}}-{Q}_{2}^{\mathrm{T}}{Q}_{1}{{\Vert}}_{2}\end{equation} \tag{ 14 }$$

with proof in appendix C. Achieving superdiffusion by minimizing equation (14) is in general dauntingly difficult, even for undirected graphs of symmetrical Laplacian matrices. The simultaneous diagonalization, which always leads to the commutativity of matrices Q₁ and ${Q}_{2}^{\mathrm{T}}$, i.e. ${Q}_{1}{Q}_{2}^{\mathrm{T}}={Q}_{2}^{\mathrm{T}}{Q}_{1}$, was posted as an open problem by Hiriart–Urruty [27] more than a decade ago. Extensive research is performed due to the close association with quadratically constrained quadratic programming and certain advances are made for real symmetrical matrices [28]. However, results on nonsymmetrical matrices are extremely scarce. Nonetheless, we propose a construction model to generate a directed multilayer network by constructing one matrix (e.g. Q₂) from a given matrix (e.g. Q₁), such that superdiffusion occurs.

We decompose the construction model into two parts: (i) constructing a graph G₂ from the replication of graph G₁ but after reversing the direction of directed links; (ii) preserving the out-degree of each node in graph G₁ such that graph Laplacian Q₂ has the same diagonal elements as Q₁. To ensure the second smallest real part of the Laplacian eigenvalue to be nonzero, we assume a single requirement of graph G₁ to be a strongly connected directed graph. To this end, we propose to reverse links in a cycle fashion, i.e. simultaneously reversing the direction of all links involved in a cycle. In addition, the number of links reversed in a cycle fashion is closely associated with the magnitude $\text{Re}\enspace {\lambda }_{2}\left(Q\right)-{\lambda }_{2}\left(\frac{{Q}_{1}+{Q}_{2}}{2}\right)$ of superdiffusion. For two directed graphs G₁ and G₂, we define a variable q as direction overlap, the fraction of links having the same direction and calculated by $q=\frac{\vert \mathcal{E}\left({G}_{1}\right)\cap \mathcal{E}\left({G}_{2}\right)\vert }{\vert \mathcal{E}\left({G}_{1}\right)\cup \mathcal{E}\left({G}_{2}\right)\vert }$ or in a matrix form as

$$\begin{equation}q=\frac{{\sum }_{i,j}{\left({A}_{1}\right)}_{ij}{\left({A}_{2}\right)}_{ij}}{{\sum }_{i,j}{\left({A}_{1}\right)}_{ij}+{\sum }_{i,j}{\left({A}_{2}\right)}_{ij}-{\sum }_{i,j}{\left({A}_{1}\right)}_{ij}{\left({A}_{2}\right)}_{ij}}.\end{equation} \tag{ 15 }$$

The denominator normalizes the direction overlap, in which q = 0 corresponds to ${Q}_{2}={Q}_{1}^{\mathrm{T}}$ and q = 1 corresponds to Q₁ = Q₂. The magnitude of superdiffusion can thus be tuned by reversing different fraction of links, parameterized by 1 − p.

Figure 4 shows diffusive behaviors on constructed multiplex networks consisting of (a) a circulant directed graph and the duplicate digraph with 1 − q reversed links and (b) a real-world genetic and protein interactions network of the Epstein–Barr virus [29, 30] and the duplicate digraph with 1 − q reversed links. The coupled digraph and its duplicate without link reversing (q = 1) displays a non-superdiffusive behavior. For 0 ⩽ q < 1, the constructed multiplex built upon both synthetic and real-world directed graph shows superdiffusion regimes. In addition, a higher fraction of reversed links (a smaller q) results in a higher magnitude of superdiffusion.

Zoom In Zoom Out Reset image size

In this paper, we study the diffusive dynamics on multiplex networks consisting of at least one directed layers. Though multilayer structure of undirected layers enhances diffusion, it is restricted to a theoretical limit by the aggregation of each layer. We show that ubiquitous property of directionality in each layer could break this limit and achieve a unique superdiffusion regime, wherein diffusion is even faster than the corresponding aggregated system. We analytically and numerically uncover that the unique superdiffusion is driven by the directionality and the underlying structure of layers rather than how strongly layers are coupled. In particular, diffusive behaviors are associated with the non-normality level of the integrated and the overall system: when both the integrated and the overall system are normal, it is assured that the multiplex exhibits only a non-superdiffusive regime; when the integrated system is normal or close to normal while the whole system is deviated from normal, it exhibits the unique superdiffusive regime. Additionally, we provide a model to construct multiplex networks, achieving superdiffusion with a tunable level of magnitude. We show that directionality induces a unique superdiffusion, a regime not existed in the counterpart of undirected multilayer, and suggest the key role of network directionality in analyzing diffusive dynamics on real-world systems which oftentimes are structurally directed and multilayered.

XW acknowledges the project 62003156 supported by NSFC and project ‛PCL Future Greater-Bay Area Network Facilities for Large-scale Experiments and Applications (LZC0019)'. AT acknowledges partial support from NSF Grants EAR-181190 and the UK Research and Innovation Global Challenges Research Fund Living Deltas Hub Grant NES0089261. YM acknowledges partial support from the Government of Aragón, Spain through a grant to the group FENOL (E36-20R), by MINECO and FEDER funds (Grant FIS2017-87519-P) and by Intesa Sanpaolo Innovation Center. The funders had no role in study design, data collection and analysis, or preparation of the manuscript.

All data that support the findings of this study are included within the article (and any supplementary files).

The convergence to the fixed point is characterized by ${\Delta}x{\left(t\right)}^{\mathrm{T}}={\left(x\left(t\right)-{x}^{{\ast}}\right)}^{\mathrm{T}}$. The rate of convergence is regulated by $\dot {{\Delta}x}{\left(t\right)}^{\mathrm{T}}=\frac{d{\left(x\left(t\right)-{x}^{{\ast}}\right)}^{\mathrm{T}}}{\mathrm{d}t}=-\left({\Delta}x{\left(t\right)}^{\mathrm{T}}+{x}^{{\ast}\mathrm{T}}\right)Q=-{\Delta}x{\left(t\right)}^{\mathrm{T}}Q$. The norm ∥Δx(t)^T∥ of the difference vector is bounded by

$$\begin{equation}{\Vert}{\Delta}x{\left(t\right)}^{\mathrm{T}}{\Vert}={\Vert}{\Delta}x{\left(0\right)}^{\mathrm{T}}\enspace {\mathrm{e}}^{-Qt}{\Vert}{\leqslant}{\Vert}{\Delta}x{\left(0\right)}^{\mathrm{T}}{\Vert}{\Vert}\enspace {\mathrm{e}}^{-Qt}{\Vert}.\end{equation} \tag{ A.1 }$$

If one mainly concerns with the convergence behavior for the decrease of ∥Δx(t)∥, the following holds for a bounded operator Q in a Hilbert space

$$\begin{equation}\underset{t\to \infty }{\mathrm{lim}}\enspace \frac{\mathrm{log}\enspace {\Vert}{\mathrm{e}}^{-Qt}{\Vert}}{t}=-{\lambda }_{\mathrm{min}}\end{equation} \tag{ A.2 }$$

where λ_min is the minimum eigenvalue of the Laplacian Q. For a connected directed graph, λ_min equals to zero which means a guaranteed convergence to the steady state x* calculated by equation (3).

In the appendix, we prove that the metric quantifying the level of directionality is normalized to 1. Rewrite the definition [16] of NDI as

$$\begin{equation}\text{NDI}=\frac{\langle {\Delta}{d}_{ij}\rangle }{\langle {d}_{ij}\rangle }\end{equation} \tag{ B.1 }$$

or equivalently,

$$\begin{equation}\text{NDI}=\frac{{\sum }_{i=1}^{N}{\sum }_{j{ >}i}{\Delta}{d}_{ij}}{{\sum }_{i=1}^{N}{\sum }_{j{ >}i}\frac{{d}_{i\to j}+{d}_{j\to i}}{2}}\end{equation} \tag{ B.2 }$$

where Δd_ij = |d_i→j − d_j→i |.

Theorem 1. For a directed cycle network consisting of N nodes and N directed links in a clockwise or counterclockwise direction, the NDI can be simplified as

$$\begin{equation}\text{NDI}=\frac{2{\sum }_{j=1}^{N}\vert {d}_{i\to j}-{d}_{j\to i}\vert }{{\sum }_{j=1}^{N}\left({d}_{i\to j}+{d}_{j\to i}\right)}\end{equation} \tag{ B.3 }$$

where ${\sum }_{j=1}^{N}\vert {d}_{i\to j}-{d}_{j\to i}\vert /N$ is the average path asymmetry and ${\sum }_{j=1}^{N}\left({d}_{i\to j}+{d}_{j\to i}\right)/2N$ is the average path length starting from an arbitrary node i and arriving at each of the remaining N − 1 nodes.

Proof. For a directed circulant graph of N nodes and N directed links in a clockwise direction, the shortest directed path from node i to node j and back to i forms a closed cycle of length N, regardless of choices of pairs of nodes (i ≠ j). Thus, we have that

$$\begin{equation}\sum\limits _{i=1}^{N}\sum\limits _{j{ >}i}\frac{{d}_{i\to j}+{d}_{j\to i}}{2}=\frac{N}{4}\sum\limits _{j=1}^{N}\left({d}_{i\to j}+{d}_{j\to i}\right).\end{equation} \tag{ B.4 }$$

Applying the symmetry of Δd_ij = Δd_ji into the numerator of (B.2), we have that

$$\begin{equation}\sum\limits _{i=1}^{N}\sum\limits _{j{ >}i}{\Delta}{d}_{ij}=\frac{1}{2}\sum\limits _{i=1}^{N}\sum\limits _{j=1}^{N}{\Delta}{d}_{ij}.\end{equation} \tag{ B.5 }$$

For a directed cycle consisting of N nodes and N directed links in a clockwise direction, the sum ${\sum }_{j=1}^{N}{\Delta}{d}_{ij}$ is the same for each node i and we arrive at

$$\begin{equation}\sum\limits _{i=1}^{N}\sum\limits _{j=1}^{N}{\Delta}{d}_{ij}=N\sum\limits _{j=1}^{N}{\Delta}{d}_{ij}.\end{equation} \tag{ B.6 }$$

Substituting (B.4) and (B.6) into the definition of NDI establishes theorem 1.(figure B1) □

Zoom In Zoom Out Reset image size

Theorem 2. Given a directed cycle consisting of N nodes and N directed links in a clockwise direction, adding a single or multiple directed links between any pair of nodes i and j generates shortcuts either in forward (i → j) or backward (j → i) direction between that pair of nodes and thus forms k (k ⩾ 1) sub-cycles with each consisting of c_k (c_k ⩽ N) nodes and c_k links. The NDI' after adding directed links is smaller than or equal to the NDI of the original directed cycle.

Proof. Adding a directed link connecting node j to i forms a sub-cycle k consisting of c_k nodes and c_k directed links in a clockwise direction, as shown in figure B2. For a pair of nodes i and j in the sub-cycle k, the difference of path asymmetry Δd_ij before and after adding directed links reads

$$\begin{equation}{\Delta}{d}_{ij}^{\prime }-{\Delta}{d}_{ij}=N-2\enspace \mathrm{min}\left({d}_{i\to j},{d}_{j\to i}\right)-\left({c}_{k}-2\enspace \mathrm{min}\left({d}_{i\to j},{d}_{j\to i}\right)\right)=N-{c}_{k}\end{equation} \tag{ B.7 }$$

and the difference of d_i→j + d_j→i before and after link addition reads

$$\begin{equation}{d}_{i\to j}^{\prime }+{d}_{j\to i}^{\prime }-\left({d}_{i\to j}+{d}_{j\to i}\right)=N-{c}_{k}.\end{equation} \tag{ B.8 }$$

Summing up all the node pairs starting from node i to the remaining N − 1 nodes yields

$$\begin{equation}{\text{NDI}}^{\prime }=\frac{2{\sum }_{j=1}^{N}\left(\vert {d}_{i\to j}-{d}_{j\to i}\vert -\left(N-{\left(m\right)}_{j}\right)\right)}{{\sum }_{j=1}^{N}\left(\left({d}_{i\to j}+{d}_{j\to i}\right)-\left(N-{\left(m\right)}_{j}\right)\right)}.\end{equation} \tag{ B.9 }$$

where m is a column vector of length N with element ${\left(m\right)}_{j}={c}_{k}$ if node j belongs to a sub-cycle of c_k nodes and c_k directed links in a clockwise direction, which can be further simplified as

$$\begin{equation}{\text{NDI}}^{\prime }=\frac{2{\sum }_{j=1}^{N}\vert {d}_{i\to j}-{d}_{j\to i}\vert -{\sum }_{j}\left(N-{\left(m\right)}_{j}\right)}{{\sum }_{j=1}^{N}\left({d}_{i\to j}+{d}_{j\to i}\right)-{\sum }_{j}\left(N-{\left(m\right)}_{j}\right)}.\end{equation} \tag{ B.10 }$$

It has been shown [16] that a directed cycle of N nodes has NDI < 1. Subtracting the same number ${\sum }_{j}\left(N-{\left(m\right)}_{j}\right)$ both in the numerator and denominator of (B.2) does not increase the NDI of the original cycle consisting of N nodes. □

Since a directed cycle of N nodes has NDI ⩽ 1 (theorem 1) and adding directed links does not increase the NDI (theorem 2), we prove that NDI ⩽ 1.

Zoom In Zoom Out Reset image size

Given the equality of the Hermitian part of the coupled directed layers, we show that the close to normal condition of the integrated system is translated to the minimization of equation (14) in the main text. The definition of normality of the integrated system reads $\left({Q}_{1}+{Q}_{2}\right){\left({Q}_{1}+{Q}_{2}\right)}^{\mathrm{T}}={\left({Q}_{1}+{Q}_{2}\right)}^{\mathrm{T}}\left({Q}_{1}+{Q}_{2}\right)$, which can be rewritten as

$$\begin{equation}{Q}_{1}{Q}_{1}^{\mathrm{T}}-{Q}_{1}^{\mathrm{T}}{Q}_{1}+{Q}_{2}{Q}_{2}^{\mathrm{T}}-{Q}_{2}^{\mathrm{T}}{Q}_{2}={Q}_{2}^{\mathrm{T}}{Q}_{1}-{Q}_{1}{Q}_{2}^{\mathrm{T}}+{Q}_{1}^{\mathrm{T}}{Q}_{2}-{Q}_{2}{Q}_{1}^{\mathrm{T}}.\end{equation} \tag{ C.1 }$$

Left multiplying Q₁ of the equality of Re Q₁ = Re Q₂ yields

$$\begin{equation}{Q}_{1}{Q}_{1}+{Q}_{1}{Q}_{1}^{\mathrm{T}}={Q}_{1}{Q}_{2}+{Q}_{1}{Q}_{2}^{\mathrm{T}}.\end{equation} \tag{ C.2 }$$

Additionally, right multiplying Q₁ yields

$$\begin{equation}{Q}_{1}{Q}_{1}+{Q}_{1}^{\mathrm{T}}{Q}_{1}={Q}_{2}{Q}_{1}+{Q}_{2}^{\mathrm{T}}{Q}_{1}.\end{equation} \tag{ C.3 }$$

Subtracting equation (C.3) from equation (C.2) results in

$$\begin{equation}{Q}_{1}{Q}_{1}^{\mathrm{T}}-{Q}_{1}^{\mathrm{T}}{Q}_{1}={Q}_{1}{Q}_{2}-{Q}_{2}{Q}_{1}+{Q}_{1}{Q}_{2}^{\mathrm{T}}-{Q}_{2}^{\mathrm{T}}{Q}_{1}.\end{equation} \tag{ C.4 }$$

Analogously, left and right multiplying Q₂ establishes

$$\begin{equation}{Q}_{2}{Q}_{2}^{\mathrm{T}}-{Q}_{2}^{\mathrm{T}}{Q}_{2}={Q}_{2}{Q}_{1}-{Q}_{1}{Q}_{2}+{Q}_{2}{Q}_{1}^{\mathrm{T}}-{Q}_{1}^{\mathrm{T}}{Q}_{2}.\end{equation} \tag{ C.5 }$$

Substituting the equations (C.4)–(C.5) into the definition of normality leads to

$$\begin{equation}{Q}_{1}{Q}_{2}^{\mathrm{T}}-{Q}_{2}^{\mathrm{T}}{Q}_{1}+{Q}_{2}{Q}_{1}^{\mathrm{T}}-{Q}_{1}^{\mathrm{T}}{Q}_{2}=0,\end{equation} \tag{ C.6 }$$

from which we have that ${Q}_{1}{Q}_{2}^{\mathrm{T}}-{Q}_{2}^{\mathrm{T}}{Q}_{1}$ is a skew symmetric matrix and the real part of all eigenvalues are 0. Minimizing the normality level of the integrated system can be quantified and analyzed by the minimization of equation (14).