Generalized patterns from local and non local reactions

TL;DR

This paper introduces a novel pattern formation mechanism in networked systems with local and non-local reactions, analyzing stability and symmetry breaking without relying on classical Turing conditions.

Contribution

It provides analytical insights into how homogeneous states become unstable through a new reactive Laplacian spectrum, enabling alternative pattern formation models.

Findings

Homogeneous states can become unstable under external disturbances.

A new reactive Laplacian spectrum influences stability.

Patterns can emerge without classical Turing conditions.

Abstract

A class of systems is considered, where immobile species associated to distinct patches, the nodes of a network, interact both locally and at a long-range, as specified by an (interaction) adjacency matrix. Non local interactions are treated in a mean-field setting which enables the system to reach a homogeneous consensus state, either constant or time dependent. We provide analytical evidence that such homogeneous solution can turn unstable under externally imposed disturbances, following a symmetry breaking mechanism which anticipates the subsequent outbreak of the patterns. The onset of the instability can be traced back, via a linear stability analysis, to a dispersion relation that is shaped by the spectrum of an unconventional reactive Laplacian. The proposed mechanism prescinds from the classical Local Activation and Lateral Inhibition scheme, which sits at the core of the Turing recipe for diffusion driven instabilities. Examples of systems displaying a fixed-point or a limit cycle, in their uncoupled versions, are discussed. Taken together, our results pave the way for alternative mechanisms of pattern formation, opening new possibilities for modeling ecological, chemical and physical interacting systems.

Equations

