Turing patterns mediated by network topology in homogeneous active systems

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 features in social contexts. 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.

Figure 1

Evolution of Turing patterns on a ER network with $N=100$. (a) The concentration $u_i$ of node $i$ (sorted in decreasing order of degree $k_i$) immediately after perturbing the uniform background. (b) The network with nodes colored according to concentration showing a single population. (c) The concentrations after the system has evolved to its steady state, indicating two distinct populations of nodes corresponding to concentrations lying above and below the initial uniform concentration (dashed black line). (d) The corresponding network representation.

Figure 2

(a) Eigenvalue distributions of four ER networks with $N={10}^3$ and varying average degrees (indicated by color) plotted on the dispersion relation curve $\lambda \left({\mathrm{\Lambda }}_{\alpha }\right)$ [Eq. (4)], shown as dashed line, corresponding to the Mimura-Murray model with diffusion coefficients $\varepsilon =5\times {10}^{-3}$ and $\sigma =30$. The instability range is marked by ${\mathrm{\Lambda }}_{{\alpha }_1},{\mathrm{\Lambda }}_{{\alpha }_2}$ and is fixed by the diffusion coefficients. The extent of overlap of the growth factors with the instability regime is controlled by $\langle k\rangle$ which shows a nonmonotonic trend. [(b)–(e)] The corresponding steady-state concentrations for each of the considered networks: $\langle k\rangle =50,200,520,720$; the black dashed line represents the initial uniform concentrations. Maximum differentiation occurs at $\langle k\rangle =200$ with an amplitude $A=97.52$, which corresponds to the growth factors being distributed around the peak of the instability curve.

Figure 3

Network topology and eigenvector localization. The average IPR $\langle P\rangle$ [Eq. (7)] as a function of degree fluctuations ${\sigma }_k=\sqrt{\langle k^2\rangle -{\langle k\rangle }^2}$ for BA, ER, and PL networks with $N={10}^3$ and $\langle k\rangle \approx 6.0$. While the ER network has a binomial degree distribution, the BA network goes as $p_k\sim k^{-3}$ in its tail and the PL network has a distribution of the form $p_k=k^{-\beta }/\zeta (\beta ,k_{\min })$, where $\zeta$ is the Hurwitz zeta function and $\beta =2.1,k_{\min }=2$. The inset shows the average of data points indicating a clear monotonic trend of the eigenvector localization with increasing degree heterogeneity.

Figure 4

Eigenvalue distributions of the star and dense graphs with $N={10}^3$, plotted on the dispersion relation curve $\lambda \left({\mathrm{\Lambda }}_{\alpha }\right)$ (shown as dashed line). The diffusion coefficients are $\varepsilon =2\times {10}^{-3}$ and $\sigma =20$. The instability range is marked by ${\mathrm{\Lambda }}_{{\alpha }_1},{\mathrm{\Lambda }}_{{\alpha }_2}$. The star graph has only one positive growth factor, whereas the dense graph has 999 growth factors.

Figure 5

Steady-state chemical concentrations $u_i$ for each node $i$ in the star graph (a) and dense graph (c). Nodes are arranged in decreasing order of degree and the dynamical parameters are the same as Fig. 4. The dashed black line represents the initial concentrations. [(b) and (d)] The contributions of the eigenvectors associated with positive growth factors [Eq. (8)] for each of the networks. Nodes are colored according to their final concentration values. The final Turing pattern for the star graph is determined exclusively by the highly localized leading eigenvector corresponding to the highest degree node $\left(\langle P_{{\alpha }_{+}}\rangle =1\right)$. There appears little correlation between the graph topology and the final Turing pattern in the dense graph, with nodes differentiating randomly $\left(\langle P_{{\alpha }_{+}}\rangle =1.3\times {10}^{-2}\right)$.

Figure 6

Steady-state chemical concentrations $u_i$ for each node $i$ in the power-law graph with parameters $\beta =2.1,k_{\min }=2$ (a), power-law graph with parameters $\beta =3.1,k_{\min }=4$ (c) and ER network (e). The diffusion coefficients are $\varepsilon =1.1\times {10}^{-1},\sigma =17$. All networks have $\langle k\rangle \approx 7$ and nodes are arranged in decreasing order of degree. The dashed black line represents the initial concentrations. [(b), (d), and (e)] The contributions of the eigenvectors associated with positive growth factors [Eq. (8)] for each of the networks. Nodes are colored according to their final concentration values. The localization of the eigenvectors associated with positive growth factors $\langle P_{{\alpha }_{+}}\rangle$ are from (top to bottom) $2.2\times {10}^{-1},6.1\times {10}^{-2},1.7\times {10}^{-2}$.

Figure 7

Effect of diffusion on Turing patterns. (a) The IPR [Eq. (6)] as a function of the Laplacian eigenvalues of a PL network with $\beta =2.1,k_{\min }=2$, for different values of the diffusion coefficient: Range 1 corresponds to $\varepsilon =1.1\times {10}^{-1},\sigma =17$ and Range 2 corresponds to the same value for $\varepsilon$, but $\sigma =50$. Positive growth factors on the dispersion curve (shown as inset) are colored purple and orange, and the localization of the associated eigenvectors decreases from Range 1 $\left(\langle P_{{\alpha }_{+}}\rangle =2.2\times {10}^{-1}\right)$ to Range 2 $\left(\langle P_{{\alpha }_{+}}\rangle =1.6\times {10}^{-1}\right)$. (b) The localization of the unstable modes plotted as a function of $\sigma$ showing a clear delocalization trend (and therefore more stochasticity in pattern formation). (c) The steady-state concentrations of the nodes corresponding to the wider instability range, showing a markedly different pattern than seen in Fig. 6. (d) The contributions from the eigenvectors associated with the positive growth factor. Nodes colored according to concentration.

Figure 8

Averaged steady-state concentrations $\langle u_i\rangle$ Eq. (10) over multiple realizations of the dynamics: (a) The PL network from Fig. 7 for Range 1. (b) The same network over Range 2. (c) The ER network with the same parameters corresponding to Range 1. Node indices are sorted in decreasing order of degree, and a blowup of the first 50 nodes is shown as inset. [(d)–(f)] The steady-state concentrations for a single realization shown at the network level (ordered the same way as the panel above). Nodes are colored by concentration level. [(g)–(i)] The steady-state concentrations, now averaged over a 100 realizations. We see that the differentiation pattern in the PL network for Range 1 is stable over the space of initial perturbations, whereas the patterns get washed out, either when $\sigma$ is increased or when degree heterogeneities vanish as in the ER network.