Turing patterns mediated by network topology in homogeneous active systems

TL;DR

This paper explores how Turing patterns can form in reaction-diffusion systems on various network topologies, revealing the influence of network structure and diffusion parameters on pattern stability and emergence.

Contribution

It provides a general mathematical framework for Turing pattern formation on networks, linking network topology, diffusion, and eigenvector localization to pattern stability.

Findings

Turing instability can occur in any network topology with appropriate diffusion.

Network degree fluctuations influence pattern stability and stochasticity.

Patterns are shaped by the interplay of diffusion coefficients and eigenvector properties.

Abstract

Mechanisms of pattern formation---of which the Turing instability is an archetype---constitute an important class of dynamical processes occurring in biological, ecological and chemical systems. Recently, it has been shown that the Turing instability can induce pattern formation in discrete media such as complex networks, opening up the intriguing possibility of exploring it as a generative mechanism in a plethora of socioeconomic contexts. Yet, much remains to be understood in terms of the precise connection between network topology and its role in inducing the patterns. Here, we present a general mathematical description of a two-species reaction-diffusion process occurring on different flavors of network topology. The dynamical equations are of the predator-prey class, that while traditionally used to model species population, has also been used to model competition between antagonistic ideas in social systems. We demonstrate that the Turing instability can be induced in any network topology, by tuning the diffusion of the competing species, or by altering network connectivity. The extent to which the emergent patterns reflect topological properties is determined by a complex interplay between the diffusion coefficients and the localization properties of the eigenvectors of the graph Laplacian. We find that networks with large degree fluctuations tend to have stable patterns over the space of initial perturbations, whereas patterns in more homogenous networks are purely stochastic.

Equations

$$\frac{du(\vec{x},t)}{dt}$$

$$\frac{dv(\vec{x},t)}{dt}$$

$$L_{ij}=A_{ij}-k_i{\delta }_{ij},$$

$$\frac{d{u}_i}{dt}$$

$$\frac{d{v}_i}{dt}$$

$${\lambda }_{\alpha }^2+b({\Lambda }_{\alpha }){\lambda }_{\alpha }+c({\Lambda }_{\alpha })=0,$$

$$\delta u_i(t_{\text{steady}})\approx \sum _{{\lambda }_{\alpha }>0}\exp \left({\lambda }_{\alpha }{t}_{\text{steady}}\right){\varphi }_i^{(\alpha )}.$$

$$P(\text{V})=\frac{{\sum }_{i=1}^N{V}_i^4}{({\sum }_{i=1}^N{V}_i^2{)}^2}.$$

$$⟨P⟩=\frac{1}{N}\sum _{\alpha =1}^{N}P({\text{V}}_{\alpha }),$$

$$⟨u_i⟩=\frac{1}{L}\sum _{l=1}^L{u}_i^l,$$

$$\frac{d{u}_i}{dt}=f(u_i,v_i)+D_{act}\sum _{j=1}^N{L}_{ij}{u}_j$$

$$\frac{d{v}_i}{dt}=g(u_i,v_i)+D_{inh}\sum _{j=1}^N{L}_{ij}{v}_j,$$

$$f(u,v)=\left(\frac{a+bu-u^2}{c}-v\right)u;g(u,v)=\left(u-(1+dv)\right)v.$$

$$\left(\begin{array}{c}\frac{d{u}_i}{dt}\\ \frac{d{v}_i}{dt}\end{array}\right)=J{|}_{u_0,v_0}\left(\begin{array}{c}{u}_i\\ v_i\end{array}\right)+\sum _{j=1}^N{L}_{ij}\left(\begin{array}{cc}{D}_{act}& 0\\ 0& D_{inh}\end{array}\right)\left(\begin{array}{c}{u}_i\\ v_i\end{array}\right)$$

$$J{|}_{u_0,v_0}=\left(\begin{array}{cc}{f}_u& f_v\\ g_u& g_v\end{array}\right),$$

$$\text{D}=D{|}_{u_0,v_0}=\left(\begin{array}{cc}{D}_{act}& 0\\ 0& D_{inh}\end{array}\right)=\left(\begin{array}{cc}ϵ& 0\\ 0& \gamma ,\end{array}\right)$$

$$\sum _{i=1}^N{L}_{ij}{\varphi }_j^{(\alpha )}={\Lambda }_{\alpha }{\varphi }_j^{(\alpha )}\text{where}\alpha =1\cdots N.$$

$$\frac{d\delta u_i}{dt}=f_u\delta u_i+f_v\delta v_i+\sum _{j=1}^N{L}_{ij}\text{D}\delta u_j$$

$$\frac{d{v}_i}{dt}=g_u\delta u_i+g_v\delta v_i+\sum _{j=1}^N{L}_{ij}\text{D}\delta v_j$$