$$\{\begin{array}{rl}\dot{x}=& f(x,y)\\ \dot{y}=& g(x,y).\end{array}$$

$$\{\begin{array}{rl}{\dot{x}}_i=& f_0(x_i,y_i)+\frac{1}{k_i^A}\sum _j{A}_{ij}{\tilde{f}}_1(x_i,x_j,y_i,y_j)\\ {\dot{y}}_i=& g_0(x_i,y_i)+\frac{1}{k_i^A}\sum _j{A}_{ij}{\tilde{g}}_1(x_i,x_j,y_i,y_j).\end{array}$$

$$\{\begin{array}{rl}{\dot{x}}_i=& f(x_i,y_i)+\sum _j{L}_{ij}^{\mathrm{R}}{\tilde{f}}_1(x_i,x_j,y_i,y_j)\\ {\dot{y}}_i=& g(x_i,y_i)+\sum _j{L}_{ij}^{\mathrm{R}}{\tilde{g}}_1(x_i,x_j,y_i,y_j)\end{array}$$

$$\left(\begin{array}[]{ c }\dot{\bm{u}}\\ \dot{\bm{v}}\end{array}\right)=\bm{J}\otimes\mathbb{1}\left(\begin{array}[]{ c }\bm{u}\\ \bm{v}\end{array}\right)+\bm{J}_{R}\otimes\bm{{\operatorname{\mathcal{L}}}}^{\rm R}\left(\begin{array}[]{ c }\bm{u}\\ \bm{v}\end{array}\right)$$

$$\{\begin{array}{rl}{\dot{x}}_i=& -d{x}_i+\beta \frac{c_1}{k_i}{y}_i\sum _j{A}_{ij}{x}_j+(1-\beta )\frac{c_1}{k_i}{x}_i\sum _j{A}_{ij}{y}_j\\ {\dot{y}}_i=& r{y}_i-\frac{s}{k_i}{y}_i\sum _j{A}_{ij}{y}_j-\frac{c_2}{k_i}{y}_i\sum _j{A}_{ij}{x}_j\end{array}$$

$${\dot{x}}_i$$

$${\dot{y}}_i$$

$$\{\begin{array}{rl}{\dot{u}}_i=& {\partial }_x{f}^{\ast }{u}_i+{\partial }_y{f}^{\ast }{v}_i+\sum _j{L}_{ij}^{\mathrm{R}}[{\partial }_2{f}_R^{\ast }{u}_j+{\partial }_4{f}_R^{\ast }{v}_j]\\ {\dot{v}}_i=& {\partial }_x{g}^{\ast }{u}_i+{\partial }_y{g}^{\ast }{v}_i+\sum _j{L}_{ij}^{\mathrm{R}}[{\partial }_2{g}_R^{\ast }{u}_j+{\partial }_4{g}_R^{\ast }{v}_j]\end{array}$$

$$\left(\begin{array}[]{ c }\dot{u_{i}}\\ \dot{v_{i}}\end{array}\right)=\bm{J}\left(\begin{array}[]{ c }u_{i}\\ v_{i}\end{array}\right)+\sum_{j}\operatorname{\mathcal{L}}_{ij}^{\rm R}\bm{J}_{R}\left(\begin{array}[]{ c }u_{j}\\ v_{j}\end{array}\right)$$

$$\bm{J}_{R}=\left(\begin{array}[]{ c c}\partial_{2}f_{R}^{*}&\partial_{4}f_{R}^{*}\\ \partial_{2}g_{R}^{*}&\partial_{4}g_{R}^{*}\end{array}\right).$$

$$\left(\begin{array}[]{ c }\dot{\bm{u}}\\ \dot{\bm{v}}\end{array}\right)=\bm{J}\otimes\mathbb{1}\left(\begin{array}[]{ c }\bm{u}\\ \bm{v}\end{array}\right)+\bm{J}_{R}\otimes\bm{{\operatorname{\mathcal{L}}}}^{\rm R}\left(\begin{array}[]{ c }\bm{u}\\ \bm{v}\end{array}\right)$$

$$u_i=\sum _{\alpha }{c}_{\alpha }{\varphi }_i^{(\alpha )},v_i=\sum _{\alpha }{d}_{\alpha }{\varphi }_i^{(\alpha )}.$$

$$\sum_{\alpha}\left(\begin{array}[]{ c }\dot{c_{\alpha}}\\ \dot{d_{\alpha}}\end{array}\right)\bm{\phi}^{(\alpha)}=\bm{J}\otimes\mathbb{1}\sum_{\alpha}\left(\begin{array}[]{ c }c_{\alpha}\\ d_{\alpha}\end{array}\right)\bm{\phi}^{(\alpha)}+\bm{J}_{R}\otimes\bm{{\operatorname{\mathcal{L}}}}^{\rm R}\sum_{\alpha}\left(\begin{array}[]{ c }c_{\alpha}\\ d_{\alpha}\end{array}\right)\bm{\phi}^{(\alpha)},$$

$$\sum_{\alpha}\left(\begin{array}[]{ c }\dot{c_{\alpha}}\\ \dot{d_{\alpha}}\end{array}\right)\bm{\phi}^{(\alpha)}=[\bm{J}+\Lambda^{(\alpha)}\bm{J}_{R}]\otimes\mathbb{1}\sum_{\alpha}\left(\begin{array}[]{ c }c_{\alpha}\\ d_{\alpha}\end{array}\right)\bm{\phi}^{(\alpha)}.$$

$$\left(\begin{array}[]{ c }\dot{c_{\alpha}}\\ \dot{d_{\alpha}}\end{array}\right)=[\bm{J}+\Lambda^{(\alpha)}\bm{J}_{R}]\left(\begin{array}[]{ c }c_{\alpha}\\ d_{\alpha}\end{array}\right).$$

$$\lambda_{\pm}(\alpha)=\frac{1}{2}\biggl{[}\operatorname{Tr}(\bm{J}_{\alpha})\pm\sqrt{\operatorname{Tr}(\bm{J}_{\alpha})^{2}-4\det(\bm{J}_{\alpha})}\biggr{]}.$$

$$\Lambda_{\pm}=\frac{1}{2}\biggl{[}-\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})\pm\sqrt{[\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})]^{2}-4\det(\bm{J}^{-1}\bm{J}_{R})}\biggr{]}.$$

$$\left(\begin{array}[]{ c }u_{i}(t)\\ v_{i}(t)\end{array}\right)=\sum_{\alpha}\left(\begin{array}[]{ c }c_{\alpha}(0)\\ d_{\alpha}(0)\end{array}\right)e^{\lambda_{\pm}(\alpha)t}\phi_{i}^{(\alpha)}.$$

Full text

I Introduction

Reaction-diffusion systems are universally invoked as a reference modeling tool, owing to their inherent ability to sustain the spontaneous generation of – spatially extended – patterned motifs. The theory of pattern formation was laid down by A. Turing Turing (1952) and later on fertilized in a cross-disciplinary perspective to eventually become a pillar to explain self-organization in nature Strogatz (2001); Pecora et al. (1997); Pismen (2006). In the original Turing setting, two chemicals, here termed species, can relocate in space following Fickean diffusion. In general, a form of local feedback is needed: this activates the short range production of a given species, which should be, at the same time, inhibited at long ranges. In standard reaction-diffusion systems, this combination is accomplished via autocatalytic reaction loops and having the inhibitors diffusing faster than activators. The Local Activation and Lateral Inhibition (LALI) paradigm provides hence the reference frame for the onset of self-organized patterns of the Turing type.

Both in the continuous and discrete (lattice or network based  Othmer and Scriven (1974); Nakao and Mikhailov (2010); Asllani et al. (2014)) versions, diffusion is a key ingredient for the Turing paradigm to apply. Assume in fact that a stable fixed-point (or, alternatively a time dependent equilibrium, e.g. a limit cycle) exists for the scrutinized reaction scheme in its a-spatial local version. Then, the spatial counterpart of the inspected model admits a homogeneous equilibrium. This follows trivially by recalling that by definition the Laplacian, the operator that implements the diffusive transport, returns zero when applied to a uniform density background. Tiny erratic disturbances, that happen to perturb the uniform solution, may therefore activate the diffusion term, and consequently feed on the reaction part, in this way promoting a self-sustained instability which sits at the core of the Turing mechanism.

Natural systems display however an innate drive towards self-organization, which certainly transcends the realms of validity of Turing’s ideas. In many cases of interest, the quantities to be monitored are anchored to the nodes of a virtual graph. The information to be processed on site is carried across the edges of the graph. Nevertheless, the spreading of information springing from nearby nodes is not necessarily bound to obey linear diffusion, as in the spirit of the original Turing formulation. At variance, it can in principle comply with a large plethora of distinct dynamical modalities. In ecology, the number of nodes reflects biodiversity: inter-species interactions are usually epitomized by a quadratic response function, which accounts for competitive predator-prey or symbiotic dependences May (1972); Thébault and Fontaine (2010). In genetic networks, inputs between neighboring genes are shaped by a sigmoidal profile, a non-linear step function often invoked to mimic threshold activation processes Becskei and Serrano (2000). At a radically different scale, the dynamics of opinions in sizable social systems are customarily traced back to pairwise exchanges, which can echo physical encounters or be filtered by device-assisted long-ranged interactions Castellano et al. (2009). Similarly, global pandemic events mirror complex migration patterns, which cannot be explained as resorting to the simplistic archetype of diffusion Pastor-Satorras and Vespignani (2001).

A first insight into the study of pattern formation for a system made of immobile species was offered in  Bullara and De Decker (2015). In this latter paper the colored motifs on the skin of a zebrafish were reproduced by postulating an experimentally justified long ranged regulatory mechanism, in absence of cell motility. The differential growth mechanism identified by Bullara and De Decker Bullara and De Decker (2015) yields a LALI paradigm alternative to that stemming from conventional reaction diffusion schemes.

Motivated by these observations, we here take one leap forward by proposing a generalized theory of pattern formation which is not bound to satisfy the LALI conditions. Interactions are assumed to be mediated by local (on site) and non local (distant, possibly long-ranged) pairings among species. More specifically, we will assume that non local interactions, as specified by a off-diagonal elements of a binary or weighted adjacency matrix, can be scaled by the node degree, following a mean-field working ansatz. The elemental units of the system being inspected (species) are permanently linked to the nodes they are bound to, and therefore classified as immobile. Local interactions will be also referred to as reactions, to establish an ideal bridge with standard reaction-diffusion systems.

II The scheme of interaction

Consider a reactive system made of two, mutually interacting species and label with $x$ and $y$ their respective concentrations. Assume:

[TABLE]

where $f(\cdot,\cdot)$ and $g(\cdot,\cdot)$ are generic non-linear functions. Suppose that the system admits an equilibrium solution, $(x^{*},y^{*})$, which can be either a fixed-point or a limit cycle. Linearizing equations (1) around ${\bm{w}^{*}}=(x^{*},y^{*})$ returns a (linear) first order system for the evolution of the imposed perturbation $\delta{\bm{w}}=(\delta x,\delta y)^{T}$, where $T$ stands for the transpose operation. In formulae, one gets $\delta\dot{\bm{w}}=\bm{J}(\bm{w}^{*})\delta{\bm{w}}$, where \bm{J}(\bm{w}^{*})=\left(\begin{array}[]{ c c}\partial_{x}f^{*}&\partial_{y}f^{*}\\ \partial_{x}g^{*}&\partial_{y}g^{*}\end{array}\right) is the Jacobian matrix of the system evaluated at equilibrium. When the equilibrium is a fixed-point, the stability of ${\bm{w}^{*}}$ is set by the spectrum of $\bm{J}$: if the largest real part of the eigenvalues of $\bm{J}$ is negative the imposed perturbation fades away exponentially and the fixed-point is deemed stable. This amounts to requiring $\operatorname{Tr}(\bm{J})<0$ and $\det(\bm{J})>0$, a pre-requisite condition for our analysis to apply. When the equilibrium solution is a limit cycle, the Jacobian exhibits a periodic dependence on time and the stability can be assessed by evaluating the Floquet exponents Grimshaw (2017); Challenger et al. (2015); Lucas et al. (2018). In the following, we develop our reasoning by assuming a constant Jacobian. Our analysis and conclusions, however, readily extend to a limit cycle setting, provided one replaces $\bm{J}$ with the time-independent Floquet matrix.

Starting from a local formulation of the reactive model, we move forward and account for non local interactions. We begin by replicating the examined system on $N$ different patches, defining lumps of embedding space or specific criteria to effectively group the variables involved. The collection of $N$ isolated spots where the local interaction scheme is replicated are the nodes of a network, identified by the index $i=1,...,N$. The edges of the network provide a descriptive representation of the non-local interaction scheme: $A_{ij}=1$, if nodes $i$ and $j$ share a link, i.e. if they are bound to mutually interact (the formalism extends readily to account for weighted interactions). For the sake of simplicity we limit the analysis to undirected networks, $A_{ij}=A_{ji}$. We note that, by definition, we have $A_{ii}=1$. The connectivity of node $i$ is hence $k_{i}^{A}=\sum_{j}A_{ij}$. Rewrite the non-linear function $f$ as $f(x,y)=f_{0}(x,y)+f_{1}(x,y)$, where $f_{0}$ speaks for the interaction terms that are genuinely local, i.e. those not modulated by distant couplings. Conversely, $f_{1}(x,y)$ will be replaced by a reaction function, which encodes for the remote interaction. In the spirit of a nearest neighbors mean-field approximation, we shall substitute $f_{1}(x,y)$ with its averaged counterpart in the heterogeneous space of the interaction network. Similarly, $g$ will be split as $g(x,y)=g_{0}(x,y)+g_{1}(x,y)$. Then we set:

[TABLE]

where $\tilde{f}_{1}$ (resp. $\tilde{g}_{1}$) reduces to $f_{1}$ (resp. $g_{1}$) when the system is in a homogeneous state. Mathematically, we require $\tilde{f}_{1}(x_{i},x_{j},y_{i},y_{j})\delta_{ij}=f_{1}(x_{i},y_{i})$ and $\tilde{g}_{1}(x_{i},x_{j},y_{i},y_{j})\delta_{ij}=g_{1}(x_{i},y_{i})$. As an illustrative example, consider a local quadratic interaction of the type $f_{1}(x_{i})=x_{i}^{2}$. This latter can be turned into $\tilde{f}_{1}(x_{i},x_{j})=x_{i}x_{j}$, a choice that matches the above constraints. This would yield $(x_{i}/k_{i}^{A})\sum_{j}A_{ij}x_{j}=x_{i}\langle x\rangle_{i}$ in the first of equations (2), where the average $\langle\cdot\rangle_{i}$ runs over the nodes adjacent to $i$. In practical terms, as anticipated above, we replace the local reaction term $f_{1}(x_{i},y_{i})$ to include interactions at distance, where the average runs over the $k_{i}^{A}$ terms defining the neighborhood of $i$. The remote interaction is assumed to depend only on the species that populate the selected pair of nodes.

By adding and subtracting $f_{1}(x_{i},y_{i})$ (resp. $g_{1}(x_{i},y_{i})$ ) in the first (resp. second) of equations (2) and making use of the definition of $f(\cdot,\cdot)$ and $g(\cdot,\cdot)$ yields:

[TABLE]

where $\operatorname{\mathcal{L}}_{ij}^{\rm R}=\frac{A_{ij}}{k_{i}^{A}}-\delta_{ij}$. Observe that the above equations are formally identical to system (1) when $\bm{A}\equiv\mathbb{1}_{N}$, where $\mathbb{1}_{N}$ stands for the $N\times N$ identity matrix. In this case, in fact, the $N$ nodes are formally decoupled and system (3) reduces to (1), for each node of the ensemble. System (3) admits an obvious homogeneous solution. This latter is found by assigning $(x_{i},y_{i})=(x^{*},y^{*})$ $\forall i$. Can one make the homogeneous equilibrium unstable to tiny non homogeneous perturbations by acting on the network of interactions, as encoded in the operator $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$? Answering this question amounts to develop a general theory for pattern formation for interacting system, in absence of species diffusion. To expand along this line one needs to characterize the operator $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$. As shown in the Appendix, it is easy to prove that $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$ displays a non-positive spectrum and that has an eigenvector $(1,\dots,1)^{T}$ associated to the [math] eigenvalue. This justifies referring to $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$ as to a reactive Laplacian. In the literature it is also known as the consensus Laplacian Ghosh et al. (2014); Lambiotte et al. (2011); Olfati-Saber et al. (2007); Krause (2009). It is important to notice that this is different from the diffusive Laplacian ($L^{diff}_{ij}=A_{ij}-k_{j}^{A}\delta_{ij}$) and its normalized version ($L^{RW}_{ij}=\frac{A_{ij}}{k_{j}^{A}}-\delta_{ij}$) which describes random walk. It is not possible to reduce the reactive Laplacian to the Laplacian of diffusion, unless the complex network reduces to a simple regular lattice.

III Linear stability analysis

The stability of the homogeneous equilibrium $(x_{i},y_{i})=(x^{*},y^{*})$, $\forall i$ can be analytically probed Pecora and Carroll (1998). Introduce a small inhomogeneous perturbation, $x_{i}=x^{*}+u_{i}$, $y_{i}=y^{*}+v_{i}$, and expand the governing equations up to the leading (linear) order to yield:

[TABLE]

where $\bm{u}=\left(u_{1},u_{2},...,u_{N}\right)$, $\bm{v}=\left(v_{1},v_{2},...,v_{N}\right)$, \bm{J}_{R}=\left(\begin{array}[]{ c c}\partial_{2}\tilde{f}_{1}^{*}&\partial_{4}\tilde{f}_{1}^{*}\\ \partial_{2}\tilde{g}_{1}^{*}&\partial_{4}\tilde{g}_{1}^{*}\end{array}\right). The symbols $\partial_{2}$ and $\partial_{4}$ indicate that the derivatives are computed with respect to the second and fourth arguments of the function to which they are applied. The symbol $\otimes$ refers to the Kronecker product. More details on the derivation of eq. (4) are given in the Appendix. To solve the above linear system, we introduce the eigenvectors ($\phi_{j}^{(\alpha)}$) and eigenvalues ($\Lambda^{(\alpha)}$) of the reactive Laplacian $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$, i.e. $\sum_{j}\operatorname{\mathcal{L}}_{ij}^{\rm R}\phi_{j}^{(\alpha)}=\Lambda^{(\alpha)}\phi_{i}^{(\alpha)}$ with $\alpha=1,...,N$. Expanding the perturbations on the basis of $\phi_{j}^{(\alpha)}$, $\alpha=1,...,N$, we obtain (see Appendix) the $2\times 2$ modified Jacobian $\bm{J}_{\alpha}=\bm{J}+\Lambda^{(\alpha)}\bm{J}_{R}$, which determines the stability of the system. For each choice of $\Lambda^{(\alpha)}$, one needs to compute $\lambda(\alpha)$, the largest real part of the eigenvalues of $\bm{J}_{\alpha}$. The discrete collection of $\lambda(\alpha)$ vs. $-\Lambda^{(\alpha)}$ defines the so-called dispersion relation. The homogeneous equilibrium is stable if $\lambda(\alpha)<0$, $\forall\alpha$. Conversely, the perturbation grows if at least one $\lambda(\alpha)$ turns out to be positive. In this latter case patterns can emerge, thus breaking the symmetry of the unperturbed initial condition. The condition of instability can be expressed in a compact form (provided $\bm{J}^{-1}$ exists), as we prove in the Appendix. Specifically, we find that the instability sets in if one of the following conditions holds: (i) if $\operatorname{Tr}(\bm{J}_{R})<0$ and provided $\Lambda^{(\alpha)}<-\frac{\operatorname{Tr}(\bm{J})}{\operatorname{Tr}{\bm{J}_{R}}}$, for at least one choice of $\alpha$; (ii) if $\det(\bm{J}_{R})<0$ and provided $\Lambda^{(\alpha)}<\Lambda_{-}$, for at least one choice of $\alpha$; (iii) if $\det(\bm{J}_{R})>0$, $\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})>0$ and provided $\Lambda_{-}<\Lambda^{(\alpha)}<\Lambda_{+}$, for at least one choice of $\alpha$. Here, \Lambda_{\pm}=\frac{1}{2}\biggl{[}-\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})\pm\sqrt{[\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})]^{2}-4\det(\bm{J}^{-1}\bm{J}_{R})}\biggr{]}. The above conclusions hold for a two species model that displays a constant homogeneous equilibrium. As such, it generalizes the renowned Turing instability conditions for reaction-diffusion models, to systems made of immobile species subject to local and non local (remote) interactions. The analysis extends to the setting where the homogeneous state is a synchronous collection of oscillators performing in unison. As we shall also remark in the following, the novel route to pattern formation overcomes the LALI constraint.

