A universal route to pattern formation

TL;DR

This paper introduces a new framework for pattern formation that overcomes the limitations of Turing's classical theory, allowing for pattern emergence in systems with any diffusion ratio and many cells, applicable across various natural systems.

Contribution

The authors propose a novel, more general mechanism for pattern formation that does not require large diffusion constant ratios, expanding the applicability of self-organization models.

Findings

Patterns can emerge with any diffusivity ratio.

Pattern formation is possible in systems with sufficiently many cells.

Resulting patterns are robust and can be oscillatory or steady-state.

Abstract

Self-organization, the ability of a system of microscopically interacting entities to shape macroscopically ordered structures, is ubiquitous in Nature. Spatio-temporal patterns are abundantly observed in a large plethora of applications, encompassing different fields and scales. Examples of emerging patterns are the spots and stripes on the coat or skin of animals, the spatial distribution of vegetation in arid areas, the organization of the colonies of insects in host-parasitoid systems and the architecture of large complex ecosystems. Spatial self-organization can be described following the visionary intuition of Alan Turing, who showed how non-linear interactions between slow diffusing activators and fast diffusing inhibitors could induce patterns. The Turing instability, as the mechanism described is universally referred to, was raised to paradigm status in those realms of investigations where microscopic entities are subject to diffusion, from small biological systems to large ecosystems. Requiring a significant ratio of the assigned diffusion constants however is a stringent constraint, which limited the applicability of the theory. Building on the observation that spatial interactions are usually direction biased, and often strongly asymmetric, we here propose a novel framework for the generation of short wavelength patterns which overcomes the limitation inherent in the Turing formulation. In particular, we will prove that patterns can always set in when the system is composed by sufficiently many cells - the units of spatial patchiness - and for virtually any ratio of the diffusivities involved. Macroscopic patterns that follow the onset of the instability are robust and show oscillatory or steady-state behavior.

Figures

Equations

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

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

$$\frac{d}{dt}\left(\begin{array}{c}\delta \text{u}\\ \delta \text{v}\end{array}\right)$$

$$\det \left(\begin{array}{c}{\text{J}}_{\alpha }-{\lambda }_{\alpha }\text{I}\end{array}\right)$$

$${\lambda }_{\alpha }=\frac{1}{2}\left[{\left(\mathrm{tr}{\text{J}}_{\alpha }\right)}_{Re}+\gamma \right]+\frac{1}{2}\left[{\left(\mathrm{tr}{\text{J}}_{\alpha }\right)}_{Im}+\delta \right]\iota$$

$$\gamma$$

$$\delta$$

$$A$$

$$B$$

$$\sqrt{z}=\pm \left(\sqrt{\frac{a+|z|}{2}}+\mathrm{sgn}(b)\sqrt{\frac{-a+|z|}{2}}\iota \right).$$

$$S_2({\Lambda }_{Re}^{(\alpha )})[{\Lambda }_{Im}^{(\alpha )}{]}^2\le -S_1({\Lambda }_{Re}^{(\alpha )}),$$

$$S_1(x)$$

$$S_2(x)$$

$$C_{14}$$

$$C_{13}$$

$$C_{12}$$

$$C_{11}$$

$$C_{10}$$

$$C_{22}$$

$$C_{21}$$

$$C_{20}$$

$$\text{C}=\left(\begin{array}{ccccc}{c}_0& c_{n-1}& \dots & c_2& c_1\\ c_1& c_0& c_{n-1}& \dots & c_2\\ ⋮& c_1& c_0& \ddots & ⋮\\ c_{n-2}& & \ddots & \ddots & c_{n-1}\\ c_{n-1}& c_{n-2}& \dots & c_1& c_0\end{array}\right).$$

Full text

A universal route to pattern formation

Malbor Asllani1, Timoteo Carletti2, Duccio Fanelli3, Philip K. Maini4

1MACSI, Department of Mathematics and Statistics, University of Limerick, Limerick V94 T9PX, Ireland

