Cross-diffusion on multiplex networks

During the past decades, pattern formulation with reaction–diffusion systems has attracted great research interest. Complex networks, from single-layer networks to more complicated multiplex networks, have made great contribution to the development of this area, especially with emergence of Turing patterns. While among vast majority of existing works on multiplex networks, they only take into account the simple case with ordinary diffusion, which is termed as self-diffusion. However, cross-diffusion, as a significant phenomenon, reveals the direction of species’ movement, and is widely found in chemical, biological and physical systems. Therefore, we study the pattern formulation on multiplex networks with the presence of both self-diffusion and cross-diffusion. Of particular interest, heterogeneous patterns with abundant characteristics are generated, which cannot arise in other systems. Through linear analysis, we theoretically derive the Turing instabilities region. Besides, our numerical experiments also generate diverse patterns, which verify the theoretical prediction in our work and show the impact of cross-diffusion on pattern formulation on multiplex networks.


Introduction
Dynamical processes coupled in space can spontaneously generate macroscopic phenomena, such as Turing patterns, via the microcosmic interactions between the species in systems. This is widely observed in nature, and considered as a crucial frontier in science nowadays. It also has comprehensive potential applications in many fields, such as, but not limited to, biology, chemistry and physics [1]. Therefore, the investigation of patterns, arising from dynamical processes, is of great importance, and at a cutting edge of academic researches. Usually, dynamical processes are modeled by reaction-diffusion (RD) systems, which have been widely employed over decades to characterize the evolution in dynamical systems, in which local elements diffuse and interact with a certain way. Many efforts [2][3][4][5][6][7][8][9][10][11] have been carried out to analyze heterogeneous RD processes, since Turing [12] showed that the formation of a periodic spatial pattern can be inspired by perturbation around the equilibrium.
Furthermore, the space can also be characterized as discrete media, i.e. network, where nodes represent the residences of the species and links act as paths of movement [13]. Othmer and Scriven [14] first brought networks into instability analysis, under the situation of biological reaction in cell. They have put forward a universal theoretical structure for pattern formulation in networks, which is subsequently extended by plenty of scholars [15][16][17][18][19][20][21][22]. In 2010, It is remarkable that Nakao and Mikhailov [23] have filled the gap of the instability analysis in terms of large-scale networks. They found that the dynamics of reaction-diffusion systems organized in large network had obvious differences from those coupled in small networks, yielding distinct non-uniform patterns. Inspired by these advances, great efforts have been made in the exploration of such network instability [24][25][26][27][28][29][30][31]. Moreover, Turing instabilities have also been explored in more complicated systems, such as [32,33]. In reference [33], Asllani et al carry out the analysis on multiplex networks by means of a perturbative approach, resulting in self-organized patterns, which are instead prevented in monolayer networks. Kouvaris et al [32] propose an approximation mathematical framework for inducing Turing patterns in networks consisting of two layers, proving that topology of multiplex networks could inspire heterogeneous instabilities in RD systems.
Nevertheless, the exploration of RD systems organized on multiplex networks is restricted to ordinary diffusion, which is termed as self-diffusion. It is a very common scene that local species have the tendency to diffuse close to or away from the high concentrations of other kinds of species, which is, apart from self-diffusion, called cross-diffusion. This kind of diffusion has been subsequently studied [34][35][36][37][38][39][40] since Shigesada et al [41] brought the concept 'cross-diffusion' to the world in 1979. Despite the fact that much attention has been taken to cross-diffusion, but almost no one performs the analysis of cross-diffusion which happens on multiplex networks. Hence, there is still a blank area for the exploration of cross-diffusion on multiplex networks, which is exactly the work we have done in this paper. In our framework, different species reside in different layers; they diffuse in their own layers via the inner-layer links, but react across the layer through the inter-layer links; the diffusion direction of activators is influenced by the inhibitors' density, i.e. cross-diffusion, which is detailed explained in section 2. The patterns observed in our scheme have striking differences from traditional frameworks. Moreover, the experiment results are in accordance with the deduction derived from theoretical analysis. In general, our work presents a vital extension to the modeling of RD systems and a framework to analyze such systems mathematically, which will shed light on the exploration of dynamical processes in nature.
The rest of this paper is arranged as follows: section 2 illustrates the model of cross-diffusion on multiplex networks proposed in our work; section 3 conducts the Turing instability analysis; experiment results including patterns on multiplex networks with different compositions are shown in section 4; in the end, the whole paper comes to a conclusion.