Before turning to discussing the applications, we remark that generalized transport schemes have been considered in the past which enable one to relax the LALI request. This includes for instance accounting for cross-diffusion Vanag and Epstein (2009); Madzvamuse et al. (2015). It should be emphasized however that within the current scheme, species are immobile, hence permanently bound to their reference nodes. It is in fact the combination of local and non-local interaction rules that materializes in the discrete non-local operator, different from the Laplacian of diffusion, which rules the dynamics of the system in its linear approximation.

IV Applications and discussions

As a first example, we consider the reaction scheme known as the Volterra model McKane and Newman (2005), $f(x,y)=c_{1}xy-dx$ and $g(x,y)=ry-sy^{2}-c_{2}xy$. The variable $x$ identifies the concentration of predators, while $y$ stands for the prey. The parameters are positive definite. The Volterra equations admit a fixed-point, $x^{*}=\frac{c_{1}r-sd}{c_{1}c_{2}}$, $y^{*}=\frac{d}{c_{1}}$, which is meaningful and stable provided $c_{1}r-sd>0$. We now turn to considering $N$ replicas of the model, by making the variables to depend on the node index $i$, associated to different ecological niches. Following the above lines of reasoning, we assume that species can sense the remote interaction with other communities populating the neighboring nodes. For instance, the competition of prey for food and resources can be easily extended so as to account for a larger habitat which embraces adjacent patches. At the same time, predators can benefit from a coordinated action to hunt in team. Assume that the linear terms in $f(x,y)$ and $g(x,y)$ stay local (death term for predators and birth term for prey), i.e. $f_{0}(x,y)=-dx$, $g_{0}(x,y)=ry$. We generalize the terms $f_{1}(x,y)=c_{1}xy$ and $g_{1}(x,y)=-sy^{2}-c_{2}xy$ to the case of remote interactions as:

[TABLE]

where $\beta\in(0,1)$ modulates the relative strength of the non local terms that drive the evolution of the predators. By simple manipulations, we cast the previous equations so as to bring into evidence the dependence on the reactive Laplacian $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$. In formulae, one gets

[TABLE]

The system can then be studied by resorting to the strategy developed above. The interaction parameters can be set so as to make the dispersion relation positive, as displayed in Fig. 1(c). Red symbols refer to the generalized Volterra model organized on a close lattice made of $N=30$ nodes (panel (a) of Fig. 1), while the blue symbols are obtained when the sites are connected by an Erdős-Rényi (ER) network with $\langle k\rangle=3.3$ (panel (b)). The colors of the nodes stand for the asymptotic stationary state of the system across the networks, showing the formation of distinct patterns for predators (outer ring) and prey (inner ring), in the two cases. While the lattice yields a regular, equally spaced, pattern, a disordered distribution of populated patches is instead obtained for the ER network, thus reflecting the randomness of the connections. We remark that in both cases, the extinction in a patch of one of the species is typically associated to the disappearance of the other one in the same node, resulting in a global habitat with populated patches alternated with uninhabited ones. These latter play the role of natural barriers separating colonized niches. The emerging patterns of coexistence differ from those observed from standard reaction-diffusion systems, where predators and prey tend to cluster in distinct sites. The mutually exclusive distribution of species, as obtained in classical reaction diffusion scheme, bears the imprint of the LALI assumption. This latter does not apply for the patterns reported in Fig 1. As a matter of fact, the associated Jacobian matrix displays non positive diagonal entries. As a side result, we notice that in some cases, Fig. 1(d), the perturbation may trigger the emergence of regular oscillations, typical of a predator-prey cycle. This behavior cannot be reproduced for a deterministic Volterra model in absence of remote couplings McKane and Newman (2005).