$${\lambda }_{\alpha }\left(\begin{array}{c}1\\ B_{\alpha }\end{array}\right)=\left(\begin{array}{cc}{f}_u+ϵ{\Lambda }_{\alpha }& f_v{\Lambda }_{\alpha }\\ g_u{\Lambda }_{\alpha }& g_v+\gamma {\Lambda }_{\alpha }\end{array}\right)\left(\begin{array}{c}1\\ B_{\alpha }\end{array}\right)$$

$${\lambda }_{\alpha }^2-[Tr(J)+(ϵ+\gamma ){\Lambda }_{\alpha }]{\lambda }_{\alpha }+[ϵg_v+\gamma f_u]{\Lambda }_{\alpha }+Det(J)+ϵ\gamma {\Lambda }_{\alpha }^2=0,$$

$${\lambda }_{\alpha }^2+b({\Lambda }_{\alpha }){\lambda }_{\alpha }+c({\Lambda }_{\alpha })=0$$

$$b({\Lambda }_{\alpha })$$

$$c({\Lambda }_{\alpha })$$

$${\lambda }_{{\alpha }_1}$$

$${\lambda }_{{\alpha }_2}$$

$${\lambda }_{\alpha }({\Lambda }_{\alpha })=\frac{-b({\Lambda }_{\alpha })+\sqrt{b({\Lambda }_{\alpha }{)}^2-4c({\Lambda }_{\alpha })}}{2}$$

$$c({\Lambda }_{\alpha })=ϵ\gamma {\Lambda }_{\alpha }^2+[ϵg_v+\gamma f_u]{\Lambda }_{\alpha }+Det(J)<0.$$

$${\Lambda }_{{\alpha }_1}$$

$${\Lambda }_{{\alpha }_2}$$

Full text

I Introduction

Pattern formation is a fundamental natural phenomenon that has abundant examples in biological, ecological and chemical processes Castets et al. (1990); Ouyang and Swinney (1991); Kondo and Asai (1995); Hans Meinhardt (2000); Sayama et al. (2003); Maini et al. (2006). The first mathematical description was proposed by Alan Turing Turing (1952) where he demonstrated the spontaneous emergence of periodic spatial patterns from a homogeneous equilibrium, driven by the reaction and diffusion of two chemical species corresponding to activators and inhibitors Baurmann et al. (2007); Rietkerk and van de Koppel (2008); Harris et al. (2005).

The existence of pattern formation on discrete media was first introduced by Othmer and Scriven in the context of morphogens diffusing over a network of intercellular connections Othmer and Scriven (1971). Plane-wave driven Turing instability was studied on one and two dimensional lattices, where the wave-functions and wave-numbers in continuous media are replaced by their discrete analogs corresponding to the eigenvalue and eigenvectors of the Laplacian matrix. In recent years, Turing patterns have been shown to exist in reaction-diffusion processes occurring on complex networks, where nodes in the graph are assigned an initial concentration of chemical species, and diffusion occurs along the edges connecting the nodes Ide et al. (2016); Wolfrum (2012). A small perturbation to the uniform state triggers the growth of Turing patterns above a critical threshold, corresponding to the ratio of the diffusion constants of the respective species. The patterns in this context, correspond to distinct populations of nodes differentiated by their levels of chemical concentrations, and exhibit properties quite different from the classical case. For instance, multiple coexisting stationary states can occur, and hysteresis effects are present, indicating that the patterns are not particularly robust over the space of initial perturbations Nakao and Mikhailov (2010).

Turing patterns have also been observed in more exotic variants of network topologies, such as directed, multiplex and time-evolving networks Asllani et al. (2014a, b); Petit et al. (2017). Indeed, given the widespread prevalence of networks across a plethora of socioeconomic, biological and technological systems Newman (2010); Albert and Barabási (2002); Pastor-satorras and Vespignani (2004); Cohen and Havlin (2010), the developed framework is being used to model interesting dynamical phenomena; in Vidal-Franco et al. (2017), predator-prey Vanag and Epstein (2009) dynamics is used along with cross-diffusion (in addition to ordinary diffusion), to model the evolution of two competing languages on a scale free network Caldarelli (2010). In Hata et al. (2014) three species interactions are used to model ecological meta-populations where concentrations correspond to the population densities of species in a food web.

