# Pattern formation in multiplex networks

## Abstract

The advances in understanding complex networks have generated increasing interest in dynamical processes occurring on them. Pattern formation in activator-inhibitor systems has been studied in networks, revealing differences from the classical continuous media. Here we study pattern formation in a new framework, namely multiplex networks. These are systems where activator and inhibitor species occupy separate nodes in different layers. Species react across layers but diffuse only within their own layer of distinct network topology. This multiplicity generates heterogeneous patterns with significant differences from those observed in single-layer networks. Remarkably, diffusion-induced instability can occur even if the two species have the same mobility rates; condition which can never destabilize single-layer networks. The instability condition is revealed using perturbation theory and expressed by a combination of degrees in the different layers. Our theory demonstrates that the existence of such topology-driven instabilities is generic in multiplex networks, providing a new mechanism of pattern formation.

Distributed active media support a variety of self-organized patterns, such as stationary and oscillatory structures, spiral waves, and turbulence turing (); murray-book (); mik-book-I (). Such media are often described by reaction-diffusion systems and consist of elements obeying an activator-inhibitor dynamics with local coupling. In his pioneering paper turing (), Turing showed that a uniform steady state can be spontaneously destabilized, leading to a spontaneous formation of a periodic spatial pattern, when reacting species diffuse with different mobilities. It was later proposed by Gierer and Meinhardt MEIN00 () that an activator-inhibitor chemical reaction is a typical example achieving Turing’s scenario. Turing instability is a classical mechanism of self-organization far from equilibrium, and plays an important role in biological morphogenesis. It has been extensively studied in biological MEIN00 (); Sick2006 (); Maini06 () and chemical Ouyang1991 () systems, as well as real ecosystems Rietkerk2008 (); Liu2013 ().

The active elements can also be coupled in more complicated ways, forming complex networks Karlsson2002 (); Bignone2001 (). Complex networks are ubiquitous in nature barrat-book-08 (); two typical examples are epidemics spreading over transportation systems COL06 () and ecological systems where distinct habitats communicate through dispersal connections Hanski1998 (); Urban2001 (); Fortuna2006 (); HOL08 (). Theoretical studies of reaction-diffusion processes on complex networks have recently attracted much attention barrat-book-08 (); Bocaletti2006 (); Arenas2008 (); KOU12 (); KOU14a (). Othmer and Scriven oth71 (); oth74 () developed the general mathematical framework to describe Turing instability in networks, and provided several examples of small regular lattices. Afterwards, Turing patterns were explored in small networks of chemical reactors hor04 (); MOO05 (). More recent work in this area includes detailed studies of Turing bifurcation and related hysteresis phenomena in large complex networks Nakao2010 (); Wolfrum2012a (), and oscillatory Turing patterns in multi-species ecological networks Hata2014 ().

In nature, the active elements of a system can communicate through different types of pathways with different architecture. Such a system with multiple types of links can be represented as a special type of complex network called a multiplex network Boccaletti2014 (). Recent theoretical studies have shown that the spectral properties of multiplex networks are significantly different from those of single-layer networks Gomez2013 (); SOL13 (); dedomenico2013 (); Radicchi2013 (); Boccaletti2014 (), and that these differences affect the diffusion processes occurring on the network Gomez2013 (); SOL13 (). Consequently, the emergent dynamics can exhibit new kinds of patterns. Examples include the breathing synchronization of cross-connected phase oscillators Louzada2013 () and the emergence of a metacritical point in epidemic networks, where diffusion of awareness is able to prevent infection and control the spreading of a disease Granel2013 (). Turing patterns have also been discussed in the context of multiplex networks ASL14 ().

It has been reported that many man-made networks and real ecosystems are spatially fragmented in such a way that different species can migrate using different paths in separate layers Stegeman2002 (); Fang2009 (); Xuan2013 (); Cardillo2013a (); Buono2014 (). In studies of classical swine fever, for example, it was found that an individual can spread the infection by different types of contacts characterized by different infection rates Stegeman2002 (). Moreover, the role of different but overlapping transportation networks was considered in a study exploring the diffusion pattern of severe acute respiratory syndrome near Beijing Fang2009 ().