As a second example, we consider a set of interacting Stuart-Landau oscillators Bordyugov et al. (2010); Teramae and Tanaka (2004); Selivanov et al. (2012); Panaggio and Abrams (2015); Zakharova et al. (2014). The Stuart-Landau equation takes the form $\dot{w}=\sigma w-\beta w|w|^{2}$, where $w=w_{Re}+iw_{Im}$ is a complex number, as well as the parameters $\sigma,\beta$. The above equation yields a stable limit cycle which can be cast in the explicit form $w^{*}(t)=\sqrt{\frac{\sigma_{Re}}{\beta_{Re}}}\exp(i\omega t)$, where $\omega=\sigma_{Im}-\beta_{Im}\sigma_{Re}/\beta_{Re}$. To determine the stability of the time dependent periodic solution $w^{*}$, one can introduce a perturbation of $w^{*}$ in the form $w^{*}\left(1+u\right)\exp(v)$, where $u$ and $v$ are small quantities. Expanding at the linear order of approximation yields a Jacobian matrix \bm{J}=\left(\begin{array}[]{ c c}-2\sigma_{Re}&0\\ -2\beta_{Im}\sigma_{Re}/\beta_{Re}&0\end{array}\right). Interestingly, even if $w^{*}$ is periodic in time, $\bm{J}$ is constant, owing to the specificity of the Stuart-Landau equation and to the form of the imposed perturbation. This observation allows us to proceed in the analysis without resorting to the Floquet machinery. The request of a stable limit cycle $w^{*}$ translates into $\sigma_{Re}>0$, which consequently implies $\beta_{Re}>0$ (since $\sigma_{Re}/\beta_{Re}>0$ for the solution to exist). Stuart-Landau oscillators are often coupled diffusively, either on a regular lattice or on a complex network, so as to result in the celebrated Ginzburg-Landau equation. For suitable choices of the parameters, the oscillators evolve in unison, giving rise to a fully synchronized homogenous state. This latter can however break apart, upon injection of a tiny disturbance, if the parameters are assigned in the complementary reference domain. This modulational instability, known in the literature as the Benjamin-Feir instability, follows from the subtle interplay between reaction and diffusion terms Di Patti et al. (2016); Benjamin and Feir (1967). Inspired by the analysis carried out above, here we consider a collection of non diffusive Stuart-Landau oscillators, interacting via their linear reaction term. This choice materializes in a variant of the traditional Ginzburg-Landau equation, ensuing a novel instability of the Benjamin-Feir type. Consider an ensemble made of $i=1,...,N$ oscillators at the nodes of a network, characterized by the complex amplitude $w_{i}$, with $i=1,...,N$. The family of Stuart-Landau oscillators evolves therefore as: $\dot{w}_{i}=\frac{\sigma}{k_{i}}\sum_{j}A_{ij}w_{j}-\beta w_{i}|w_{i}|^{2}$ or, equivalently, $\dot{w}_{i}=\sigma w_{i}-\beta w_{i}|w_{i}|^{2}+\sigma\sum_{j}\operatorname{\mathcal{L}}_{ij}^{\rm R}w_{j}$. A linear stability analysis which moves from the ansatz $w_{i}=w^{*}\left(1+u_{i}\right)\exp(v_{i})$ (with $u_{i},v_{i}<<1$) and implements the above calculation strategy leads to the $\bm{J}_{\alpha}=\bm{J}+\Lambda^{(\alpha)}\bm{J}_{R}$, where \bm{J}_{R}=\left(\begin{array}[]{ c c}\sigma_{Re}&-\sigma_{Im}\\ \sigma_{Im}&\sigma_{Re}\end{array}\right). As $\operatorname{Tr}(\bm{J}_{R})=2\sigma_{Re}>0$, the homogeneous solution turns unstable only if $\det(\bm{J}_{\alpha})<0$. It is easy to prove that $\det(\bm{J}_{\alpha})<0$ given $\Lambda^{(\alpha)}>2\frac{\sigma_{Re}^{2}+\sigma_{Re}\sigma_{Im}\beta_{Im}/\beta_{Re}}{\sigma_{Re}^{2}+\sigma_{Im}^{2}}$, for at least one choice of $\alpha$. Since the eigenvalues of the reactive Laplacian are by definition negative and $\sigma_{Re}>0$, the instability request materializes in the compact necessary condition $\sigma_{Re}+\sigma_{Im}\frac{\beta_{Im}}{\beta_{Re}}<0$, which constitutes the revisitation of the Benjamin-Feir instability to the present setting. Fig. 2(a) displays the region of instability in the reference plane ($\sigma_{Im},\beta_{Im}$) (when setting $\sigma_{Re}$ and $\beta_{Re}$ to unit). The dashed white line marks the transition to the region of instability. Emerging patterns, and their associated dispersion relations, $\lambda(\alpha)$ vs. $-\Lambda^{(\alpha)}$, are depicted for the Stuart-Landau oscillators distributed on, respectively, a lattice with non-local couplings and a Watts-Strogatz graph. In the Appendix, we show that long-range couplings promotes the stability of the system.