While a lot has been studied regarding the application and existence of Turing patterns in different flavors of large networked systems, there remains much to be understood regarding the precise role of network topology. For instance, it is known that heterogeneities in the connectivity structures (degree-distributions) of networks play a major role in dynamical processes such as percolation, diffusion and epidemic spreading among others Callaway et al. (2000); Pastor-Satorras and Vespignani (2001); Colizza and Vespignani (2007); Dorogovtsev et al. (2008); Ghoshal and Barabási (2011). Although a connection has been established between degree-fluctuations and the localization of the Laplacian eigenvectors that play a role in pattern-formation Nakao and Mikhailov (2010); Hata and Nakao (2017), the space of topologies in which this has been studied is rather limited. Previous investigation on large networks have typically interpolated between the Erdős-Rényi (ER) random graph Erdős and Rényi (1959) (where fluctuations vanish in the thermodynamic limit) and Barabási-Albert networks that provide a very limited range for probing heterogenous structure Barabasi and Albert (1999). Furthermore, the properties of the patterns have primarily been studied at the critical threshold (the point at which the Turing pattern emerges) and a rich set of potential dynamical behavior above threshold has not been studied in detail.

Here, we present an exhaustive analysis of the connection between network topology and the dynamical parameters, uncovering new insights on the nature of Turing patterns in networks. We find that in certain instances, network topology plays a role not only in instigating the spontaneous differentiation of nodes, but also provides clues to the steady state concentrations. Specifically we show, that while the onset of the pattern can be controlled by tuning the average connectivity of the network, the nature of the final patterns is determined by a complex interplay between degree fluctuations and the diffusion coefficients. In those networks with peaked degree distributions, the patterns are purely stochastic, being determined primarily by perturbations to the initial uniform state. On the other hand, networks with large degree fluctuations induce patterns that are quite robust to initial perturbations—the structure of the network is correlated strongly with the final concentrations of the nodes. This correlation gradually vanishes as a function of increasing the diffusion coefficients such that the patterns are indistinguishable from those formed by homogenous network topologies. We end by discussing the implications of our results in applying reaction-diffusion networks to model real world systems.

II Reaction-Diffusion Model and Turing Instabilities on Networks

Pattern-forming, reaction-diffusion systems, in continuous media, are typically described by a set of two partial differential equations of the form

[TABLE]

Here $u(\vec{x},t),v(\vec{x},t)$ correspond to the local concentrations of two chemical species, $f(u,v),g(u,v)$ specify their local dynamics and $\epsilon,\gamma$ are the corresponding diffusion coefficients Prigogine and Lefever (1968); Yin et al. (2013); Lou et al. (1998). Typically $u$ corresponds to an activator, that grows through autocatalytic growth, and $v$ an inhibitor that suppresses $u$. The system is initially considered to be at an uniform steady state ($u_{0}$, $v_{0}$) where $f(u_{0},v_{0})=g(u_{0},v_{0})=0$. The Turing instability spontaneously emerges above a critical threshold of the ratio of diffusion constants, $\sigma=\gamma/\epsilon$ and corresponds to alternating spatial regions of high- and low-concentrations of $u$.

A similar set of equations can be used if the system is not continuous but instead is composed of $N$ independent nodes that interact via diffusive transport over $m$ edges Ide et al. (2016); Mccullen and Wagenknecht (2016); Nakao and Mikhailov (2010); Vidal-Franco et al. (2017). The analog of the operator $\nabla^{2}$ is now the Laplacian matrix,

[TABLE]

where $A_{ij}$ is the symmetric adjacency matrix, whose elements are 1 if there is an edge between nodes $i$ and $j$ and 0 otherwise; $k_{i}=\sum_{j}A_{ij}$ is the degree of node $i$ (number of connections) and $\langle k\rangle=2m/N$ is the average degree of the network. The diffusive transport of chemical species to node $i$ is the sum of all the incoming fluxes from its neighbors $j$ and proportional to their concentration differences. With this modification, Eq. (1) is transformed into a set of $N$ differential equations ($i=1\dots N$) thus,

[TABLE]

For the purposes of our analysis, we use the Mimura-Murray model which has been used to model predator-prey populations Mimura and Murray (1978). Our choice is motivated by its ubiquity in studying reaction-diffusion dynamics in networked systems, however our results are broadly applicable to a wider range of dynamics.

The Turing instability is examined through a linear stability analysis of the uniform stationary state with respect to non-uniform perturbations. While in continuous media, uniform perturbations are decomposed into a set of spatial Fourier modes representing plane waves with different wave-numbers, in networks their analogs are the eigenvectors $\phi_{i}^{(\alpha)}$ and eigenvalues $\Lambda_{\alpha}$ of the Laplacian matrix, where $\alpha=1,\ldots N$ corresponds to the eigenmode. The eigenvalues $\Lambda_{\alpha}$ are sorted in decreasing order $\Lambda_{1}>\Lambda_{2}\dots>\Lambda_{N}$ and the first eigenvalue is always zero ($\Lambda_{1}=0$). Introducing small perturbations $(\delta u_{i},\delta v_{i})$, substituting into Eq. (3), and expanding over the set of the Laplacian eigenvectors, the linear growth rate $\lambda_{\alpha}$ for each node is calculated from a polynomial equation of the form