2Department of Mathematics and naXys, Namur Institute for Complex Systems, University of Namur, rempart de la Vierge 8, B 5000 Namur, Belgium

3Dipartimento di Fisica e Astronomia, Università di Firenze, INFN and CSDC, Via Sansone 1, 50019 Sesto Fiorentino, Firenze, Italy

4Mathematical Institute, University of Oxford, Woodstock Rd, OX2 6GG Oxford, UK

Self-organization, the ability of a system of microscopically interacting entities to shape macroscopically ordered structures, is ubiquitous in Nature. Spatio-temporal patterns are abundantly observed in a large plethora of applications, encompassing different fields and scales. Examples of emerging patterns are the spots and stripes on the coat or skin of animals Kondo and Asai (1995); Murray (1988), the spatial distribution of vegetation in arid areas Klausmeier (1999), the organization of the colonies of insects in host-parasitoid systems Maron and Harrison (1997) and the architecture of large complex ecosystems Rietkerk and van de Koppel (2008). Spatial self-organization can be described following the visionary intuition of Alan Turing, who showed how non-linear interactions between slow diffusing activators and fast diffusing inhibitors could induce patterns Nicolis and Prigogine (1977); Murray (2001). The Turing instability, as the mechanism described is universally referred to, was raised to paradigm status in those realms of investigations where microscopic entities are subject to diffusion, from small biological systems to large ecosystems. Requiring a significant ratio of the assigned diffusion constants however is a stringent constraint, which limited the applicability of the theory. Building on the observation that spatial interactions are usually direction biased, and often strongly asymmetric Shaw et al. (2007); Kim et al. (2005); Murray (2001); Dallon and Othmer (1997); Asllani et al. (2018), we here propose a novel framework for the generation of short wavelength patterns which overcomes the limitation inherent in the Turing formulation. In particular, we will prove that patterns can always set in when the system is composed by sufficiently many cells – the units of spatial patchiness – and for virtually any ratio of the diffusivities involved. Macroscopic patterns that follow the onset of the instability are robust and show oscillatory or steady state behaviour.

In the early $1950^{\prime}s$, Alan Turing laid down, in a seminal paper Turing (1952), the mathematical basis of pattern formation, the discipline that aims at explaining the richness and diversity of forms displayed in Nature. Turing’s idea paved the way for a whole field of investigation and fertilized a cross-disciplinary perspective to yield a universally accepted paradigm of self-organization Ball (1999). The onset of pattern formation on a bound spatial domain originates from the loss of stability of an unpatterned equilibrium. To start with, Turing proposed a minimal model composed of at least two chemicals, hereby termed species. The species were assumed to diffuse across an ensemble of cells, adjacent to each other and organized in a closed ring, as depicted in Figure 1 a). One of the species should trigger its own growth, acting therefore as a self-catalyst. This is opposed by the competing species, which therefore promotes an effective stabilization of the underlying dynamics. The emergence of the ensuing spatial order relies on these contrasting interactions and requires, as an unavoidable constraint, a marked difference (as measured by the ratio) of the diffusion constants associated with the interacting species Murray (2001). This symmetry-breaking mechanism is at the core of a general principle, widely known as the Turing instability Nicolis and Prigogine (1977): small inhomogeneous perturbations from a uniform steady state initiate the instability and individuals, in a quest for space and resources, organise in spatially extended, regular motifs 111The loss of stability of an homogeneous equilibrium, as triggered by an external perturbation, will be referred to as a Turing instability, as an natural extension of the original model discussed in Turing (1952). This original idea was built upon by Meinhardt Gierer and Meinhardt (1972), who proposed the notion of activators and inhibitors, so that the Turing patterning principle could be conceptualized as arising through short-range activation (slow diffusion), long-range inhibition (fast diffusion).