Summing up, we have here introduced a novel mechanism for the emergence of coherent patterns for systems made of interacting, although immobile, species. Remote interactions, treated in a mean-field approximation, lead to a reactive Laplacian, whose spectrum ultimately sets the conditions for the instability. Taken all together, our theory paves the way for investigating novel mechanisms of pattern formation, by providing a flexible alternative to traditional schemes that implements the LALI constraint.

.1 Characterizing the spectrum of the reactive Laplacian.

Recalling that $A_{ii}=1$, let us introduce matrix $\bm{{\mathcal{A}}}=\bm{A}-\mathbb{1}$, the adjacency matrix of the same network, after removal of the self-loops. Hence, $k_{i}^{\mathcal{A}}=\sum_{j}\mathcal{A}_{ij}=k_{i}^{A}-1$ for all $i$. Define then $L_{ij}^{A}=A_{ij}-k_{i}^{A}\delta_{ij}$ and $L_{ij}^{\mathcal{A}}=\mathcal{A}_{ij}-k_{i}^{\mathcal{A}}\delta_{ij}$. We can straightforwardly prove that $L_{ij}^{A}=L_{ij}^{\mathcal{A}}$. Indeed, $L_{ii}^{A}=A_{ii}-k_{i}^{A}=1-k_{i}^{A}=-k_{i}^{\mathcal{A}}=L_{ii}^{\mathcal{A}}$ and $\forall i\neq j$ $L_{ij}^{A}=A_{ij}=\mathcal{A}_{ij}=L_{ij}^{\mathcal{A}}$. $L_{ij}^{\mathcal{A}}$ is by definition a Laplacian matrix, hence characterised by a non-positive spectrum. Based on the above the same conclusion applies to $L_{ij}^{A}$. In particular, $(1,\dots,1)^{T}$ is the eigenvector of $\bm{L}^{A}$ associated to the [math] eigenvalue. We now turn to considering the reactive Laplacian $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$ and set to prove that it also displays a non-positive spectrum, thus providing an a posteriori justification to the name assigned to the operator. Denote by $\bm{D}$ the diagonal matrix of the connectivities of $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$, $D_{ii}=k_{i}^{A}$. Then, introduce $\bm{L}^{sym}=\bm{D}^{-1/2}\bm{L}^{A}\bm{D}^{-1/2}$. By using the definition of $\bm{L}^{A}$ we get $L^{sym}_{ij}=A_{ij}/\sqrt{k_{i}^{A}k_{j}^{A}}-\delta_{ij}$ from which it immediately follows that $\bm{L}^{sym}$ is symmetric. We can then show that $\bm{L}^{sym}$ is non-positive definite. Take any $\bm{x}\in\mathbb{R}^{N}\setminus\{0\}$, $N$ standing for the dimension of the matrices, i.e. the size of the network. Then $(\bm{x},\bm{L}^{sym}\bm{x})=(\bm{x},\bm{D}^{-1/2}\bm{L}^{A}\bm{D}^{-1/2}x)=(\bm{D}^{-1/2}\bm{x},\bm{L}^{A}\bm{D}^{-1/2}\bm{x})$ where we used the property of the scalar product $(\bm{x},\bm{M}\bm{x})=(\bm{M}^{T}\bm{x},\bm{x})$ and the fact the $\bm{D}^{-1/2}$ is diagonal. By setting $\bm{y}=\bm{D}^{-1/2}\bm{x}$, we get $(\bm{x},\bm{L}^{sym}\bm{x})=(\bm{y},\bm{L}^{A}\bm{y})\leq 0\quad\forall\bm{x}$ where the last inequality follows from the fact that $\bm{L}^{A}$ is non-positive definite. Hence $\bm{L}^{sym}$ shares the same property. We are finally in a position to show that $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$ has the characteristic spectrum of a Laplacian matrix. Observe that $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}=\bm{D}^{-1}\bm{L}^{A}=\bm{D}^{-1/2}\left(\bm{D}^{-1/2}\bm{L}^{A}\bm{D}^{-1/2}\right)\bm{D}^{1/2}=\bm{D}^{-1/2}\bm{L}^{sym}\bm{D}^{1/2}$. Hence, $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$ is similar to $\bm{L}^{sym}$ and, thus, the two operators display the same non-positive spectrum. Finally, by the very first definition, it follows that $(1,\dots,1)^{T}$ is the eigenvector associated to the [math] eigenvalue of $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$, which ends the proof.