This literature leads us to consider a new class of dynamical systems, multiplex reaction networks, where reacting species are transported over their own networks in distinct layers, but can react with each other across the inter-layer connections. This paper provides a general framework for multiplex reaction networks and constructs a theory for self-organized pattern formation in such networks. As a typical example, we investigate a diffusively-coupled activator-inhibitor system where Turing patterns can develop.

## I Multiplex reaction networks

We consider multiplex networks of activator and inhibitor populations, where the different species occupy separate network nodes in distinct layers. Species react across layers according to the mechanism defined by the activator-inhibitor dynamics, and diffuse to other nodes in their own layer through connecting links (see figure 1). Such a process can be described by the equations

(1a) | |||||

(1b) |

where and are the densities of activator and inhibitor species in nodes and of layers and , respectively. The superscripts and refer to activator and inhibitor. The activator nodes are labeled by indices in order of decreasing connectivity. The same index ordering is applied to the inhibitor layer. The functions and specify the activator-inhibitor dynamics. The Laplacian matrices and describe diffusion processes in the two layers, and the constants and are the corresponding mobility rates (see details in the Methods section).

As a particular example we consider the Mimura-Murray ecological model Mimura1978 () on a multiplex network consisting of two scale-free layers. In the absence of diffusive coupling, such that and , the multiplex system relaxes to a uniform state, i.e. for all . The homogeneous densities are determined by (see Methods). Under certain conditions, which we present here, non-uniform patterns can evolve from an instability driven by the multiplex structure.

## Ii Linear stability of the uniform state

In simplex networks, where , the uniform state may undergo a Turing instability as the ratio increases and exceeds a certain threshold. The instability leads to the spontaneous emergence of stationary patterns consisting of nodes with high or low densities of activators Nakao2010 (). Such diffusion-induced instability can also take place in multiplex reaction networks (1). This phenomenon can be explained through a linear stability analysis with non-uniform perturbations. We introduce small perturbations, and , to the uniform steady state, as follows: . We then substitute the perturbed state into equations (1) to obtain a set of coupled linearized differential equations. Finally, by means of an approximation technique described fully in the Methods section, we obtain a characteristic equation for the growth rate of the perturbations for each pair of nodes.

The onset of the instability occurs when for some pair of nodes and . The instability condition is fulfilled when these nodes possess a combination of degrees and such that, the equation

(2) |

is satisfied. Here, and are partial derivatives at the uniform steady state. Condition (2) implies that a sufficiently large value of brings about instability, in the same manner as the Turing instability. However, an alternative scenario of the instability is revealed by Eq. (2). This can happen by increasing , even if the mobilities are equal (). This instability occurs in a strikingly different regime from classical diffusion-induced instabilities.

Figure 2a shows the linear stability of system (1) for varying , holding fixed. We clearly see that the uniform steady state is always a solution of the multiplex system. It is linearly stable (green line) for small values of . But at some critical value of which satisfies equation (2), the system undergoes a transcritical bifurcation (red point) and becomes unstable (magenta line). Two new branches of solutions arise from the transcritical bifurcation. The unstable branch (magenta line) undergoes a second bifurcation (blue point), this time a saddle-node, giving rise to a new branch of stable solutions (green line) different from the uniform steady state. Figure 2b shows the transcritical (red line) and the saddle-node bifurcation (blue line) in the - plane. The curve of the transcritical bifurcation is given by equation (2), while the curve of the saddle-node bifurcation has been derived by numerical continuation. One can see from Eq. (2) that by increasing , the boundary curve (red line) asymptotically approaches . This indicates that the instability can be observed with sufficiently large , if the mean degree of the activator layer is less than this value. This fact reveals an important difference from the classical Turing instability, which always takes place by increasing irrespective of Nakao2010 ().

The diffusion-induced instability occurs on the transcritical bifurcation. However, non-uniform patterns can also develop after the saddle-node bifurcation. In other words, we find that multiplex systems exhibit multistability in the area between these two bifurcations (cyan), where a branch of stable solutions coexists with the uniform steady state.

## Iii Pattern formation arising from the instability