Pattern formation for systems evolving on cellular arrays was further analyzed by Othmer and Scriven Othmer and Scriven (1971) under the assumption of symmetric diffusion. Agregates of cells yield macroscopic tissues which, in general, can be schematized, to a reasonable approximation, by regular lattices Murray (2001). Branching architectures, or coarse-grained models of compartimentalized units, justify invoking the generalized notion of a spatial network. In this case, the nodes stand for individual cellular units, linked via a heterogeneous web of intertangled connections, as exemplified by the network structure. The study of pattern formation for reaction-diffusion systems anchored on symmetric networks was developed in Nakao and Mikhailov (2010). Diffusion instigates a spatial segregation of species, a counterintuitive outcome of Turing analysis, which holds true both in its continuous and discrete (lattice or network based) versions.

Notwithstanding its unquestionable conceptual relevance, it is still unclear if the Turing scheme, in its original formulation, can be successfully applied to describe real life systems. The large ratio of (inhibitor vs. activator) diffusivities can be attained in the laboratory, under controlled operating conditions Castets et al. (1990). An unequivocal proof of the existence of Turing patterns in non-artificial, biological or ecological, contexts is however still lacking, despite several candidate systems that have been thoroughly scrutinized and challenged in this respect Madzvamuse et al. (2015); Levine and Rappel (2005).

Asymmetries in flows, as displayed in real systems via chemical or electrical gradients Shaw et al. (2007), can possibly provide the key to bridge the gap between theory and applications. For example, osmosis Shaw et al. (2007); Kim et al. (2005) in the cell membrane or chemotaxis in the motility of cells Murray (2001); Dallon and Othmer (1997) are examples of asymmetric transport. The case of Dictyostelium is also worth mentioning: this is a multi-celled eukaryotic bacterivore, which develops pseudopodia under externally induced chemotaxis, so triggering a directed biased motion Postma et al. (2003); Dallon and Othmer (1997). Rovinsky & Menzinger Rovinsky and Menzinger (1992, 1993) proved that an unbalance in the directions of the flows of the species can destabilise the spatially homogeneous state of the system. The differential-flow-induced instability yields oscillatory patterns for a sufficiently pronounced degree of asymmetry. This mechanism suffers, however, from the same limitation that applies to the standard Turing theory: diffusion constants need to be significantly different for the pattern to emerge. It is however crucial to realize that asymmetric diffusion can eventually yield a richer pattern formation dynamics, as intuited by Othmer and Scriven in their pioneering work Othmer and Scriven (1971). Turing theory for systems evolving on directed networks was recently cast on rigorous grounds in  Asllani et al. (2014) and asymmetry in the diffusion shown to produce a consistent enlargement of the region of Turing-like patterning.

Starting from these premises, we prove that patterns can develop for virtually any ratio of the diffusivities, larger or smaller than unity, assuming asymmetric transport on (i) an array of sufficiently many cells and (ii) with adequately large values of the nominal diffusion constants (while keeping the sought ratio fixed). To anticipate some of the technical aspects to be addressed in the following, we will show that tuning the diffusivity of just one species – the ratio of the two diffusivities being frozen to a constant strictly different from (but arbitrarily close to) one – amounts to performing a homothetic transformation of the spectrum of the generalized Laplacian operator in the complex plane. The Laplacian eigenvalues consequently move along straight lines, the slope being set by the number, $\Omega$, of cells that define the lattice space. By modulating these latter quantities, one can always get the eigenvalues to protrude into the region of instability. The patterns that follow the onset of the instability are either oscillatory or stationary. A pictorial representation of the proposed scheme is provided in panels c) and d) of Figure 1. Panels a) and b) refer instead to the conventional scenario which assumes undirected transport.

Let us begin by considering a generic two species model of the reaction-diffusion type. The concentration of the two species are, respectively, labelled $u_{i}$ and $v_{i}$, where the index $i$ refers to the hosting cell and $i$ runs from $1$ to $\Omega$. The governing equations can be cast in the form:

[TABLE]