.2 The dispersion relation and the conditions for instability

In the following we provide a detailed procedure to characterize the linear stability of the homogeneous equilibria $(x^{*},y^{*})$ for systems in the general form of eqs. (2) in the main text.

First of all, we introduce a small inhomogeneous perturbation, $x_{i}=x^{*}+u_{i}$, $y_{i}=y^{*}+v_{i}$, and expand the governing equations up to the linear order to eventually get:

[TABLE]

where the symbols $\partial_{2}$ and $\partial_{4}$ indicate that the derivatives are computed with respect to the second and fourth arguments of the function to which they are applied. In matrix form,

[TABLE]

where

[TABLE]

Eq. (7) can be written in a more compact form by making use of the Kronecker product, $\otimes$:

[TABLE]

Introduce the eigenvectors and eigenvalues of the Laplacian $\bm{{\operatorname{\mathcal{L}}}}^{\rm R}$, as $\sum_{j}\operatorname{\mathcal{L}}_{ij}^{\rm R}\phi_{j}^{(\alpha)}=\Lambda^{(\alpha)}\phi_{i}^{(\alpha)}$. Expand the perturbations on the basis of the reactive Laplacian to eventually yield:

[TABLE]

We can therefore write equation (9) as

[TABLE]

and by diagonalizing the reactive Laplacian:

[TABLE]

The above linear system is ruled by a $2\times N$ block-diagonal Jacobian matrix. Each block is the $2\times 2$ reduced Jacobian $\bm{J}_{\alpha}\equiv\bm{J}+\Lambda^{(\alpha)}\bm{J}_{R}$. By making use of the fact that the eigenvectors are linearly independent one eventually gets system

[TABLE]

Starting from this setting we can elaborate on the conditions that drive the system unstable. The eigenvalues of the reduced Jacobian $\bm{J}_{\alpha}$ are given by:

[TABLE]

It is straightforward to verify that, in absence of couplings, i.e. for $\Lambda^{(\alpha)}=0$, the real part of $\lambda_{\pm}$ is always negative, as it should be because of the stability of the uncoupled equilibrium. We now seek to find the condition for which there exists at least one eigenvalue of the reactive Laplacian, that corresponds to a positive entry of the dispersion relation, i.e. $\operatorname{Re}(\lambda_{\pm})>0$, so signaling the instability. This condition is clearly met if $\operatorname{Tr}(\bm{J}_{\alpha})>0$ or if $\det(\bm{J}_{\alpha})<0$.

The first condition can be written as $\operatorname{Tr}(\bm{J})+\Lambda^{(\alpha)}\operatorname{Tr}(\bm{J}_{R})>0$. By recalling that the Laplacian eigenvalues are negative and that $\operatorname{Tr}(\bm{J})<0$, one finds that the latter condition is matched only if $\operatorname{Tr}(\bm{J}_{R})<0$, for $\Lambda^{(\alpha)}<-\frac{\operatorname{Tr}(\bm{J})}{\operatorname{Tr}{\bm{J}_{R}}}$.