Suppose that the multiplex system starts almost in the uniform steady state with small perturbations. Equation (2) allows us to identify pairs of nodes (, ) where the small perturbations will be amplified, so that these nodes leave the uniform state, triggering the formation of a non-uniform stationary pattern. Such a pattern cannot develop from pairs of nodes possessing degrees in the grey area of Figure 2b, where only the uniform state exists. However, pairs of nodes with degrees in the yellow area, beyond the transcritical bifurcation, are unstable. Under small perturbations they can leave the uniform state, yielding the formation of a stationary non-uniform pattern. The cyan area between the two bifurcations indicates that the system exhibits multistability, where the uniform steady state coexists with a branch of solutions corresponding to non-uniform patterns.

We verify this scenario for a multiplex network where both layers, and , are scale-free. Figure 3a displays the actual degree combination (, ) for each pair of nodes , (orange points) of this network in the - plane, together with the bifurcation curves. Three pairs of nodes, the critical ones, denoted by stars, have degrees exceeding the instability threshold. Thus, a non-uniform pattern starts to grow from these nodes. The critical node denoted by the red star is the first to spontaneously leave the uniform state, as shown in figure 3b. Next, figures 3c,d show that the critical nodes denoted by the green and blue stars rapidly differentiate from the uniform state. Finally, triggered by these growing perturbations, all nodes leave the steady state to establish a non-uniform pattern (Figure 3e).

Multistability corresponding to the cyan are of Figure 2b has been studied via numerical simulations. Figure 4a shows the amplitude of the observed patterns (see Methods section), averaged over different simulations. Each point of the diagram is the average of ten different implementations of with the same mean degree ; is fixed. We clearly see that the amplitude is zero; i.e., the uniform state is the only stable attractor of the system, for smaller than a critical threshold . However, a more detailed look in the vicinity of this transition reveals that a number of different stationary patterns could be identified for the same parameter values. As an example, figure 4b shows the amplitudes in three simulations where different perturbations have been applied to the same sequence of multiplex networks. Starting from the uniform state with small perturbations, the instability occurs at some critical threshold, resulting in a small abrupt increase of amplitude. Different perturbations result in different values for the instability threshold.

Obviously, different values lead to patterns of different amplitudes. Figure 3e shows a pattern for , close to the transition. However, patterns where more nodes leave the uniform state can also develop far from the transition. Figures 5a-f show the evolution of small perturbations in the uniform state and the formation of a non-uniform pattern in a multiplex network with scale-free layers of nodes, and mean degrees and . Under the influence of small perturbations, some critical nodes differentiate rapidly from the uniform steady state. Afterwards, nonlinear effects (which are not described by our theory) drive the multiplex system to self-organize into a stationary pattern with two separate group of nodes. The separation between nodes of low and high activator densities is more pronounced in nodes with small degrees , while nodes with large tend to sustain their initial state. Figures 6a,b show this pattern in the activator and inhibitor layers respectively, whereas figure 6c shows the actual multiplex pattern.

## Iv Discussion

We have proposed a new class of dynamical systems, multiplex reaction networks, where each reacting species occupies its own network layer and reacts with the other species using cross-layer contacts. As a demonstration of this new reaction scheme, we investigate pattern formation induced by diffusive transport in a multiplex network with two reacting species. Our theory, based on linear stability analysis with perturbations around the uniform steady state, correctly predicts the instability threshold observed in numerical simulations of the multiplex network.

If the different layers have the same architecture, i.e. , then this multiplex diffusion-induced instability reduces to the well-known Turing instability which may occur when the inhibitor diffuses much faster than the activator. Our theory (2) predicts that the analogous instability can also appear in multiplex reaction networks by increasing the inhibitor diffusion rate. However, a significantly different kind of instability can occur in multiplex reaction networks, even if the two species have the same mobilities (). This new instability is related to the degree combination (, ) of a pair of nodes. The basic condition for any given pair of nodes , to undergo instability is that their degrees and satisfy Eq. (2). Indeed, the instability always takes place for any which is less or equal to the value calculated from Eq. (2), for a given large .

Similar to simplex networks, multiplex systems exhibit multistability. The onset of pattern formation can occur even before the instability described by Eq. (2). The minimal condition for developing non-uniform patterns is that in a pair of nodes , the degree is less than or equal to the value on the saddle-node bifurcation curve that corresponds to . In the multistability regime, different stationary patterns can coexist with the uniform steady state for the same parameter values.