Cross-diffusion model on multiplex networks
For a better comprehending of the advances in our framework, this section is started with a review of RD models with ordinary diffusion organized in discrete space. A two species RD system on monolayer networks can be cast by a set of two ordinary differential equations, which is shown as follow: where u and v indicate the activators and inhibitors, u i and v i correspond to the local concentrations of these reactors who occupy the node i. Generally, the reaction dynamics is specified by f (u i , v i ) and g(u i , v i ), d 11 and d 12 represent the self-diffusion coefficients. L ij = A ij − k i δ ij , where L ij is the element of network Laplacian matrix L. δ ij illustrates the Kronecker product, whose value will be equal to 1, if and only if i = j. Otherwise, δ ij = 0. A denotes the adjacent matrix, the element A ij of which equals to 1 if there exists a link between node i and node j, if not, A ij = 0. In the aforementioned model, different species reside in the same network where the reaction takes place in the nodes, diffusion is conducted through the links. To extend the model showed by equation (1), Asllani et al studied the Brusselator model as one classical RD system on multiplex networks [33], where they found that the diffusion between layers can seed instabilities, leading to diverse patterns which are obstructed with the limit of uncoupled networks. Soon after, A qualitative research by Kouvaris et al put forward a new framework for pattern formulation on multiplex networks [32], which could be summarized as follow: where L (u) ij denotes the Laplacian of activators' layer, similarly, L (v) ij indicates the Laplacian of inhibitors' layer. The meanings of other notations are the same with model (1). In model (2), different particles occupy distinct nodes in separate layers, they react across the layers within the node pairs and diffuse in their own layers as shown in figure 1.
Compared to self-diffusion, cross-diffusion has also attracted much attention over the past years. However, to our knowledge, the explorations of cross-diffusion on multiplex networks are rare. Hence, we introduce a framework of RD systems with cross-diffusion organized on multiplex networks, which is described by the equations: where d 12 represents the cross-diffusion coefficient. When d 12 > 0, the activators travel towards the direction of lower concentration of inhibitors. Whereas, when d 12 < 0, the activators will diffuse along the direction to higher concentration of inhibitors. As a representative example of RD systems, Brusselator model, which was first established by Prigogine and Lefever [42], is considered in our work. The reaction dynamics in Brusselator can be cast as follow: where a and b are constant parameters. Brusselator model, as showed in equation (4), is quite rich in dynamical behavior, which is also the reason we decide to consider it in our work.

Turing instability analysis
Small perturbations δu i and δv i are then added to the fixed stable point (a, b/a), yielding that: Substitute equations (5) and (4) to model (3), and conduct the taylor expansion on equation (3), we can get: where f u , f v , g v and g u are partial derivative of dynamics function at the stable fixed point, the specific calculations are given as: Transform equation (6) to matrix form: where ω = (δu 1 , . . . , δu N , δv 1 , . . . , δv N ) is the perturbation vector; the elements of matrix P and S in equation (8) are indicated as follows: Here, I denotes the N × N identity matrix; L (u) and L (v) are the Laplacian matrix of activators' layer and inhibitors' layer respectively. Decomposing S as S = S A − S D , where: Assuming that both layers of reactors are dense enough so that k u 1 and k v 1. It can be clearly seen that the values of elements in S A are much smaller than those in S D . Hence, simillar to reference [32], S A can be neglected in the formula S = S A − S D , yielding S ≈ −S D . Inserting this approximation and equations (10) into (8), we can obtain: ∂ω According to equation (11), the growth rate λ can be cast through: Where k u and k v represent the degrees of nodes in activators' layer and inhibitors' layer respectively. Simplifying equation (12), we can obtain the characteristic equation of growth rate λ: Solving equation (13), we can get: Where: Turing instabilities require that Re(λ) < 0 for all λ when there is no diffusion, which derives that: Following this, the range of constant parameters will be work out by substituting equation (7) into (16), deriving that: When it comes to the case with the presence of diffusion, Turing patterns will arise when there exists λ that satisfy Re(λ) > 0. Examing formula (14), it is found that Δ < 0 will meet the requirements. According to equation (15), the demand Δ < 0 is equivalent to: Therefore, equations (17) and (18)