where $D_{u}$ and $D_{v}$ stand for the diffusion constants and $f(\cdot,\cdot)$, $g(\cdot,\cdot)$ are the nonlinear reaction terms that model the local (on site) dynamics of the species. $A_{ij}$ are the binary entries of the adjacency matrix that specifies the topology of the embedding spatial support. For the case at hand (panel c) of Fig. 1), the adjacency matrix $\boldmath{A}=\{A_{ij}\}$ is circulant (i.e. invariant for translation). The transport operator introduced here represents the straightforward generalization of standard diffusion, to the case where the spatial arrangements of the cells is supposed heterogeneous, either in terms of physical links, connecting individual patches of the collection, or in terms of their associated weights. Both are viable strategies for imposing the sought asymmetry which sits at the core of the mechanism upon which we shall hereby elaborate. In short, the hypothesized spatial coupling implements a local balance of incoming and outgoing fluxes, as seen from the observation cell $i$. We introduce the Laplacian operator ${\mathbf{L}}$ with entries $L_{ij}=A_{ij}-k_{i}^{out}\delta_{ij}$ where $\delta_{ij}$ denotes the Kronecker function and $k_{i}^{out}=\sum_{j=1}^{\Omega}A_{ji}$ quantifies the outgoing degree of cell $i$. The contributions that relate to inter-cell couplings in Eqs. (1) can be, respectively, rewritten in the equivalent, more compact form $\sum_{j=1}^{\Omega}D_{u}{L}_{ij}u_{j}$ and $\sigma\sum_{j=1}^{\Omega}D_{u}{L}_{ij}v_{j}$, where $\sigma$ represents the ratio of the diffusion constants $D_{u}$ and $D_{v}$. In the following, we will denote with the symbol $\bm{\mathcal{L}}$ the Laplacian $\mathbf{L}$, modulated by the multiplicative constant $D_{u}$, namely $\bm{\mathcal{L}}\equiv D_{u}\mathbf{L}$. As a key observation for what follows, we emphasize that the magnitude of the eigenvalues of $\bm{\mathcal{L}}$ can be freely controlled by the value of $D_{u}$. Changing $D_{u}$ (while for instance keeping $\sigma$ frozen) implies performing a homothetic transformation of the spectrum of $\mathbf{L}$, as we shall hereafter clarify.

Assume now that the reaction dynamics admits a stable fixed point $({u}^{*},{v}^{*})$, namely that $f({u}^{*},{v}^{*})=g({u}^{*},{v}^{*})=0$. Hence, the spatially extended system (1) possesses a homogeneous equilibrium solution $(\textbf{u}^{*},\textbf{v}^{*})$ which is obtained by replicating on each of the $\Omega$ cells the solution $({u}^{*},{v}^{*})$. This conclusion can be readily derived by noticing that, by definition, $\sigma\sum_{j=1}^{\Omega}\mathcal{L}_{ij}K=0$, for any constant $K$ (recall that the system is hosted on a lattice with periodic boundary conditions). The homogeneous fixed point can, in principle, become unstable upon injection of a tiny heterogeneous perturbation, as in the spirit of the original Turing mechanism. The conditions for the emergence of the instability are determined by a linear stability analysis that follows the procedure outlined in  Asllani et al. (2014) and that we develop in the annexed Methods. The calculation yields the compact inequality (3) which can be graphically illustrated in the complex plane \textbf{z}=\big{(}\Lambda_{Re},\Lambda_{Im}\big{)}, where the eigenvalues $\Lambda^{(\alpha)}$ of the Laplacian operator $\mathbf{L}$ reside (see Methods). In fact, inequality (3) enables one to delimit a model-dependent region of instability, which is depicted in panel a) of Fig. 2. When drawing the domain of interest, we assumed an abstract setting, without insisting on the specific details that stem from a particular reaction model. For what will follow, it is only important to appreciate that condition (3) makes it possible to define a portion of the parameter space, by construction symmetric with respect to the horizontal (real) axis, that is eventually associated with the onset of the instability. More specifically, if a subset of the spectrum of the Laplacian $\bm{\mathcal{L}}$ falls inside the region outlined above, the instability can take place (red stars in panel a) of Fig. 2). This conclusion is general and independent of the reaction scheme employed. Notice that when the two regions merge and incorporate a finite part of the real axis, the outbreak of the instability is also possible on a symmetric support (when i.e. the eigenvalues of the Laplacian operator are real).