Although the observed properties of the stationary patterns are similar to those found in simplex networks Nakao2010 (); Wolfrum2012a (), the cause of destabilization of the uniform state is different. This cause is only characteristic of multiplex networks and lies in the relationship between and for a pair of nodes. Therefore, the purposeful design of nonequilibrium patterns should be possible by tuning the architecture of the multiplex structure. Recently, new algorithms for building multiplex networks with positive or negative degree correlations across the layers have been proposed Lee2012 (); Nicosia2013b (); Kim2013 (). Using these algorithms, we can design multiplex networks where the onset of instability is controlled by tuning the degrees and , and the source of instability can be located at any desired pair of nodes , .

Multiplex networks can be used to represent different types of interaction Stegeman2002 (); Granel2013 (); Buono2014 () or different transportation lines Fang2009 (); Kaluza2010 (); Cardillo2013a () between discrete nodes. In ecological multiplex networks, for example, pairs of nodes might represent separate habitat patches which communicate through dispersal connections. However, prey and predators may use different connections (such as forest paths, rivers and tributaries or various transportation systems) to move among the fragmented habitats. Often, predators have more choices to move; in our representation their layer is more densely connected than the prey’s layer. This is exactly the sort of situation that favors the formation of instability and the subsequent establishment of non-uniform patterns. Considering that self-organized patterns can be found in real ecosystems Rietkerk2008 (); Liu2013 () it is possible that such patterns can also be observed in natural ecological systems for which the multiplex structure is innate.

## V Methods

### Layer architecture.

In the numerical simulations, each layer is a scale-free network constructed by the preferential attachment algorithm Albert2002 (). The network structure is determined by a symmetric adjacency matrix , whose elements are if there is a link connecting nodes and , and otherwise. The degree, i.e. number of links, of node is defined as . The network Laplacian matrix is given by the expression .

The activator’s network was constructed with mean degree . The same network was used throughout all numerical simulations. Each simulation uses a different realization of the inhibitor’s network , whose mean degree is varied between simulations. The superscripts and refer to activator and inhibitor. For convenience, the indices of nodes in the layer are assigned in order of decreasing degrees : that is, . The nodes in the layer follow the ordering of their counterpart in , so for example node in the inhibitor network (the most highly connected node) always corresponds to node in the activator network, but the latter may or may not be highly connected.

### Multiplex networks.

The multiplex networks used in our numerical simulations consist of two separate layers and two different types of links, intra-layer and inter-layer links. Intra-layer links are described by the adjacency matrices and limit the diffusional mobility of the species. Inter-layer links connect every node of layer to its counterpart in layer . They represent the reaction dynamics defined in the functions and .

### Activator-inhibitor dynamics.

We choose the Mimura-Murray model Mimura1978 () as an example of an activator-inhibitor system. In this model the dynamics are given by the functions and , where correspond to the densities of activator and inhibitor respectively. The chosen parameters are , yielding the linearly stable fixed point . This requires the networks to satisfy and , where is the Jacobian matrix , and , , and are partial derivatives.

### Linear stability analysis.

The linear stability analysis is performed using a perturbation method. We introduce small perturbations to the uniform steady state , as . Substituting into equations (1), we obtain the linearized differential equations and . Alternatively, the linearized differential equations can be written as , where is the perturbation vector, and ; is the identity matrix. For the linear stability analysis, the perturbation vector should be expanded over the set of eigenvectors of the matrix . It is, however, difficult to calculate them for different network topologies, i.e. different Laplacian matrices and . Here we propose an approximation technique to analyze the linear stability of the system. Matrix is split into , where and . The matrices and are the adjacency matrices of layers and , respectively. The matrices and are the corresponding degree matrices, which have the nodes degrees in the main diagonal and are zero elsewhere. Then, matrix can be rewritten as , where . Examining matrices and , the first has elements with values of order or , while the second has elements with values of order or . If both layers are dense enough that and , we can clearly see that the elements of matrix have larger values than those of matrix , so that can be neglected. This approximation yields the approximate linearized equation . The characteristic equation for the eigenvalues is then given by

and is the same for each pair of nodes , .

This approximation neglects entirely the matrix , which is associated with the precise architectures of the layers. Instead, each node is characterized only by its degree. This is quite similar to the powerful mean-field methods used for analyzing Turing patterns in single-layer networks Nakao2010 (); Wolfrum2012a (), and is always valid for multiplex networks consisting of layers with small diameters.