[TABLE]

where $b(\Lambda_{\alpha}),c(\Lambda_{\alpha})$ are functions of the diffusion coefficients as well as $f,g$. The Turing instability occurs when at least one of the modes becomes unstable, indicated by Re $\lambda_{\alpha}>0$ which happens when $c(\Lambda_{\alpha})<0$ (See Appendix A for details of the calculation). In Fig. 1 we show an example of Turing patterns in an ER network of 100 nodes. Panel a shows the initial concentration $u_{i}$ for each node $i$ (sorted in decreasing order of degree) and in b we show the network with the nodes colored according to the concentrations. After perturbations, spontaneous differentiation occurs and the system evolves to its final state as shown in panels c and d.

Above threshold, the instability range is parameterized by the roots $\Lambda_{\alpha_{1}},\Lambda_{\alpha_{2}}$ of the parabola $c(\Lambda_{\alpha})$ which takes on negative values in the range $\Lambda_{\alpha_{1}}\leq c(\Lambda_{\alpha})\leq\Lambda_{\alpha_{2}}$ (See Appendix A). The emergence of positive growth rates $\lambda_{\alpha}$ corresponds to the overlap of the Laplacian eigenvalues with the instability regime. Given that $Tr(L)=\sum_{i}k_{i}=\sum_{\alpha}\Lambda_{\alpha}$, we have that $\langle k\rangle=\langle\Lambda\rangle$. Thus for a fixed set of diffusion coefficients and dynamical parameters, Turing instabilities can be triggered in any flavor of network topology by tuning the average connectivity of the system. Conversely, for fixed network topology, the diffusion coefficients can be tuned (within physically meaningful bounds) to precipitate the instabilities.

Indeed, it has been pointed out that Turing patterns can be enhanced, as measured by the amplitude of separation $A=[\sum_{i=1}^{N}{(u_{i}-u_{0})^{2}+(v_{i}-v_{0})}]^{1/2}$, by increasing the connectivity of the network Vidal-Franco et al. (2017). However the effect is far more nuanced as shown in Fig. 2a where we plot the growth factors as a function of the Laplacian eigenvalues for four different ER networks of varying connectivity. The instability range (delineated by $\Lambda_{\alpha_{1}},\Lambda_{\alpha_{2}}$) is fixed by the diffusion coefficients. For a network of $\langle k\rangle=50$, we see that none of the growth factors are positive and consequently there is no change in the initial concentrations (Fig. 2b). As the average connectivity increases to 200, we see that all growth factors are positive and spontaneous differentiation occurs, with a large separation between the concentrations (Fig. 2c). As the connectivity is increased even further, there is a mixture of positive and negative growth factors, leading to less pronounced differentiation (Fig. 2d). Finally, as connectivity increases further, all growth factors are negative, and no differentiation is apparent (Fig. 2e). Note that the maximum of the amplitude $A$ occurs near the peak of the instability curve. Thus, while $\langle k\rangle$ can be tuned to trigger the Turing instability, it is the location of the eigenmodes in the instability regime that determines the trend and amplitude of the differentiation.

III Eigenvector Localization

Having established the role of the average connectivity in triggering the Turing instability, we next investigate the connection with degree heterogeneities. To do so, it is instructive to outline the role that the eigenmodes and eigenvectors play in the evolution of the system. After perturbing the system, temporal evolution progresses via the change in chemical concentrations until they reach their steady state. The instantaneous evolution is characterized by $u_{i}(t)=u_{0}+\delta u_{i}(t)$, where $\delta u_{i}(t)=\sum_{\alpha=1}^{N}a_{\alpha}\exp\left(\lambda_{\alpha}t\right)\phi_{i}^{(\alpha)}$ for some constants $a_{\alpha}$. For negative growth factors, $\lambda_{\alpha}$, the perturbations to the initial state vanish, implying that the evolution of the system can be characterized by considering contributions to the sum for only those modes where $\lambda_{\alpha}>0$.

Without loss of generality, we set $a_{\alpha}=1$ (given the negligible contribution to the sum), such that in the steady state we have

[TABLE]

Thus, differentiation of nodes is associated with the properties of the eigenvectors corresponding to positive growth factors. The extent to which the eigenvectors are connected to network topology depends on its localization properties McGraw and Menzinger (2008); Pastor-Satorras and Castellano (2016); Jalan et al. (2011). A localized vector has the majority of its normalization weight concentrated on a small subset of its components, whereas a delocalized vector has its weights distributed relatively uniformly across all of its components. An ideal way to quantify the localization of a vector V of size $N$, is to measure its Inverse Participation ratio (IPR) defined as

[TABLE]

The IPR lies in the range $1/N\leq P\left(\textbf{V}\right)\leq 1$ with the upper limit corresponding to all weights localized on a single component, whereas the lower limit represents the situation where the weights are uniformly distributed over all components.