Turing patterns on multiplex networks
In this section, the constant parameters a and b are assigned with the value of 4 and 16.9 respectively, which fits all the requirements showed by equation (17). Some other fixed parameters are set as follows: According to equation (18), the Turing boundaries of the degree combinations of node pairs are shown in figure 2. The area below the boundary is defined as Turing region (TR), the area above the boundary is named as unable Turing region (UTR). The expression 'degree combination' means the combination of degrees for each node pair, where nodes are located in separated layers yet have one-to-one correspondence. If the degree combination of a node pair is located in TR, this pair of nodes will satisfy the demand of Turing instability described by equation (18). On the contrary, if the degree combination of a node pair is located in UTR, this pair of nodes will not meet the conditions. The value of cross-diffusion coefficient d 12 varies from −2 to 2 with interval 1, which is reflected by different curves in figure 2. It can be clearly seen that, the larger the value of d 12 is, the vaster the TR is. Hence, we can theoretically conclude that the increasing of the value of cross-diffusion coefficient will help generate Turing instabilities. However, when d 12 < 0, the TR does not contain positive degrees for nodes in activator layer, which is not realistic. Therefore, the situation d 12 < 0 is not considered in the following experiments.
To verify the findings aforementioned, we conduct the experiments on multiplex networks. Species react in the nodes through inter-layer links (red imaginary lines in figure 1), but diffuse in their own layers through inner-layer links (blue and green solid lines in figure 1). Both layers, Q(u) and Q(v), are of same size (number of nodes N = 3600) and scale-free but have different topology structures, corresponding to different mean degrees k u = 20, k v = 30. As a matter of fact, the topology differences between Q(u) and Q(v) are consistent with the fact that two different species usually have distinct movement routes in nature.  figure 3(m), all degree combinations are in TR which is contrary to the situation that d 12 = 0.3. Therefore, the ideal Turing instabilities (presented by figures 3(n)-(p)) will be generated, in which the density of species in any node is eventually dynamically stable but far away from the equilibrium. In short, it can be observed from the top row to the bottom row that, the value increasing of cross-diffusion coefficient d 12 has facilitative effect on Turing instability arising.  Consequently, some supplementary proofs need to be put forward. Following the analysis in figure 3 and figure 4 shows a more comprehensive study of Turing patterns in terms of amplitude. The amplitude is defined as: The amplitude evolutions of all node pairs are shown in figure 4(a), different lines correspond to different values of d 12 , of which the range is the same as the setting in figure 3. We can see that all four lines ultimately become stable. But the line corresponding to larger value of d 12 has the stable amplitude with higher value. The bigger stable amplitude means that Turing instabilities are stronger. Therefore, from the aspect of amplitude, cross-diffusion also acts as a catalyst of Turing instabilities. The bigger d 12 is, the stronger Turing instabilities are.
After the investigation of amplitude for all node pairs, there may exists doubt that is there any amplitude difference between the node pairs in TR or UTR? To address this confusion, figure 4(b) exhibits the amplitude evolution of different kind of node pairs. Red lines represent the amplitude evolutions of node pairs in UTR, while the gray lines denote the amplitude evolutions of node pairs in TR. Different shapes of lines correspond to different d 12 . Here, the values of d 12 are 0.6 and 1.5, of which the node pairs classifications are the same as those in figures 3(e) and (i) respectively. Therefore, figure 4(b) provides a more visual presentation of the difference of node pairs in TR and UTR. One can well perceives that, regardless of the value of d 12 , the stable amplitude of the node pairs in UTR is always smaller than that of the node pairs in TR. Just as said before, the bigger stable amplitude means that Turing instabilities are stronger. So it can be concluded that the node pairs in TR contribute more to the arising of Turing instabilities, compared to those in UTR. In figures 3 and 4, it should be noted that node 1 in activators' layer Q(u) is always paired with node 1 in inhibitors' layer Q(v). The indexes of the nodes in Q(u) are descending sorted by degree, which means k 1 u k 2 u · · · k N u . But the indexes of nodes in Q(v) are randomly assigned.
We then explore the pattern formulations of Turing instabilities on lattice network. Except that d 12 = 3.3, the other parameters are set the same as the aforementioned figures. For simplicity, the lattice network with mean degree 3.9 is abbreviated as LA3.9, other lattice networks with different mean degrees are abbreviated in the similar way. Figure 5 shows heterogeneous Turing patterns of activators arising from the multiplex networks which consist of lattice layers. From the first row to the fourth row, the layers of activators are LA3.9, LA7.7, LA11.2 and LA21.6 respectively. From the left column to the right column, the layers of inhibitors are LA3.9, LA7.7 and LA11.2 separately. The specific color represents the corresponding density of activators according to color bar, density difference in space is reflected by color difference. We can see that with the presence of cross-diffusion, diverse patterns have been generated by adjusting the mean degree of the layers of activators and inhibitors. Those patterns include but not limited to spot patterns, strip patterns and labyrinth patterns. See figure 5 to enjoy the beautiful patterns generated by Brusselator model. In figure 5, node 1 in activators' layer Q(u) is always coupled with node 1 in inhibitors' layer Q(v). The indexes of nodes in Q(u) and Q(v) are both assigned according to their positions in lattice, which is different from the marking rules in figures 3 and 4.

