A universal route to pattern formation
Malbor Asllani, Timoteo Carletti, Duccio Fanelli, Philip K. Maini
AA universal route to pattern formation
Malbor Asllani , Timoteo Carletti , Duccio Fanelli , Philip K. Maini MACSI, Department of Mathematics and Statistics,University of Limerick, Limerick V94 T9PX, Ireland Department of Mathematics and naXys, Namur Institute for Complex Systems,University of Namur, rempart de la Vierge 8, B 5000 Namur, Belgium Dipartimento di Fisica e Astronomia, Universit`a di Firenze,INFN and CSDC, Via Sansone 1, 50019 Sesto Fiorentino, Firenze, Italy and Mathematical Institute, University of Oxford, Woodstock Rd, OX2 6GG Oxford, UK
Self-organization, the ability of a system of mi-croscopically interacting entities to shape macro-scopically ordered structures, is ubiquitous in Na-ture. Spatio-temporal patterns are abundantlyobserved in a large plethora of applications, en-compassing different fields and scales. Examplesof emerging patterns are the spots and stripes onthe coat or skin of animals [1, 2], the spatial distri-bution of vegetation in arid areas [3], the organi-zation of the colonies of insects in host-parasitoidsystems [4] and the architecture of large complexecosystems [5]. Spatial self-organization can bedescribed following the visionary intuition of AlanTuring, who showed how non-linear interactionsbetween slow diffusing activators and fast diffus-ing inhibitors could induce patterns [6, 7]. TheTuring instability, as the mechanism described isuniversally referred to, was raised to paradigmstatus in those realms of investigations where mi-croscopic entities are subject to diffusion, fromsmall biological systems to large ecosystems. Re-quiring a significant ratio of the assigned diffusionconstants however is a stringent constraint, whichlimited the applicability of the theory. Buildingon the observation that spatial interactions areusually direction biased, and often strongly asym-metric [7–11], we here propose a novel frameworkfor the generation of short wavelength patternswhich overcomes the limitation inherent in theTuring formulation. In particular, we will provethat patterns can always set in when the system iscomposed by sufficiently many cells – the units ofspatial patchiness – and for virtually any ratio ofthe diffusivities involved. Macroscopic patternsthat follow the onset of the instability are robustand show oscillatory or steady state behaviour.
In the early 1950 s , Alan Turing laid down, in a seminalpaper [12], the mathematical basis of pattern formation,the discipline that aims at explaining the richness and di-versity of forms displayed in Nature. Turing’s idea pavedthe way for a whole field of investigation and fertilizeda cross-disciplinary perspective to yield a universally ac-cepted paradigm of self-organization [13]. The onset ofpattern formation on a bound spatial domain originatesfrom the loss of stability of an unpatterned equilibrium.To start with, Turing proposed a minimal model com- posed of at least two chemicals, hereby termed species.The species were assumed to diffuse across an ensembleof cells, adjacent to each other and organized in a closedring, as depicted in Figure 1 a). One of the species shouldtrigger its own growth, acting therefore as a self-catalyst.This is opposed by the competing species, which there-fore promotes an effective stabilization of the underlyingdynamics. The emergence of the ensuing spatial orderrelies on these contrasting interactions and requires, asan unavoidable constraint, a marked difference (as mea-sured by the ratio) of the diffusion constants associatedwith the interacting species [7]. This symmetry-breakingmechanism is at the core of a general principle, widelyknown as the Turing instability [6]: small inhomoge-neous perturbations from a uniform steady state initi-ate the instability and individuals, in a quest for spaceand resources, organise in spatially extended, regular mo-tifs [14]. This original idea was built upon by Mein-hardt [15], who proposed the notion of activators and in-hibitors, so that the Turing patterning principle could beconceptualized as arising through short-range activation(slow diffusion), long-range inhibition (fast diffusion).Pattern formation for systems evolving on cellular ar-rays was further analyzed by Othmer and Scriven [16]under the assumption of symmetric diffusion. Agregatesof cells yield macroscopic tissues which, in general, can beschematized, to a reasonable approximation, by regularlattices [7]. Branching architectures, or coarse-grainedmodels of compartimentalized units, justify invoking thegeneralized notion of a spatial network. In this case,the nodes stand for individual cellular units, linked viaa heterogeneous web of intertangled connections, as ex-emplified by the network structure. The study of pat-tern formation for reaction-diffusion systems anchoredon symmetric networks was developed in [17]. Diffusioninstigates a spatial segregation of species, a counterintu-itive outcome of Turing analysis, which holds true bothin its continuous and discrete (lattice or network based)versions.Notwithstanding its unquestionable conceptual rele-vance, it is still unclear if the Turing scheme, in its orig-inal formulation, can be successfully applied to describereal life systems. The large ratio of (inhibitor vs. acti-vator) diffusivities can be attained in the laboratory, un-der controlled operating conditions [18]. An unequivocalproof of the existence of Turing patterns in non-artificial, a r X i v : . [ n li n . PS ] J un biological or ecological, contexts is however still lacking,despite several candidate systems that have been thor-oughly scrutinized and challenged in this respect [19, 20].Asymmetries in flows, as displayed in real systems viachemical or electrical gradients [8], can possibly providethe key to bridge the gap between theory and applica-tions. For example, osmosis [8, 9] in the cell membrane orchemotaxis in the motility of cells [7, 10] are examples ofasymmetric transport. The case of Dictyostelium is alsoworth mentioning: this is a multi-celled eukaryotic bac-terivore, which develops pseudopodia under externallyinduced chemotaxis, so triggering a directed biased mo-tion [10, 21]. Rovinsky & Menzinger [22, 23] proved thatan unbalance in the directions of the flows of the speciescan destabilise the spatially homogeneous state of thesystem. The differential-flow-induced instability yieldsoscillatory patterns for a sufficiently pronounced degreeof asymmetry. This mechanism suffers, however, fromthe same limitation that applies to the standard Turingtheory: diffusion constants need to be significantly dif-ferent for the pattern to emerge. It is however crucial torealize that asymmetric diffusion can eventually yield aricher pattern formation dynamics, as intuited by Oth-mer and Scriven in their pioneering work [16]. Turingtheory for systems evolving on directed networks was re-cently cast on rigorous grounds in [24] and asymmetry inthe diffusion shown to produce a consistent enlargementof the region of Turing-like patterning.Starting from these premises, we prove that patternscan develop for virtually any ratio of the diffusivities,larger or smaller than unity, assuming asymmetric trans-port on (i) an array of sufficiently many cells and (ii) withadequately large values of the nominal diffusion constants(while keeping the sought ratio fixed). To anticipate someof 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 con-stant strictly different from (but arbitrarily close to) one– amounts to performing a homothetic transformationof the spectrum of the generalized Laplacian operatorin the complex plane. The Laplacian eigenvalues conse-quently move along straight lines, the slope being set bythe number, Ω, of cells that define the lattice space. Bymodulating these latter quantities, one can always getthe eigenvalues to protrude into the region of instability.The patterns that follow the onset of the instability areeither oscillatory or stationary. A pictorial representa-tion of the proposed scheme is provided in panels c) andd) of Figure 1. Panels a) and b) refer instead to the con-ventional scenario which assumes undirected transport.Let us begin by considering a generic two species modelof the reaction-diffusion type. The concentration of thetwo species are, respectively, labelled u i and v i , wherethe index i refers to the hosting cell and i runs from 1 to Ω. The governing equations can be cast in the form: du i dt = f ( u i , v i ) + D u Ω X j =1 A ij u j − Ω X j =1 A ji u i dv i dt = g ( u i , v i ) + D v Ω X j =1 A ij v j − Ω X j =1 A ji v i (1)where D u and D v stand for the diffusion constants and f ( · , · ), g ( · , · ) are the nonlinear reaction terms that modelthe local (on site) dynamics of the species. A ij arethe binary entries of the adjacency matrix that speci-fies the topology of the embedding spatial support. Forthe case at hand (panel c) of Fig. 1), the adjacencymatrix A = { A ij } is circulant (i.e. invariant for transla-tion). The transport operator introduced here representsthe straightforward generalization of standard diffusion,to the case where the spatial arrangements of the cells issupposed heterogeneous, either in terms of physical links,connecting individual patches of the collection, or interms of their associated weights. Both are viable strate-gies for imposing the sought asymmetry which sits at thecore of the mechanism upon which we shall hereby elab-orate. In short, the hypothesized spatial coupling imple-ments a local balance of incoming and outgoing fluxes, asseen from the observation cell i . We introduce the Lapla-cian operator L with entries L ij = A ij − k outi δ ij where δ ij denotes the Kronecker function and k outi = P Ω j =1 A ji quantifies the outgoing degree of cell i . The contribu-tions that relate to inter-cell couplings in Eqs. (1) canbe, respectively, rewritten in the equivalent, more com-pact form P Ω j =1 D u L ij u j and σ P Ω j =1 D u L ij v j , where σ represents the ratio of the diffusion constants D u and D v . In the following, we will denote with the symbol L the Laplacian L , modulated by the multiplicative con-stant D u , namely L ≡ D u L . As a key observation forwhat follows, we emphasize that the magnitude of theeigenvalues of L can be freely controlled by the value of D u . Changing D u (while for instance keeping σ frozen)implies performing a homothetic transformation of thespectrum of L , as we shall hereafter clarify.Assume now that the reaction dynamics admits astable fixed point ( u ∗ , v ∗ ), namely that f ( u ∗ , v ∗ ) = g ( u ∗ , v ∗ ) = 0. Hence, the spatially extended system (1)possesses a homogeneous equilibrium solution ( u ∗ , v ∗ )which is obtained by replicating on each of the Ω cellsthe solution ( u ∗ , v ∗ ). This conclusion can be readily de-rived by noticing that, by definition, σ P Ω j =1 L ij K = 0,for any constant K (recall that the system is hosted on alattice with periodic boundary conditions). The homoge-neous fixed point can, in principle, become unstable uponinjection of a tiny heterogeneous perturbation, as in thespirit of the original Turing mechanism. The conditionsfor the emergence of the instability are determined by alinear stability analysis that follows the procedure out-lined in [24] and that we develop in the annexed Meth-ods. The calculation yields the compact inequality (3) Turing pattern D v ≫ D u a)c) D v ∼ D u b) activatorinhibitor (-) (-) (+)(+) diffusion diffusiondiffusiondiffusion cell i cell i +1 cell i − d) activatorinhibitor (-) (-) (+)(+) diffusiondiffusion cell i cell i +1 cell i − diffusiondiffusion Asymmetry pattern
FIG. 1.
Conventional Turing instability vs. the asymmetry-driven model: a schematic representation. (a)
InTuring’s original model the onset of pattern formation is studied for a two species model reacting and diffusing on a collection ofcells, arranged so as to form a 1 D ring. (b) Turing instability requires breaking the symmetry among the species. In particular,feedback loops (positive for the activators and negative for the inhibitors) are proposed. Further, the inhibitor should relocatein space faster than the activator, D v (cid:29) D u . The diffusion between neighboring cells is assumed symmetric. (c) In theasymmetry-induced instability instead the system is made up of a larger number of cells and the diffusion is asymmetric, asschematised by the counterclockwise arrow. (d)
The instability is triggered by increasing the number of cells, for virtually anyratio of the diffusion constants, provided the latter take sufficiently large values. which can be graphically illustrated in the complex plane z = (cid:0) Λ Re , Λ Im (cid:1) , where the eigenvalues Λ ( α ) of the Lapla-cian operator L reside (see Methods). In fact, inequality(3) enables one to delimit a model-dependent region ofinstability, which is depicted in panel a) of Fig. 2. Whendrawing the domain of interest, we assumed an abstractsetting, without insisting on the specific details that stemfrom a particular reaction model. For what will follow, itis only important to appreciate that condition (3) makesit 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 ofthe instability. More specifically, if a subset of the spec-trum of the Laplacian L falls inside the region outlinedabove, the instability can take place (red stars in panela) of Fig. 2). This conclusion is general and independentof the reaction scheme employed. Notice that when thetwo regions merge and incorporate a finite part of the real axis, the outbreak of the instability is also possibleon a symmetric support (when i.e. the eigenvalues of theLaplacian operator are real).We are now in a position to elaborate on the univer-sal mechanism which drives a reaction system unstable,when bound to performing asymmetric diffusion on a dis-crete collection of lattice sites. As we shall realize, con-sidering a finite, although large, ensemble of mutuallyconnected cells is one of the key ingredients that insti-gates the instability in an activator-inhibitor system forvirtually any reaction parameters and any ratio of thediffusion coefficients (including σ < σ = 1). We begin byobserving that the Laplacian matrix L , associated witha closed directed ring as assumed in Fig. 1 c), is a circu-lant matrix. This simplified geometrical arrangement issolely assumed for illustrative purposes. Our conclusionshold in general and include those systems that undergo a) Asymptote b) Turing classicalregionAsymmetric diffusionextended region bc θ ∗ 𝛬 Re 𝛬 Im FIG. 2.
On the mechanism of pattern formation with asymmetric diffusion. a)
The region colored in green denotes thedomain of the complex plane z = (cid:0) Λ Re , Λ Im (cid:1) where the instability can eventually take place. This is a pictorial representationof a general situation that is always found, irrespectively of the specific choice of the reaction model. The blue empty symbolsstand for the spectrum of the Laplacian operator L obtained from a directed lattice of the type depicted in Fig. 1 c). Theeigenvalues are distributed on the unitary circle, as discussed in the Methods section. When increasing Ω, the number of cellsthat compose the examined lattice, the Laplacian spectrum displays eigenvalues with progressively larger (but still negative)component Λ Re and still lying on the circle. We denote by θ the inclination of the line (dot dashed in the figure) that connectsthe selected eigenvalue (empty cross) to the origin of the complex plane. This latter eigenvalue can be freely moved along theaforementioned line, by modulating the diffusion coefficient D u while keeping the ratio σ fixed. Stated differently, the spectrumof the matrix L is an homothetic transformation (with centre 0 and ratio D u ) of the spectrum of L . If θ is sufficiently small,i.e. if the dashed line is steeper than the solid one (the instability threshold), it is always possible (by increasing D u ) to movethe associated eigenvalue inside the region (green) of the instability (see the filled red star). If the inclination θ is larger thanthe critical one, θ ∗ , the eigenvalues that slide on the corresponding dashed line are permanently confined outside the domain ofinstability (filled blue star). The vertical dashed line is an asymptote and sets the leftmost boundary of the instability domain,as described in the annexed Methods. b) To provide a quantitative illustration of the method, we consider the Brusselatormodel, f ( u, v ) = 1 − ( b + 1) u + cu v and g ( u, v ) = bu − cu v , where b and c are free control parameters. In the plane ( b, c ), weisolate the domain where the conventional Turing instability takes place (small blue domain), for σ = 1 .
4. The region shadedin red identifies the domain of instability that is found when the system is evolved on a directed lattice made of Ω = 1000 cellsand assuming D u = 100. The ratio σ is kept constant to the reference value 1 . asymmetric diffusion on lattices in arbitrary dimensions,subject to periodic boundary conditions. The eigenvaluesof the Laplacian L (see Methods) are complex and fallon the unitary circle centered at ( − ,
0) – empty starsin Fig. 2, panel a). When making Ω larger, one pro-gressively reduces the spectral gap, the relative distancebetween the two Laplacian eigenvalues that display thelargest real parts. Recalling that the largest eigenvalueof the Laplacian is by definition zero, this implies thatthe second eigenvalue of 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Ω is increased [25]. This, in turn, implies that we cancontrol at will the inclination (a measure complementaryto the slope) θ of the line (dashed in the figure) that con-nects the second largest eigenvalues to the origin of thecomplex plane. Similar considerations apply to the othereigenvalues that follow in the ranking. In Fig. 2) a) welabel by θ ∗ the critical inclination of the solid line thatintersects tangentially the domain of instability traced according to inequality (3). Clearly, by construction, asymmetric line always exists with inclination − θ ∗ , thatis tangent to the region of instability in the lower com-plex semi-plane. Assume now that Ω is sufficiently largeso that θ < θ ∗ . Then recall that for the instability toemerge, at least one eigenvalue of the rescaled Laplacian L has to protrude into the region where the instabilityis bound to occur. On the other hand, the eigenvaluesof L are a homothetic transformation, with centre 0 andratio D u , of the eigenvalues L . In other words, we canforce the second eigenvalue of L to move along the lineto which it belongs by modulating the diffusion constant D u , while keeping σ constant, and so invade the regionof the instability (filled red star). If θ > θ ∗ , then theeigenvalue is instead confined to the stable portion of thecomplex plane for each choice of the scaling factor D u . Insummary, the system can be triggered unstable by simul-taneously acting on two independents “knobs”: first, bymaking Ω sufficiently large, we force a subset of eigenval-ues into the vicinity of the origin, thus making their asso- - Λ ( α ) Re λ α -250-200-150-100-50050100 a) T ime u i ( t ) b) T ime u i ( t ) T ime u i ( t ) - Λ ( α ) Re λ α -300-250-200-150-100-50050100 e)f) - Λ ( α ) Re λ α -15-10-505 c)d) FIG. 3.
Patterns triggered in the Brusselator model by asymmetric diffusion on a D ring. Left column: real (bluestars) and imaginary (red circles) dispersion relation λ α a) and the associated pattern evolution b) . The ensuing pattern isoscillatory and organizes as a travelling wave. The parameters are set to b = 8, c = 10, D u = 100, and σ = 1 .
4. Center column:dispersion relation c) and pattern evolution d) yielding a steady state pattern, when b = 50, c = 62, D u = 10 and σ = 1 . e) and pattern evolution f ) yielding a traveling wave with b = 8, c = 10, D u = 140 and σ = 1 / . ’ .
71. In the insets in the upper panels, we zoom on the (real and imaginary) dispersion relations focusing on theportion of the curve where the instability takes place. In all cases, Ω = 100. Notice that, for the Brusselator model, oscillatoryand stationary patterns are numerically detected for σ >
1, while for σ < σ < ciated θ smaller that the critical amount θ ∗ . Secondly, byacting on D u (and consequently on D v so as to maintain σ unchanged), we make the eigenvalues slide on their cor-responding lines, until the instability threshold is eventu-ally breached. Interestingly, the instability involves longwavelengths (hence yielding macroscopic patterns) since,by construction, the real part of the eigenvalues associ-ated with the modes triggered to grow approaches zero,when Ω increases.In Fig. 2 b) we present results for a specific case study,the Brusselator model, often invoked in the literature asa paradigm nonlinear reaction scheme for studying self-organised phenomena, synchronisation, Turing patternsand oscillation death. We find the region in the param-eter space ( b, c ) for which the homogeneous steady stateis unstable for σ = 1 .
4. The Brusselator may undergoa conventional Turing instability, in the portion of theplane that is colored in blue. The asymmetry driven in-stability (for Ω = 1000 cells and assuming D u = 100)is found on a considerably larger domain, which can bemade arbitrarily large by further modulating Ω and D u .In Fig. 3, panels b), d) and f), we show patterns foundwhen integrating the Brusselator model for different pa-rameter values. Panels a), c) and e) display the real andimaginary parts of the eigenvalues λ α (the dispersion re- lation, see Methods) against − Λ ( α ) Re . When the imaginarycomponent of λ α is small as compared to the correspond-ing real part, steady state patterns are found to emerge.Otherwise, the patterns are oscillatory in nature. This isa rule of the thumb, which seems to implement a neces-sary criterion. The mechanism of pattern selection is infact heavily influenced by the nonlinearities, which be-come prominent beyond the initial stages of evolution,hence limiting the predictive ability of the linear stabilityanalysis. It is worth stressing that the patterns displayedin Fig. 3 occur for a choice of parameters for which theclassical Turing instability is not permitted and for both σ larger and smaller than one. In particular, one cananalytically show that the asymmetry driven instability,that we here introduced and characterized, can developfor any σ = 1 ± (cid:15) with 0 < (cid:15) (cid:28)
1. This conclusion holdsin general, as can be proved rigorously (see Methods).To summarise, we have developed a general theory ofpattern formation for the case of asymmetric transport ofconstituents in a large assembly of adjacent cells. Macro-scopic oscillatory or steady state patterns in an activator-inhibitor system are now found for virtually any ratio ofdiffusivities, larger or smaller than one, and any choice ofthe reaction parameters. The applications of this newlyproposed route to pattern formation are multiple and in-clude all those settings, from biology to ecology passing through neuroscience [8–10, 21, 26, 27], where asymmet-ric flows are reported to occur. [1] S. Kondo and R. Asai, Nature , 765 (1995).[2] J. D. Murray, Sci. Am. , 80 (1988).[3] C. Klausmeier, Science , 1826 (1999).[4] J. L. Maron and S. Harrison, Science , 1619 (1997).[5] M. Rietkerk and J. van de Koppel, Trends Ecol. Evol. , 169 (2008).[6] G. Nicolis and I. Prigogine, Self-organization in nonequi-librium systems: From dissipative structures to orderthrough fluctuations (J. Wiley & Sons, 1977).[7] J. D. Murray,
Mathematical biology II: Spatial models andbiomedical applications (Springer-Verlag, 2001).[8] R. S. Shaw, N. Packard, M. Schroter, and H. L. Swinney,PNAS , 9580 (2007).[9] K. S. Kim, I. S. Davis, P. A. Macpherson, T. J. Pedley,and A. E. Hill, Proc. R. Soc. A , 273 (2005).[10] J. C. Dallon and H. G. Othmer, Phil. Trans. R. Soc.Lond. B , 391 (1997).[11] M. Asllani, R. Lambiotte, and T. Carletti, Sci. Adv. ,eaau9403 (2018).[12] A. M. Turing, Phil. Trans. R. Soc. B , 37 (1952).[13] P. Ball, The self-made tapestry: Pattern formation inNature , 1st ed. (Oxford University Press, 1999).[14] The loss of stability of an homogeneous equilibrium, astriggered by an external perturbation, will be referredto as a Turing instability, as an natural extension of theoriginal model discussed in [12].[15] A. Gierer and H. Meinhardt, Kybernetik , 30 (1972).[16] H. G. Othmer and L. E. Scriven, J. theor. Biol. , 507(1971).[17] H. Nakao and A. S. Mikhailov, Nat. Phys. , 544 (2010).[18] V. Castets, E. Dulos, J. Boissonade, and P. De Kepper,Phys. Rev. Lett. , 2953 (1990).[19] A. Madzvamuse, H. S. Ndakwo, and R. Barreira, J.Math. Biol. , 709 (2015).[20] H. Levine and W.-J. Rappel, Phys. Rev. E , 061912(2005).[21] M. Postma, J. Roelofs, J. Goedhart, T. W. J. Gadella,A. J. W. G. Visser, and P. J. M. V. Haastert, Mol. Biol.Cell , 5019 (2003).[22] A. B. Rovinsky and M. Menzinger, Phys. Rev. Lett. ,1193 (1992).[23] A. B. Rovinsky and M. Menzinger, Phys. Rev. Lett. ,778 (1993).[24] M. Asllani, J. D. Challenger, F. S. Pavone, L. Sacconi,and D. Fanelli, Nature Commun. , 4517 (2014).[25] The eigenvalues comes in conjugate pairs. In what it fol-lows we shall refer to the eigenvalue in the pair that dis-plays positive imaginary part.[26] J. M. Pringle, A. M. H. Blakeslee, J. E. Byers, andJ. Roman, PNAS , 15288 (2011).[27] O. Sporns, Networks of the Brain (MIT Press, 2010).
ACKNOWLEDGMENTS
The work of M.A. was supported by a FRS-FNRS Post-doctoral Fellowship.
Methods
The dispersion relation and the conditions for in-stability
Linearising the dynamics of system (1) around the fixedpoint solution ( u ∗ , v ∗ ) we obtain: ddt (cid:18) δ u δ v (cid:19) = (cid:18) f u I + D u L f v I g u I g v + D v L (cid:19) (cid:18) δ u δ v (cid:19) = (cid:0) J + D (cid:1) (cid:18) δ u δ v (cid:19) where δ u , δ v stands for the perturbation vectors. J = (cid:18) f u I f v I g u I g v I (cid:19) is the Jacobian matrix, evaluated at theequilibrium point, stemming from the reaction terms and D = (cid:18) D u L OO D v L (cid:19) , where O is the Ω × Ω zero-valuedmatrix. Introduce the basis formed by the eigenvectorsΦ ( α ) i , with α = 1 , . . . , Ω, of the Laplacian operator L .We have P j L ij Φ ( α ) j = Λ ( α ) Φ ( α ) i where Λ ( α ) identifiesthe eigenvalues of L . The latter operator is asymmet-ric, as is the matrix A , and thus Λ ( α ) are, in principle,complex. Further, the eigenvectors form an orthonormalbasis, for the case at hand. Hence, to solve the abovelinear system, we expand the perturbation in terms ofthe basis of the eigenvectors, i.e. δu i = P Ω α =1 b α Φ ( α ) i and δv i = P Ω α =1 c α Φ ( α ) i . At this point it is straightfor-ward to show that the 2Ω ×
2Ω system reduces to a 2 × α = 1 , . . . , Ω:det (cid:0) J α − λ α I (cid:1) = det (cid:18) f u + D u Λ ( α ) − λ α f v g u g v + D v Λ ( α ) − λ α (cid:19) = 0where J α is the 2 × J α ≡ J + DΛ ( α ) with D = diag ( D u , D v ), Λ ( α ) = diag (Λ ( α ) , Λ ( α ) ).The steady state is unstable to small heterogeneous per-turbations, if λ α has a positive real part over a finiterange of modes. The dispersion relation (the largest realpart of λ α , for any given α ) can be readily computedfrom: λ α = 12 [(tr J α ) Re + γ ] + 12 [(tr J α ) Im + δ ] ι (2)where γ = s A + √ A + B δ = sgn( B ) s − A + √ A + B A = [(tr J α ) Re ] − [(tr J α ) Im ] − [(det J α ) Re ] B = 2 (tr J α ) Re (tr J α ) Im − [(det J α ) Im ] . Here sgn( · ) 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 thedefinition of a square root of a complex number. Take, z = a + bι where ι = √− √ z = ± r a + | z | b ) r − a + | z | ι ! . The instability sets in when | (tr J α ) Re | ≤ γ , a condi-tion that translates into the following inequality [24]: S (Λ ( α ) Re )[Λ ( α ) Im ] ≤ − S (Λ ( α ) Re ) , (3)where (Λ Re , Λ Im ) span the complex plane where theLaplacian eigenvalues reside. In the above relation S , S are polynomials which take the following explicit form: S ( x ) = C x + C x + C x + C x + C S ( x ) = C x + C x + C . The constants here are given by: C = σ (1 + σ ) C = (1 + σ ) ( σJ + J ) + 2tr J σ (1 + σ ) C = det J (1 + σ ) + (tr J ) σ + 2tr J (1 + σ ) ( σJ + J ) C = 2tr J (1 + σ ) det J + (tr J ) ( σJ + J ) C = det J (tr J ) and C = σ (1 − σ ) C = ( σJ + J ) (1 − σ ) C = J J (1 − σ ) . The system displays a generalized Turing instability, de-termined by the asymmetric nature of the imposed cou-pling if, after the homothetic transformation, the eigen-values D u (Λ ( α ) Re , Λ ( α ) Im ) of the Laplacian operator L enterthe region of the complex plane (Λ Re , Λ Im ) that is de-limited by inequality (3). Spectral properties of circulant matrices An n × n matrix C is circulant if takes the form C = c c n − . . . c c c c c n − . . . c ... c c . . . ... c n − . . . . . . c n − c n − c n − . . . c c . The circulant matrix C is fully specified by its firstcolumn, c = ( c , c , . . . c n − ). The other columns of C are generated as cyclic permutations of the vector c with offset equal to the column index. The normal-ized eigenvectors of a circulant matrix are given by φ i = 1 √ n (cid:0) , ω i , ω i , . . . , ω n − i (cid:1) where ω i = exp(2 πιi/n )is the n − th root of the unity and the eigenvaluesare λ i = c + c n − ω i + c n − ω i + · · · + c ω n − i where i = 0 , , . . . , n − upplementary Material A universal route to pattern formation
Malbor Asllani , Timoteo Carletti , Duccio Fanelli , Philip K. Maini MACSI, Department of Mathematics and Statistics,University of Limerick, Limerick V94 T9PX, Ireland Department of Mathematics and naXys, Namur Institute for Complex Systems,University of Namur, rempart de la Vierge 8, B 5000 Namur, Belgium Dipartimento di Fisica e Astronomia, Universit`a di Firenze,INFN and CSDC, Via Sansone 1, 50019 Sesto Fiorentino, Firenze, Italy and Mathematical Institute, University of Oxford, Woodstock Rd, OX2 6GG Oxford, UK
On the shape of the instability domain and theonset of the directed instability
The aims of this section are twofold. We will first pro-vide a detailed characterization of the region D of thecomplex plane (Λ Re , Λ Im ) where the instability eventu-ally takes place. Then, we will show that a non-zeroinclination θ ∗ exists for any value of σ strictly differentfrom one. More specifically, we shall prove that the linetangent to the instability domain displays a slope whichscales as 1 /(cid:15) , for σ = 1 ± (cid:15) and 0 < (cid:15) (cid:28)
1. This observa-tion yields an important consequence. It is in fact alwayspossible to select a minimal lattice size, Ω, so as to have θ ,the inclination of the line that connects the second eigen-value to the origin, smaller than θ ∗ >
0. This implies,in turn, that the system can be always made unstableby appropriately setting the strength of D u (and conse-quently D v ), while keeping σ fixed to a nominal valuearbitrarily close to one. The instability cannot be seededat σ = 1 ( θ ∗ , being in this case zero, or, equivalentlythe slope of the tangent to D , infinite), a zero measurelimiting condition.By using a symbolic manipulator, it is possible to showthat the fourth degree polynomial S ( x ) admits a doubleroot at x = − tr J / ( D v + D u ), which is a positive definednumber because of the necessary assumption tr J < S ( x ) = ( x + tr J ) Q ( x )where Q ( x ) = ( D v + D u )(det J +( D u J + D v J ) x + D u D v x ) , is a second order polynomial whose roots are straightfor-wardly determined as x = − ( D u J + D v J ) − p ( D u J + D v J ) − D u D v det J D u D v x = − ( D u J + D v J ) + p ( D u J + D v J ) − D u D v det J D u D v . (1)If ( D u J + D v J ) − D u D v det J >
0, then x < x ; otherwise Q ( x ) > x . In fact, det J > S ( x ) is straightforward as itis a second degree polynomial. Its roots can be written: y = − J D v and y = − J D u . (2)Without lost of generality, we can assume J > J < y > y <
0. Observe that when S van-ishes, the domain D hits the x -axis with a vertical slope(provided S is non-zero at the same time, which is gener-ically true). On the other hand, if S vanishes (but not S , which is again true in general) then the domain dis-plays a vertical asymptote (the vertical dashed line inFig. 2 in the Main Text).Collecting together the above information, we can iso-late the following possible scenarios: A: ( D u J + D v J ) − D u D v det J > D u J + D v J < x ≤ S ( x ) is always positivewhile S ( x ) is positive for x < y and negative oth-erwise. The domain D therefore displays a verticalasymptote at x = y . B: ( D u J + D v J ) − D u D v det J > D u J + D v J >
0. Here, we can identify three sub-cases,depending on the relative positions of the roots x , x and y . B.1: x < x < y (see Fig. 2). In the negativehalf-plane the domain D is composed of thefollowing parts: two unbounded regions nevertouching the x -axis and a bounded domainwith a finite intersection with the x -axis; thelatter can be responsible for the emergence ofTuring patterns on symmetric supports. S ( x )changes sign twice while S ( x ) is positive for x < y and negative otherwise. B.2: x < y < x (see Fig. 3). In the negativehalf-plane D is composed of a single domain a r X i v : . [ n li n . PS ] J un intersecting the x -axis. Hence Turing patternson symmetric supports can develop. In the in-terval ( y , x ) both S ( x ) and S ( x ) are neg-ative, hence the inequality S ( x ) y ≤ − S ( x )is satisfied for all y . B.3: y < x < x (see Fig. 4). In the negativehalf-plane D is composed again of a singledomain intersecting the x -axis. Hence Tur-ing patterns on symmetric support can possi-bly occur. In the interval ( x , x ) both S ( x )and S ( x ) are negative, hence the inequality S ( x ) y ≤ − S ( x ) is satisfied for all y . C: ( D u J + D v J ) − D u D v det J < S ( x ) is always positive while S ( x ) is positive for x < y and negative otherwise in the negative half-plane x ≤
0. The domain D displays a verticalasymptote at x = y .In the cases where Turing patterns cannot develop ona symmetric network, namely A and C (we will also con-sider case B.1 because of its peculiar shape), we are inter-ested in determining the straight line that passes throughthe origin and is tangent to the domain D . Because thelatter domain presents a clear symmetry y
7→ − y wewill limit ourselves to considering the tangency condi-tion for positive values of y . Further we shall denote by y = f ( x ) = p − S ( x ) /S ( x ) the boundary of D .We select a generic point on the boundary of D , that is,a point with coordinates ( x , f ( x )) (see Fig. 6). Then,we consider the equation of the line tangent to the do-main in x , i.e. y = f ( x )( x − x ) + f ( x ). Further, welet x vary and among all the straight lines, we select theone passing through the origin, i.e. setting x = y = 0in the previous equation. The coordinate of such a pointmust hence satisfy: 0 = − x f ( x ) + f ( x ).Using the definition of f ( x ), x must satisfy2 S ( x ) S ( x ) − x ( S ( x ) S ( x ) − S ( x ) S ( x ) = 0 ;(3)this appears to be a 6-th order polynomial. It can beshown that the coefficient of x is identically equal tozero because of the definition of C ij . We therefore havea 5-th order polynomial, that we label P ( x ). By using asymbolic manipulator, we can show that − tr J / ( D v + D u )is a root of P ( x ) [2]. We are hence finally left with P ( x ) = ( x + tr J / ( D v + D u )) P ( x ), where P ( x ) isa polynomial of 4-th degree, whose roots can be analyti-cally determined by using well known formulas.To proceed in the analysis, we write P ( x ) = e + dx + cx + bx + ax . Then the four roots read r = − b a − S + 12 r − S − p + qS (4) r = − b a − S − r − S − p + qS (5) r = − b a + S + 12 r − S − p − qS (6) r = − b a + S − r − S − p − qS , (7) where p = 8 ac − b a (8) q = b − abc + 8 a d a (9)∆ = c − bd + 12 ae (10)∆ = 2 c − bcd + 27 b e + 27 ad − ace (11) Q = " ∆ + p ∆ − / (12) S = 12 s − p + 13 a (cid:18) Q + ∆ Q (cid:19) . (13)Let us call r any negative (real) solution. Then, weneed to evaluate f ( r ), the slope of the line passingthrough the origin and tangent to the domain D . f ( r )will depend on the set of parameters involved. As antic-ipated in the beginning of this section, we aim at char-acterizing the dependence of f ( r ) on (cid:15) , when setting σ = 1 ± (cid:15) and for (cid:15) (cid:28)
1. To this end we observe that thecoefficients C i can be rewritten as: C = D u (1 + (cid:15) ) (cid:15) = ˆ C (cid:15) (14) C = D u ((1 + (cid:15) ) J + J ) (cid:15) = ˆ C (cid:15) (15) C = J J D u (cid:15) = ˆ C (cid:15) , (16)where ˆ C i are consistently defined by the above equali-ties. As a consequence, we can factor out a term (cid:15) from S ( x ). This implies that f ( x ) can be rewritten has f ( x ) = s − S ( x ) S ( x ) = 1 (cid:15) s − S ( x )ˆ S ( x ) = 1 (cid:15) ˆ f ( x ) , where ˆ S ( x ) follows the definition of the coefficients ˆ C i .We can hence apply the above analysis to the scaledfunction ˆ f ( x ), which is regular in the limit (cid:15) →
0. Inparticular, we can now select the negative root(s) ofEq. (3) (where S is replaced with ˆ S and where (cid:15) isset to zero), say ˆ r . We can then compute ˆ f (ˆ r ) andconclude that f ( r ) = ˆ f (ˆ r ) /(cid:15) . That is, the optimalslope gets steeper as (cid:15) get closer to 0. Moreover, for anygiven (cid:15) the slope is finite, as anticipated above. S ( x )
The domain D : ( D u J + D v J ) − D u D v det J > and D u J + D v J <
0. Left panel: the signs of the functions S ( x ) and S ( x ) for x ≤
0. Right panel: the corresponding domain D is depicted. Note the vertical asymptote and the fact thatthe domain has a limited extension on the horizontal axis. The domain never intercepts the x -axis. Hence, Turing patternscannot develop on a symmetric support, whose Laplacian spectrum is real. S ( x )
The domain D : ( D u J + D v J ) − D u D v det J > and D u J + D v J > , x < x < y . Left panel: thesigns of the functions S ( x ) and S ( x ) are characterized for x ≤
0. Right panel: the corresponding shape of the D domain isplotted. Note the vertical asymptote and the fact that the domain has a limited extension on the horizontal axis. The domainintercepts the x -axis over a finite interval. Turing patterns can hence develop on symmetric supports.[1] It is straightforward to prove that J and J shouldhave different signs for the generalized instability discussedhere. If J and J are both positive then the homo-geneous fixed point is unstable, a setting that we ruledout by hypothesis. Conversely, when J and J are bothnegative, there is no region in the complex plane C thatmatches the condition set by inequality (3) in the MainText. Hence, the two coupled species should be necessar-ily linked via an activator-inhibitor scheme, for the gen-eralized instability to eventually develop (i.e. J and J should display opposite signs), as in the original Turingframework. Recall however that the route to instabilityhere discussed enables one to virtually remove all con-straints on the diffusion constants, as imposed under Tur-ing’s approach.[2] In fact − tr J / ( D v + D u ) is a double root of S , hence asimple root of S . Inserting x = − tr J / ( D v + D u ) intoEq. (3) yields the identity 0 = 0. S ( x )
The domain D : ( D u J + D v J ) − D u D v det J > and D u J + D v J > , x < y < x . Left panel: the signsof the functions S ( x ) and S ( x ) for x ≤
0. Right panel: the corresponding shape of the D domain is displayed. The domainintercepts the x -axis over an extended zone. Turing patterns can hence develop on symmetric supports. S ( x )
The domain D : ( D u J + D v J ) − D u D v det J > and D u J + D v J > , y < x < x . Left panel: the signs ofthe functions S ( x ) and S ( x ) are shown for x ≤
0. Right panel: the corresponding profile of the D domain is depicted. Notethe vertical asymptote and the fact that the domain has a limited extension on the horizontal axis. The domain intercepts the x -axis over an extended zone, and Turing patterns are consequently possible on symmetric supports. S ( x )
The domain D : ( D u J + D v J ) − D u D v det J < S ( x ) and S ( x ) areshown for x ≤
0. Right panel: the corresponding shape of the D domain is given. Note the vertical asymptote and the factthat the domain displays a limited horizontal extension. The domain does not intercept the x -axis. Turing patterns cannot setin when the system is made to evolve on a symmetric support. D