There is a clear correlation between node degree and the Laplacian eigenvalues in networks Hata and Nakao (2017). In addition, eigenvectors associated with higher degree nodes are strongly localized, whereas the localization is weaker for nodes with fewer connections. To see the effects clearly in terms of network topology, we generated a series of networks using the BA model, the ER random graph and power law (PL) networks ($p_{k}\sim k^{-\beta}$) using the configuration model Molloy and Reed (1995, 1998). The networks are created with different random seeds and the relevant parameters are tuned such that they all have the same average degree $\langle k\rangle$ but differ significantly in their degree fluctuations $\sigma_{k}=\sqrt{\langle k^{2}\rangle-\langle k\rangle^{2}}$. For each generated network, we compute the average IPR

[TABLE]

and plot it as a function of $\sigma_{k}$ as shown in Fig. 3. A clear monotonic trend is seen, where the PL networks with large degree fluctuations have significantly higher $\langle P\rangle$ as compared to the ER networks. Taking the average of all points (shown as inset) shows a marked difference between the BA networks ($p_{k}\sim k^{-3}$) and the PL graphs ($p_{k}\sim k^{-2.1}$), given that $\sigma_{k}\rightarrow\infty$ for $\beta<3$ in the thermodynamic limit. Thus increasing degree heterogeneities leads to increasing localization of the Laplacian eigenvectors of the network.

Earlier work Nakao and Mikhailov (2010); Hata et al. (2014) connect the eigenvector $\phi_{i}^{(\alpha_{c})}$ associated with the critical eigenmode $\alpha_{c}$(corresponding to the onset of the Turing instability) with the differentiation patterns in early phases of the evolution, with the steady-state concentrations $\delta u_{i}(t_{\textrm{steady}})$ and $\delta v_{i}(t_{\textrm{steady}})$, being determined by non-linear effects. However, as we show next, if contributions to Eq. (5) are composed of localized eigenvectors, then a connection can be made to the resultant Turing Patterns and topological properties of the network (quantified by $\sigma_{k}$). If on the other hand, most of the sum is determined by non-localized vectors then the effect of the network topology is washed out and patterns emerge randomly, primarily as function of initial perturbations to the uniform background.

IV Numerical Results

Next, we investigate the typology of the emergent Turing patterns as a function of degree heterogeneity. We start with the limiting cases of the star- and complete-graphs and fill the spectrum, interpolating between the ER and PL networks.

IV.1 Limiting cases

A star graph consists of $N$ nodes with a central node connected to $N-1$ peripheral nodes, with no additional connection between them. The average degree $\langle k\rangle=2\left(1-1/N\right)$ and the second moment is $\langle k^{2}\rangle=N-1$. Consequently the fluctuations are extensive and of the form $\sigma_{k}^{\textrm{star}}\approx\sqrt{N}$. On the other end of the spectrum is the complete graph where all pairs of nodes are connected to each other. In this case $\langle k\rangle=N-1$ and $\langle k^{2}\rangle=(N-1)^{2}$ and therefore $\sigma_{k}^{\textrm{complete}}=0$.

We simulate Mimura-Murray dynamics on each of the networks with diffusion constants $\epsilon=2\times 10^{-3}$ and $\sigma=20$ to induce differentiation via the Turing instability. Plotting the growth factors as a function of the Laplacian eigenvalues (Fig. 4) reveals that the star graph has one positive growth factor, whereas the complete graph has $N-1$ positive growth factors. This is reflected in the final concentrations for each of the networks as shown in Fig.  5a,c. In the star graph the only node that differentiates is the central node, whereas in the complete graph two separate clusters of roughly equal size emerge. Approximating the final concentrations using Eq. (5) we see that in Fig. 5b, the Turing pattern formed in the star graph is described entirely by the eigenvector corresponding to the central node with the maximum number of connections. Conversely, as seen in Fig. 5d, there appears to be no connection with the final Turing pattern and the structure of the eigenvectors associated with differentiation. Nodes differentiate randomly as a function of the initial perturbations. This sharply contrasting behavior is strongly connected to the varying localization properties in each flavor of network. Calculating the average IPR of the eigenvectors associated with positive growth factors $\langle P_{\lambda_{\alpha^{+}}}\rangle$, we see that the star graph has a perfectly localized eigenvector $\langle P_{\lambda_{\alpha^{+}}}\rangle=1$, whereas the complete graph has $\langle P_{\lambda_{\alpha^{+}}}\rangle=4.7\times 10^{-2}$.

IV.2 Dynamics on ER and PL Networks