We are now in a position to elaborate on the universal mechanism which drives a reaction system unstable, when bound to performing asymmetric diffusion on a discrete collection of lattice sites. As we shall realize, considering a finite, although large, ensemble of mutually connected cells is one of the key ingredients that instigates the instability in an activator-inhibitor system for virtually any reaction parameters and any ratio of the diffusion coefficients (including $\sigma<1$ and except for the zero measure, limiting condition $\sigma=1$). We begin by observing that the Laplacian matrix $\mathbf{L}$, associated with a closed directed ring as assumed in Fig. 1 c), is a circulant matrix. This simplified geometrical arrangement is solely assumed for illustrative purposes. Our conclusions hold in general and include those systems that undergo asymmetric diffusion on lattices in arbitrary dimensions, subject to periodic boundary conditions. The eigenvalues of the Laplacian $\mathbf{L}$ (see Methods) are complex and fall on the unitary circle centered at $(-1,0)$ – empty stars in Fig. 2, panel a). When making $\Omega$ larger, one progressively reduces the spectral gap, the relative distance between the two Laplacian eigenvalues that display the largest real parts. Recalling that the largest eigenvalue of the Laplacian is by definition zero, this implies that the second eigenvalue of $\mathbf{L}$ (ranked in descending order, with respect to the value of their associated real parts) tends to approach the origin of the complex plane, when $\Omega$ is increased 222The eigenvalues comes in conjugate pairs. In what it follows we shall refer to the eigenvalue in the pair that displays positive imaginary part.. This, in turn, implies that we can control at will the inclination (a measure complementary to the slope) $\theta$ of the line (dashed in the figure) that connects the second largest eigenvalues to the origin of the complex plane. Similar considerations apply to the other eigenvalues that follow in the ranking. In Fig. 2) a) we label by $\theta^{*}$ the critical inclination of the solid line that intersects tangentially the domain of instability traced according to inequality (3). Clearly, by construction, a symmetric line always exists with inclination $-\theta^{*}$, that is tangent to the region of instability in the lower complex semi-plane. Assume now that $\Omega$ is sufficiently large so that $\theta<\theta^{*}$. Then recall that for the instability to emerge, at least one eigenvalue of the rescaled Laplacian $\bm{\mathcal{L}}$ has to protrude into the region where the instability is bound to occur. On the other hand, the eigenvalues of $\bm{\mathcal{L}}$ are a homothetic transformation, with centre [math] and ratio $D_{u}$, of the eigenvalues $\mathbf{L}$. In other words, we can force the second eigenvalue of $\mathbf{L}$ to move along the line to which it belongs by modulating the diffusion constant $D_{u}$, while keeping $\sigma$ constant, and so invade the region of the instability (filled red star). If $\theta>\theta^{*}$, then the eigenvalue is instead confined to the stable portion of the complex plane for each choice of the scaling factor $D_{u}$. In summary, the system can be triggered unstable by simultaneously acting on two independents “knobs”: first, by making $\Omega$ sufficiently large, we force a subset of eigenvalues into the vicinity of the origin, thus making their associated $\theta$ smaller that the critical amount $\theta^{*}$. Secondly, by acting on $D_{u}$ (and consequently on $D_{v}$ so as to maintain $\sigma$ unchanged), we make the eigenvalues slide on their corresponding lines, until the instability threshold is eventually breached. Interestingly, the instability involves long wavelengths (hence yielding macroscopic patterns) since, by construction, the real part of the eigenvalues associated with the modes triggered to grow approaches zero, when $\Omega$ increases.