### Amplitude of non-uniform patterns.

The amplitude of a non-uniform pattern is quantified as .

### Acknowledgments.

This study was supported by the EU/FP7-2012-STREP-318132 in the framework project LASAGNE. Partial support from Spanish DGICYT Grant No. FIS2012-38266-C02-02, Cross-ministerial Strategic Innovation Promotion Program and JSPS KAKENHI (Grant No. 26880033) in Japan is also acknowledged.

### Author Contributions

N.E.K. and A.D.-G. conceived and designed the study. N.E.K. performed numerical simulations. All authors carried out the analysis and wrote the article.

### References

- Turing, A. M. The Chemical Basis of Morphogenesis. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237, 37–72 (1952).
- Murray, J. D. Mathematical Biology (Springer-Verlag, Berlin Heidelberg, 2003), third edn.
- Mikhailov, A. S. Foundations of Synergetics I: Distributed Active Systems (Springer-Verlag, Berlin, 1994), second edn.
- Meinhardt, H. & Gierer, A. Pattern formation by local self-activation and lateral inhibition. BioEssays : news and reviews in molecular, cellular and developmental biology 22, 753–60 (2000).
- Sick, S., Reinker, S., Timmer, J. & Schlake, T. WNT and DKK determine hair follicle spacing through a reaction-diffusion mechanism. Science (New York, N.Y.) 314, 1447–50 (2006).
- Maini, P. K., Baker, R. E. & Chuong, C.-M. Developmental biology. The Turing model comes of molecular age. Science (New York, N.Y.) 314, 1397–8 (2006).
- Ouyang, Q. & Swinney, H. L. Transition from a uniform state to hexagonal and striped Turing patterns. Nature 352, 610 (1991).
- Rietkerk, M. & van de Koppel, J. Regular pattern formation in real ecosystems. Trends in ecology & evolution 23, 169–75 (2008).
- Liu, Q.-X. et al. Phase separation explains a new class of self-organized spatial patterns in ecological systems. Proceedings of the National Academy of Sciences of the United States of America 110, 11905–10 (2013).
- Karlsson, M. et al. Formation of geometrically complex lipid nanotube-vesicle networks of higher-order topologies. Proceedings of the National Academy of Sciences of the United States of America 99, 11573–11578 (2002).
- Bignone, F. A. Structural Complexity of Early Embryos : A Study on the Nematode Caenorhabditis elegans. Journal of Biological Physics 27, 257–283 (2001).
- Barrat, A., Barthélemy, M. & Vespignani, A. Dynamical Processes on Complex Networks (Cambridge University Press, 2008).
- Colizza, V., Barrat, A., Barthélemy, M., Vespignani, A. & Barthe, M. The role of the airline transportation network in the prediction and predictability of global epidemics. Proceedings of the National Academy of Sciences of the United States of America 103, 2015–2020 (2006).
- Hanski, I. Metapopulation dynamics. Nature 396, 41–49 (1998).
- Urban, D. & Keitt, T. Landscape Connectivity : A Graph-Theoretic Perspective. Ecology 82, 1205–1218 (2001).
- Fortuna, M. A., Gómez-Rodríguez, C. & Bascompte, J. Spatial network structure and amphibian persistence in stochastic environments. Proceedings of The Royal Society B 273, 1429–34 (2006).
- Holland, M. D. & Hastings, A. Strong effect of dispersal network structure on ecological dynamics. Nature 456, 792–4 (2008).
- Boccaletti, S., Latora, V., Moreno, Y., Chavez, M. & Hwang, D. Complex networks: Structure and dynamics. Physics Reports 424, 175–308 (2006).
- Arenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y. & Zhou, C. Synchronization in complex networks. Physics Reports 469, 93–153 (2008).
- Kouvaris, N. E., Kori, H. & Mikhailov, A. S. Traveling and Pinned Fronts in Bistable Reaction-Diffusion Systems on Networks. PLoS ONE 7, e45029 (2012).
- Kouvaris, N. E., Isele, T., Mikhailov, A. S. & Schöll, E. Propagation failure of excitation waves on trees and random networks. EPL (Europhysics Letters) 106, 68001 (2014).
- Othmer, H. G. & Scriven, L. E. Instability and dynamic pattern in cellular networks. Journal of theoretical biology 32, 507–37 (1971).
- Othmer, H. G. & Scriven, L. E. Non-linear aspects of dynamic pattern in cellular networks. Journal of theoretical biology 43, 83–112 (1974).
- Horsthemke, W., Lam, K. & Moore, P. K. Network topology and Turing instabilities in small arrays of diffusively coupled reactors. Physics Letters A 328, 444–451 (2004).
- Moore, P. K. & Horsthemke, W. Localized patterns in homogeneous networks of diffusively coupled reactors. Physica D: Nonlinear Phenomena 206, 121–144 (2005).
- Nakao, H. & Mikhailov, A. S. Turing patterns in network-organized activatorâinhibitor systems. Nature Physics 6, 544–550 (2010).
- Wolfrum, M. The Turing bifurcation in network systems: Collective patterns and single differentiated nodes. Physica D: Nonlinear Phenomena 241, 1351–1357 (2012).
- Hata, S., Nakao, H. & Mikhailov, A. S. Dispersal-induced destabilization of metapopulations and oscillatory Turing patterns in ecological networks. Scientific reports 4, 3585 (2014).
- Boccaletti, S. et al. The structure and dynamics of multilayer networks. Physics Reports (2014).
- Gómez, S. et al. Diffusion Dynamics on Multiplex Networks. Physical Review Letters 110, 028701 (2013).
- Solé-Ribalta, A. et al. Spectral properties of the Laplacian of multiplex networks. Physical Review E 88, 032807 (2013).
- De Domenico, M. et al. Mathematical Formulation of Multilayer Networks. Physical Review X 3, 041022 (2013).
- Radicchi, F. & Arenas, A. Abrupt transition in the structural formation of interconnected networks. Nature Physics 9, 717–720 (2013).
- Louzada, V. H. P., Araújo, N. A. M., Andrade, J. S., Herrmann, H. J. & Araújo, N. A. M. Breathing synchronization in interconnected networks. Scientific reports 3, 3289 (2013).
- Granell, C., Gómez, S. & Arenas, A. Dynamical Interplay between Awareness and Epidemic Spreading in Multiplex Networks. Physical Review Letters 111, 128701 (2013).
- Asllani, M. et al. Turing patterns in multiplex networks. eprint arXiv:1406.6401v1 (2014).
- Stegeman, J. A., Elbers, A. R. W., Boum, A. & De Jong, M. C. M. Rate of inter-herd transmission of classical swine fever virus by different types of contact during the 1997-8 epidemic in The Netherlands. Epidemiology and infection 128, 285–91 (2002).
- Fang, L.-Q. et al. Geographical spread of SARS in mainland China. Tropical medicine & international health : TM & IH 14 (Suppl., 14–20 (2009).
- Xuan, Q., Du, F., Yu, L. & Chen, G. Reaction-diffusion processes and metapopulation models on duplex networks. Physical Review E 87, 032809 (2013).
- Cardillo, A. et al. Emergence of network features from multiplexity. Scientific reports 3, 1344 (2013).
- Buono, C., Alvarez-Zuzek, L. G., Macri, P. A. & Braunstein, L. A. Epidemics in partially overlapped multiplex networks. PloS one 9, e92200 (2014).
- Mimura, M. & Murray, J. D. On a diffusive prey-predator model which exhibits patchiness. Journal of theoretical biology 75, 249–62 (1978).
- Lee, K.-M., Kim, J. Y., Cho, W.-K., Goh, K.-I. & Kim, I.-M. Correlated multiplexity and connectivity of multiplex random networks. New Journal of Physics 14, 033027 (2012).
- Nicosia, V., Bianconi, G., Latora, V. & Barthélemy, M. Growing Multiplex Networks. Physical Review Letters 111, 058701 (2013).
- Kim, J. Y. & Goh, K.-I. Coevolution and Correlated Multiplexity in Multiplex Networks. Physical Review Letters 111, 058702 (2013).
- Kaluza, P., Kölzsch, A., Gastner, M. T. & Blasius, B. The complex network of global cargo ship movements. Journal of the Royal Society, Interface 7, 1093–103 (2010).
- Albert, R. & Barabási, A.-L. Statistical mechanics of complex networks. Reviews of modern physics 74, 47–97 (2002).