The cases considered thus far constitute the two end points of the spectrum of network structure in terms of their degree heterogeneities. For maximal degree fluctuations, the Turing pattern is determined entirely by topological effects, whereas for no degree fluctuations, the patterns are random and have little-to-no connection with network topology. To determine the extent to which topological effects manifest themselves, we next interpolate between these two limits, starting with the ER network and then two power law networks generated via the configuration model; one with non-extensive fluctuations $\left(p_{k}\sim k^{-3.1}\right)$, and the other with extensive fluctuations $(p_{k}\sim k^{-2.1})$. The relevant parameters are tuned such that all networks have the same average degree $\langle k\rangle\approx 7$.

In Fig. 6 we show the results of our simulations for $\epsilon=1.1\times 10^{-1}$, $\sigma=17$ for all three network topologies. Given that both the instability range and the average degree is fixed, the observed variations are purely a function of $\langle k^{2}\rangle$ and the attendant localization properties. For the PL network with extensive fluctuations, the pattern is akin to that seen for the star graph in the sense that there is a close connection between the final concentrations (Fig. 6a) and the contributions from the eigenvectors associated with node differentiation (Fig. 6b), although it is not as pronounced as in the star graph, given the comparatively lower localization $\left(\langle P_{\lambda_{\alpha^{+}}}\rangle=2\times 10^{-1}\right)$. This correspondence gradually vanishes as fluctuations dissipate and the localization of the eigenvectors decrease $\left(\langle P_{\lambda_{\alpha^{+}}}\rangle=6.1\times 10^{-2}\right)$, as seen for the PL network with finite fluctuations (Figs. 6c,d). Eventually, as we get to the ER network, the behavior is much like that for the complete graph (Figs. 6c,d), whereby the final concentrations appear uncorrelated with the properties of the relevant eigenvectors $\left(\langle P_{\lambda_{\alpha^{+}}}\rangle=6.1\times 10^{-2}\right)$, and are more a function of the random perturbations to the initial state.

IV.3 Interplay between diffusion and topology

Next, we investigate the role played by the diffusion constants in the generated patterns. Our strategy thus far has been to keep the dynamical parameters fixed —therefore fixing the instability range $\left(\Lambda_{\alpha_{1}},\Lambda_{\alpha_{2}}\right)$—while probing for the effects of network topology. Now we fix the network topology (and the corresponding eigenspectrum) and tune the diffusion parameters such that we can change the number of positive growth factors that lead to differentiation. In Fig. 7a, we plot the IPR of all the eigenmodes (Eq. (6)) as a function of the eigenvalues of the power law network with degree distribution $p_{k}\sim k^{-2.1}$. The eigenvalues associated with positive growth factors corresponding to the parameters used in Fig. 6 are colored blue. Changing the values of the diffusion constants to $\epsilon=7\times 10^{-2}$ and $\sigma=55$ leads to a wider instability regime, and the corresponding unstable eigenvalues are colored green. The distribution of positive growth factors for both instances is shown as inset.

Even though the networks are identical, one sees a larger set of unstable eigenvalues associated with delocalized eigenvectors in the case of the wider instability range. The IPR decreases from $2.2\times 10^{-1}$ to $1.6\times 10^{-2}$ and its effect is reflected in the final concentrations shown in Fig. 7c which are markedly different than seen in Fig. 6a. Furthermore, the differentiation patterns are more irregular and are comparatively uncorrelated with the structure of the associated eigenvectors as seen in Fig. 7d. Indeed, the trend is more akin to seen that for the ER network, where the resultant patterns are purely stochastic and dependent on the initial perturbations to the system. Plotting the localization as a function of increasing $\sigma=\gamma/\epsilon$ (which leads to an increase in the instability range), shows a clear montonically decreasing trend of the eigenvector localization. Thus even in networks where there is a connection between the topological structure and the resultant differentiation patterns, this effect gets washed out depending on the diffusion coefficients. In the specific case of the Mimura-Murray model, this happens when the inhibitor diffuses much faster than the activator.

IV.4 Averaged patterns

A clear way to see the effect of topology on the final Turing pattern is to generate multiple realizations of the dynamical process as a function of initial perturbations to the uniform background. To probe for this, we generated different realizations of the dynamics for the network shown in Fig. 7, for the two different instability ranges, as well as the ER network shown in Fig. 6 and took the average of the resultant concentrations for each node. That is for each node $i$ in the relevant network, we computed

[TABLE]