In Fig. 2 b) we present results for a specific case study, the Brusselator model, often invoked in the literature as a paradigm nonlinear reaction scheme for studying self-organised phenomena, synchronisation, Turing patterns and oscillation death. We find the region in the parameter space $(b,c)$ for which the homogeneous steady state is unstable for $\sigma=1.4$. The Brusselator may undergo a conventional Turing instability, in the portion of the plane that is colored in blue. The asymmetry driven instability (for $\Omega=1000$ cells and assuming $D_{u}=100$) is found on a considerably larger domain, which can be made arbitrarily large by further modulating $\Omega$ and $D_{u}$.

In Fig. 3, panels b), d) and f), we show patterns found when integrating the Brusselator model for different parameter values. Panels a), c) and e) display the real and imaginary parts of the eigenvalues $\lambda_{\alpha}$ (the dispersion relation, see Methods) against $-\Lambda^{(\alpha)}_{Re}$. When the imaginary component of $\lambda_{\alpha}$ is small as compared to the corresponding real part, steady state patterns are found to emerge. Otherwise, the patterns are oscillatory in nature. This is a rule of the thumb, which seems to implement a necessary criterion. The mechanism of pattern selection is in fact heavily influenced by the nonlinearities, which become prominent beyond the initial stages of evolution, hence limiting the predictive ability of the linear stability analysis. It is worth stressing that the patterns displayed in Fig. 3 occur for a choice of parameters for which the classical Turing instability is not permitted and for both $\sigma$ larger and smaller than one. In particular, one can analytically show that the asymmetry driven instability, that we here introduced and characterized, can develop for any $\sigma=1\pm\epsilon$ with $0<\epsilon\ll 1$. This conclusion holds in general, as can be proved rigorously (see Methods).

To summarise, we have developed a general theory of pattern formation for the case of asymmetric transport of constituents in a large assembly of adjacent cells. Macroscopic oscillatory or steady state patterns in an activator-inhibitor system are now found for virtually any ratio of diffusivities, larger or smaller than one, and any choice of the reaction parameters. The applications of this newly proposed route to pattern formation are multiple and include all those settings, from biology to ecology passing through neuroscience Shaw et al. (2007); Pringle et al. (2011); Kim et al. (2005); Sporns (2010); Dallon and Othmer (1997); Postma et al. (2003), where asymmetric flows are reported to occur.

Acknowledgments

The work of M.A. was supported by a FRS-FNRS Postdoctoral Fellowship.

**Methods

**

**The dispersion relation and the conditions for instability

**Linearising the dynamics of system (1) around the fixed point solution $(\textbf{u}^{*},\textbf{v}^{*})$ we obtain:

[TABLE]

where $\delta\textbf{u}$, $\delta\textbf{v}$ stands for the perturbation vectors. $\bm{\mathcal{J}}=\begin{pmatrix}f_{u}\bm{\mathcal{I}}&f_{v}\bm{\mathcal{I}}\\ g_{u}\bm{\mathcal{I}}&g_{v}\bm{\mathcal{I}}\end{pmatrix}$ is the Jacobian matrix, evaluated at the equilibrium point, stemming from the reaction terms and $\bm{\mathcal{D}}=\begin{pmatrix}D_{u}\mathbf{L}&\bm{\mathcal{O}}\\ \bm{\mathcal{O}}&D_{v}\mathbf{L}\end{pmatrix}$, where $\bm{\mathcal{O}}$ is the $\Omega\times\Omega$ zero-valued matrix. Introduce the basis formed by the eigenvectors $\Phi_{i}^{(\alpha)}$, with $\alpha=1,\dots,\Omega$, of the Laplacian operator $\mathbf{L}$. We have $\sum_{j}L_{ij}\Phi_{j}^{(\alpha)}=\Lambda^{(\alpha)}\Phi_{i}^{(\alpha)}$ where $\Lambda^{(\alpha)}$ identifies the eigenvalues of $\mathbf{L}$. The latter operator is asymmetric, as is the matrix $\mathbf{A}$, and thus $\Lambda^{(\alpha)}$ are, in principle, complex. Further, the eigenvectors form an orthonormal basis, for the case at hand. Hence, to solve the above linear system, we expand the perturbation in terms of the basis of the eigenvectors, i.e. $\delta u_{i}=\sum_{\alpha=1}^{\Omega}b_{\alpha}\Phi_{i}^{(\alpha)}$ and $\delta v_{i}=\sum_{\alpha=1}^{\Omega}c_{\alpha}\Phi_{i}^{(\alpha)}$. At this point it is straightforward to show that the $2\Omega\times 2\Omega$ system reduces to a $2\times 2$ eigenvalue problem, for each choice of the scalar index $\alpha=1,\dots,\Omega$:

[TABLE]

where $\textbf{J}_{\alpha}$ is the $2\times 2$ modified Jacobian, i.e. $\textbf{J}_{\alpha}\equiv\textbf{J}+\textbf{D}\bm{\Lambda}^{(\alpha)}$ with $\textbf{D}=diag(D_{u},D_{v})$, $\bm{\Lambda}^{(\alpha)}=diag(\Lambda^{(\alpha)},\Lambda^{(\alpha)})$. The steady state is unstable to small heterogeneous perturbations, if $\lambda_{\alpha}$ has a positive real part over a finite range of modes. The dispersion relation (the largest real part of $\lambda_{\alpha}$, for any given $\alpha$) can be readily computed from:

[TABLE]

where

[TABLE]

and

[TABLE]

Here $\mathrm{sgn}(\cdot)$ is the sign function and $f_{Re}$, $f_{Im}$ indicate, respectively, the real and imaginary parts of the operator $f$. To further manipulate Eq. (2), we make use of the definition of a square root of a complex number. Take, $z=a+b\iota$ where $\iota=\sqrt{-1}$ is the imaginary unit, then

[TABLE]

The instability sets in when $|\left(\mathrm{tr}\textbf{J}_{\alpha}\right)_{Re}|\leq\gamma$, a condition that translates into the following inequality Asllani et al. (2014):

[TABLE]

where $(\Lambda_{Re},\Lambda_{Im})$ span the complex plane where the Laplacian eigenvalues reside. In the above relation $S_{1}$, $S_{2}$ are polynomials which take the following explicit form:

[TABLE]

The constants here are given by:

[TABLE]

and

[TABLE]

The system displays a generalized Turing instability, determined by the asymmetric nature of the imposed coupling if, after the homothetic transformation, the eigenvalues $D_{u}(\Lambda^{(\alpha)}_{Re},\Lambda^{(\alpha)}_{Im})$ of the Laplacian operator $\bm{\mathcal{L}}$ enter the region of the complex plane $(\Lambda_{Re},\Lambda_{Im})$ that is delimited by inequality (3).

**Spectral properties of circulant matrices

**An $n\times n$ matrix $C$ is circulant if takes the form

[TABLE]

The circulant matrix C is fully specified by its first column, $\textbf{c}=(c_{0},c_{1},\dots c_{n-1})$. The other columns of C are generated as cyclic permutations of the vector c with offset equal to the column index. The normalized eigenvectors of a circulant matrix are given by $\phi_{i}=\dfrac{1}{\sqrt{n}}\left(1,\omega_{i},\omega_{i}^{2},\dots,\omega_{i}^{n-1}\right)$ where $\omega_{i}=\exp(2\pi\iota i/n)$ is the $n-$th root of the unity and the eigenvalues are $\lambda_{i}=c_{0}+c_{n-1}\omega_{i}+c_{n-2}\omega_{i}^{2}+\dots+c_{1}\omega_{i}^{n-1}$ where $i=0,1,\dots,n-1$.