The second condition can be reformulated as $\det(\bm{J})\det(\mathbb{1}+\Lambda^{(\alpha)}\bm{J}^{-1}\bm{J}_{R})<0$, if $\bm{J}$ is invertible, which is equivalent to requiring $\det(\mathbb{1}+\Lambda^{(\alpha)}\bm{J}^{-1}\bm{J}_{R})={\Lambda^{(\alpha)}}^{2}\det(\bm{J}^{-1}\bm{J}_{R})+\Lambda^{(\alpha)}\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})+1\equiv H(\Lambda^{(\alpha)})<0$. By observing that $\det(\bm{J}^{-1}\bm{J}_{R})=\frac{\det(\bm{J}_{R})}{\det(\bm{J})}$ and that $\det(\bm{J})>0$, we can isolate two distinct scenarios: if $\det(\bm{J}_{R})<0$, the curve $H(\Lambda^{(\alpha)})$ draws a parabola like the one depicted with a green line in Fig 3. The sought condition is hence matched only if there exists an eigenvalue of the reactive Laplacian smaller than the value where the (continuous) parabola intercepts the $x$-axis. If instead $\det(\bm{J}_{R})>0$, the curve $H(\Lambda^{(\alpha)})$ is negative defined over a finite interval (blue curve in fig 3), provided $\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})>0$ and, at the same time, $\Delta\equiv[\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})]^{2}-4\det(\bm{J}^{-1}\bm{J}_{R})>0$. It is however easy to prove that the latter inequality is always true. Indeed, one can always write $\Delta=[\operatorname{Tr}(\bm{\Gamma})]^{2}-4\det(\bm{\Gamma})$ and, by denoting with $\mu_{1},\mu_{2}$ the eigenvalues of the generic matrix $\bm{\Gamma}$, we can write $\Delta=(\mu_{1}+\mu_{2})^{2}-4\mu_{1}\mu_{2}=(\mu_{1}-\mu_{2})^{2}$, which is always positive.

In summary, the homogeneous equilibrium $(x^{*},y^{*})$ can be unstable in the following cases:

• •

if $\operatorname{Tr}(\bm{J}_{R})<0$ and provided $\Lambda^{(\alpha)}<-\frac{\operatorname{Tr}(\bm{J})}{\operatorname{Tr}{\bm{J}_{R}}}$, for at least one choice of $\alpha$.

• •

if $\det(\bm{J_{\alpha}})<0$ which, if $\bm{J}$ is invertible, is verified in the following alternative cases:

• –

if $\det(\bm{J}_{R})<0$ and provided $\Lambda^{(\alpha)}<\Lambda_{-}$, for at least one choice of $\alpha$

• –

if $\det(\bm{J}_{R})>0$, $\operatorname{Tr}(\bm{J}^{-1}\bm{J}_{R})>0$ and provided $\Lambda_{-}<\Lambda^{(\alpha)}<\Lambda_{+}$, for at least one choice of $\alpha$,

where $\Lambda_{\pm}$ identifies the points where the curve $H(\Lambda^{(\alpha)})$ is equal to zero, namely:

[TABLE]

.3 Stationary state and Laplacian eigenvectors

The Laplacian spectrum plays a pivotal role in setting the conditions for the emergence of the patterns, as shown in the previous section. Hereafter, we will show that the spectral properties of the Laplacian convey important information on the characteristics of the patterns that are established in the asymptotic regime. In order to solve the perturbed equations we have indeed expanded the perturbations on the Laplacian eigenvectors basis, thus obtaining

[TABLE]

The perturbations are therefore projected along independent directions which are biased to reflect the peculiar structure of the network architecture, analytically associating to each direction a growth/damping rate given by the respective entry of the dispersion relation. The stationary state will be mainly shaped, at the leading order of approximation, by the most unstable mode, the one that displays the larger dispersion relation entry, as it can be appreciated from figure 4.

.4 Network topology affecting pattern formation: the case of the coupled Stuart-Landau equations.

As an addition to the analysis carried out in the main body of the paper, we here get back to considering the system of coupled Stuart-Landau oscillators and elaborate on the role of long-range interactions for the onset of the instability. To this end we start from a lattice with $k-1$ nearest neighbors and perturb it with the inclusion of remote links assigned with a given probability. In particular, the generic element $A_{ij}$ is one with probability $|i-j|^{-\gamma}$, where $|\cdot|$ measures the shortest distance on the circle between node $i$ and node $j$. The rewiring is implemented is such a way that the average connectivity is preserved. The exponent $\gamma$ quantifies the range of the interaction, $\gamma_{c}=1$ setting the threshold between short and long ranged systems, in one dimension.

The network displayed in the panel of fig. 5(a) is a ring with $k=7$, including self-edges. The associated dispersion relation is characterized by a collection of unstable modes that populate the region at small $|\Lambda^{(\alpha)}|$, where the instability localizes. The emerging patterns are plotted in Fig. 5(d). For $\gamma=2$, the process of network generation is biased to short range couplings, see Fig. 5(b), and the ensuing patterns are depicted in Fig. 5(e). When reducing $\gamma$, the spectral gap (the distance from zero of the first negative eigenvalue of the Laplacian) opens. For sufficiently small values of $\gamma$ ($\gamma=0.6<\gamma_{c}$, in Fig. 5(c)), all Laplacian eigenvalues are associated to negative entries of the dispersion relation. The fully synchronous state is hence stable and the imposed perturbation gets re-absorbed after a transient in time, see figs. 5(f).

As a further check, we consider the Stuart-Landau oscillators coupled via a weighted adjacency matrix that implements the above strategy for interpolating from a short to a long-ranged scenario. The assigned weights scale as a power law of the nodes distance, with exponent $\gamma$. Stated differently, the dependence on the inter-nodes distance lies in the intensity of the connections. As clearly illustrated in Fig. 6, the spectral gap increases when progressively strengthening the connections among distant nodes (i.e. by decreasing $\gamma$), and the system gains a stable configuration characterized by the synchronous oscillators. By adding long-range couplings, one favors a mean-field homogenization of the dynamics, the system behaving hence as a tightly connected uniform component.