where $l$ corresponds to a single instance of random perturbations, and $L=100$ is the total number of realizations of this process. In Fig. 8 we plot the results of our analysis. Panel a shows $\langle u_{i}\rangle$ for the PL network in the lower instability range. While some fluctuations are apparent, these are associated with the relatively narrow set of differentiating nodes whose concentrations are well approximated with their corresponding eigenvectors. The Turing pattern for the network in a single realization (Fig. 8d) is almost identical to the averaged concentration over a 100 different realizations (Fig. 8g). The degree heterogeneity in this case leads to an element of determinism in the final Turing pattern, that is robust to variations in the perturbations introduced to the initial state. Increasing the diffusion of the inhibitor relative to the activator for the same network, results in a very different situation as seen in Fig. 8e. Fluctuations in the average concentration dramatically increase for practically all nodes, indicating little correlation between network topology and the resultant pattern. This is reflected when comparing the Turing pattern for one realization (Fig. 8e) to the average over multiple realizations (Fig. 8h).The average of the concentrations bear no resemblance to any one instance, indicating that the patterns are purely stochastic. A similar trend is seen for the ER network with the same diffusion coefficients as the first PL network. Due to negligible degree heterogeneities, the topology has little-to-no effect on the Turing patterns (Fig. 8c,f,i).

V Discussion

Taken together, our results shed new insight on the influence of network topology in the instigation of Turing patterns and the eventual differentiation of nodes in the steady state. The conditions for instability are set by the diffusion coefficients, and the instability emerges when the eigen-spectrum of the network overlaps with the instability regime. In principle node differentiation can be induced in any flavor of network topology, by tuning the average connectivity of the system, which controls the extent of overlap between the eigen-spectrum and the instability regime; conversely for a fixed network topology, one can tune the dynamical parameters—within meaningful bounds—to once again generate Turing patterns. The nature of the eventual pattern is connected deeply with the localization properties of the eigenvectors of the associated network. Networks with large degree fluctuations have highly localized eigenvectors, which for certain regimes of the diffusion parameters correlate strongly with the steady state Turing pattern. Indeed, in the case of the star graph which has maximal degree fluctuations, the differentiation is determined entirely by the maximally localized eigenvector, associated with the central node. On the other end of the spectrum, pattern formation in networks with little-to-no degree fluctuations is entirely stochastic and essentially independent of the specific properties of the eigenvectors of the network.

Interpolating between heavy-tailed distributions with extensive fluctuations and peaked distributions with non-extensive fluctuations, we see a gradual decoupling between the emergent Turing pattern and the topological properties of the network, with the former being robust to the space of initial perturbations to the uniform background, whereas the latter being entirely influenced by it. Even for those networks with strongly localized eigenvectors, we find that its influence in pattern formation is diminished by tuning the diffusion coefficients. For the specific case of the Mimura-Murray model, the link between topology and pattern formation vanishes either when degree fluctuations are suppressed, or when the inhibitor diffuses at a much faster rate than the activator. In both cases, the extent to which topology matters, can be quantified by calculating the inverse proportionality ratio (IPR) of eigenvectors associated with unstable eigenvalues.

Our results have interesting implications for the application of Turing patterns in networks, in service of modeling real world dynamical processes. One of the key prevalent phenomenon in recent times is the polarization of opinions, and socioeconomic segregation across many social systems. In the case of opinion formation, one possibility, would be to model the evolution of mutualistic-antagonistic dynamics on social networks, which is similar in flavor to the predator-prey model. It has been suggested that the formation of echo chambers and polarized opinions, is a function of both the speed of information spread (diffusion) and the increased connectivity afforded by the prevalence of social media Barberá et al. (2015); Vicario et al. (2016). This is indeed borne out by Fig. 2, where we see an increase in the amplitude of the differentiation as connectivity in the network is increased (for that particular combination of diffusion parameters). However, perhaps encouragingly, the differentiation vanishes as the network is even more connected. To the extent that opinion formation can be modeled by the framework presented here, tuning network connectivity appears to be a proscriptive tool, in terms of preventing segregation.

Furthermore, beyond the question of whether Turing patterns exist on networks, it is worth investigating the extent to which the differentiation properties are associated with specific nodes. In other words, given a set of dynamical parameters, the extent to which a given node takes on a characteristic concentration (opinion, socioeconomic state, language adoption etc.) in the final state. As Fig. 8 indicates, for highly heterogenous network topologies and small diffusion rates, the network has a stable hierarchical structure in terms of the final concentrations. The situation persists for the case when the mutualistic opinion diffuses much faster than the antagonistic one. As the rate of antagonistic opinions spread faster relative to mutualistic ones, the hierarchical structure is lost and while multiple opinions form, they are stochastic with respect to the population on which the information is spreading. Of course, similar considerations apply when modeling ecological species Hata et al. (2014) or linguistic adoption Vidal-Franco et al. (2017).

Acknowledgements.

SM, GG and JP acknowledge support from the Global Research Network program through the Ministry of Education of the Republic of Korea and the National Research Foundation of Korea (NRF-2016S1A2A2911945). APM and MMJ are supported by the Spanish Ministerio de Economía y Competitividad and European Regional Development Fund under contract MAT2015-71119-R AEI/FEDER, UE, and by Xunta de Galicia under Research Grant No. GPC2015/014. GG and SM also thank the University of Rochester for financial support.