Conclusion
In a recent research, Asllani et al [33] showed that the connection through the inter-layer links can seed instabilities, resulting in patterns with distinct characteristics, which are yet prevented in monolayer networks. Latter, Kouvaris et al [32] put forward a new analytical framework to study the pattern formulation on multiplex networks, where the Turing region in theory is applicable for degree combinations of nodes in separate layers.
Following the great efforts made by former scholars, we extend the research on pattern formulation on multiplex networks, by introducing cross-diffusion, which is generally existing in nature. More specifically, a novel theoretical model of multiplex reaction-diffusion systems, with the presence of both self-diffusion and also cross-diffusion, is proposed in this paper. It is very common that the movement directions are influenced by other kinds of species, the phenomenon of which is termed so called 'cross-diffusion'. This is also the reason why we bring cross-diffusion in.
In our scheme, species react between network layers through the inter-layer links but diffuse in their own layers via the inner-layer links. The movement of activators is towards the lower concentration of inhibitors, which is reflected by the positive value of d 12 in equation (3). Through the linear analysis of Turing instabilities, Turing region (TR) and unable Turing region (UTR) could be theoretically obtained through the boundary derived from equation (18). Here, similar to reference [32], degree combinations of nodes in separate layers are the determinants of whether the Turing condition is satisfied. Numerical simulations are performed to validate the results inferred from the linear analysis. We find that when the degree combinations of all node pairs are in UTR, there is no phenomenon like Turing patterns. However, if there exist node pairs whose degree combinations are in TR, Turing instabilities will arise. Furthermore, with the number of node pairs in TR increasing, the amplitude will be larger which corresponds to more apparent Turing patterns. From another point of view, when other settings in our framework are fixed, the increasing of cross-diffusion coefficient d 12 will contribute to the increasing of the node pairs in TR, as shown in figure 3, thus promoting the arising of Turing instabilities. Aim to exploring the diversity of patterns, experiments on multiplex networks, each layer of which is a 60 × 60 lattice, are then conducted. Heterogeneous patterns, such as dot patterns, catena patterns and maze patterns, have been generated with different properties. We note that the patterns arising in our paper differ from those in other frameworks.
Multiplicity of patterns is of great importance in describing and characterizing dynamical processes in nature. As is known, the intrinsic reaction mechanism of reaction-diffusion systems is crucial for the formation of specific patterns; while the integration of networked population advances the further development of this area and its applications in biology, physics and chemistry systems. As for cross-diffusion, which we introduced to multiplex networks in this paper, it greatly promotes the formulation of Turing patterns in more complicated space. Thus, the present research could shed light on various domains, ranging from biology to physics, it will also inspire further studies in future, such as the exploration of Hopf bifurcation with the presence of both self-diffusion and cross-diffusion, the control of complex patterns.