Appendix A Dynamical Model and Linear Stability Analysis

The general form of the differential equations used here are,

[TABLE]

[TABLE]

where $D_{act}=\epsilon$ and $D_{inh}=\gamma$, with $\gamma=\sigma\epsilon$. For the Mimura-Murray model

[TABLE]

Setting the constants $a=35,b=16,c=9,d=2/5$ yields a fixed point concentration $(u_{0},v_{0})=(5,10)$. The equations are evaluated using the $4^{th}$ order Runge-Kutta Method. Initial non-uniform perturbations to node concentrations are $\delta u(0)=u_{0}+0.1u_{0}r$ and $\delta v(0)=v_{0}+0.1v_{0}r$ where $r$ is a uniform random number in the range $-1<r<1$.

The equations can be cast in matrix form thus,

[TABLE]

where

[TABLE]

is the Jacobian matrix and

[TABLE]

is the diffusion matrix both evaluated at the fixed point $(u_{0},v_{0})$.

The corresponding eigenvalue equation for the Laplacian matrix is:

[TABLE]

Introducing perturbations around the fixed points: $(u_{i},v_{i})=(u_{0},v_{0})+(\delta u_{i},\delta v_{i})$ and substituting into Eq. (A2) yields linearized equations for $\delta u_{i}$ and $\delta v_{i}$,

[TABLE]

[TABLE]

The perturbations can be expanded over the set of Laplacian eigenvectors as $\delta u_{i}(t)=\sum_{\alpha=1}^{N}c_{\alpha}exp[\lambda_{\alpha}t]\phi_{i}^{(\alpha)}$ and $\delta v_{i}(t)=\sum_{\alpha=1}^{N}c_{\alpha}B_{\alpha}exp[\lambda_{\alpha}t]\phi_{i}^{(\alpha)}$. Substituting these into equations (A3) and (A4) we obtain the following eigenvalue equation:

[TABLE]

The characteristic equation of this system is given by:

[TABLE]

where $Tr(J)=f_{u}+g_{v}$ and $Det(J)=f_{u}g_{v}-f_{v}g_{u}$.

More compactly this can be written as

[TABLE]

where

[TABLE]

The roots of Eq. (A7) are

[TABLE]

In the absence of diffusion, the system is stable and all eigenvalues are negative. Setting all elements of the diffusion matrix to zero yields $b(\Lambda_{\alpha})=-Tr(J)$ and $c(\Lambda_{\alpha})=Det(J)$. If both solutions are negative then, $c(\Lambda_{\alpha})>0$ and $-b(\Lambda_{\alpha})<0$, thus instability require the following conditions to be satisfied,

• •

$Det(J)=f_{u}g_{v}-f_{v}g_{u}>0$

• •

$Tr(J)=f_{u}+g_{v}<0$

From the second inequality, we can deduce that $Tr(J)+(\epsilon+\gamma)\Lambda_{\alpha}<f_{u}+g_{v}<0$ and therefore $\epsilon+\gamma>0$ (Laplacian eigenvalues are always non-positive). Then, the only positive solution is

[TABLE]

In order to induce differentiation, we need $\lambda_{\alpha}(\Lambda_{\alpha})>0$ for some $\alpha$. So, we write: $\frac{-b(\Lambda_{\alpha})+\sqrt{b(\Lambda_{\alpha})^{2}-4c(\Lambda_{\alpha})}}{2}>0$ which implies $c(\Lambda_{\alpha})<0$. This means that

[TABLE]

Since all Laplacian eigenvalues are non-positive, we have: $\epsilon g_{v}+f_{u}\gamma>0$, and the condition $\epsilon g_{v}+f_{u}\gamma-4\epsilon\gamma Det(J)>0$ guarantees negative $c(\Lambda_{\alpha})$ values (we take $\epsilon\gamma>0$ for an upward opening parabola). The roots of $c(\Lambda_{\alpha})$ are then,

[TABLE]

Thus, $c(\Lambda_{\alpha})$ is negative if some of the Laplacian eigenvalues are in the range [$\Lambda_{\alpha_{1}}$,$\Lambda_{\alpha_{2}}$] and that guarantees the existence of corresponding growth factors $\lambda_{\alpha}$ to be positive. This is shown in Fig. 9b where the eigenvalues of a PL and ER network are plotted on the curve $c(\Lambda_{\alpha})$. The ER network has all of its eigenvalues (except for $\Lambda_{1}=0$) in the instability range [$\Lambda_{\alpha_{1}}$,$\Lambda_{\alpha_{2}}$] whereas the PL network has some of its eigenvalues inside the range and some outside due to its larger degree fluctuation. In Fig. 9a,c we plot the eigenspectrum of both networks.