Turing-like patterns in an asymmetric dynamic Ising model

TL;DR

This paper introduces a finite-range, asymmetric multi-spin Ising model that can generate Turing-like patterns with tunable length scales, inspired by reaction-diffusion systems and applicable to complex interacting agent systems.

Contribution

It presents a novel Ising model with finite-range, asymmetric interactions capable of producing Turing-like patterns, bridging reaction-diffusion phenomena and spin systems.

Findings

Model can generate Turing patterns with tunable length scales

Asymmetric interactions enable pattern formation without infinite-range interactions

Applicable to complex systems beyond genetic regulation

Abstract

To investigate novel aspects of pattern formation in spin systems, we use a mapping between reactive concentrations in a reaction-diffusion system and spin orientations in a dynamic multiple-spin Ising model. While pattern formation in Ising models always rely on infinite-range interactions, this mapping allows us to design an finite-range interactions Ising model that can produce patterns observed in reaction diffusion systems including Turing patterns with a tunable typical length scale. This model has asymmetric interactions and several spin types coexisting at a site. While we use the example of genetic regulation during embryo-genesis to build our model, it can be used to study the behavior of other complex systems of interacting agents.

Figures

Tables

Table 1. Table 1: Correspondence between parameters

SRDA		ADIM
Diffusion constant	$\sigma$	$r$	Interaction range
Reaction constant	$\overline{\varphi }$	$\overline{J}$	Interaction strength
Reaction constant with gradient	$\kappa$	$k$	Coupling with linear field
Activation threshold	${\theta }_0$	$h_0$	Homogeneous external field
Degradation	$\epsilon =0.5$
Noise		$T$=1	Temperature

Equations

$${\partial }_t{c}_a=D_a{\nabla }^2{c}_a+R_a(\mathbf{c},g)$$

$$P_{ia}^{\mathrm{prod}}=\frac{1}{1+e^{-u_a({\mathbf{C}}_i,g_i)}}$$

$$u_a({\mathbf{c}}_i,g_i)=\sum _b({\varphi }_{ba}{c}_{ib})+{\kappa }_a{g}_i-{\theta }_0$$

$$c_{ia}(t^{\prime \prime })=\sum _{|j-i|<2{\sigma }_a}{G}_{{\sigma }_a}(|j-i|)c_{ja}(t^{\prime })$$

$$c_{ia}(t+1)=\epsilon c_{ia}(t^{\prime \prime })$$

$$E=-J\sum _i\sum _{j\in {\partial }_i}{S}_i\cdot S_j-\sum _i{h}_i{S}_i,$$

$$E=-\sum _i\sum _{j\in {\partial }_i}{{\mathbf{S}}_{\mathrm{j}}}^t\bar{J}{\mathbf{S}}_{\mathrm{i}}-\sum _i{h}_i\mathbf{K}{\mathbf{S}}_{\mathrm{i}}$$

$$h_{ia}^{\mathrm{eff}}(t)=\sum _b\left(\frac{J_{ba}}{V_{{\partial }_{ib}}}\sum _{j\in {\partial }_{ib}}{s}_{jb}(t)\right)+k_a{h}_i-h_0.$$

$$P(s_{ia}(t+1))=\frac{e^{-\beta h_{ia}^{\mathrm{eff}}(t)s_{ia}(t+1)}}{2\cosh (\beta h_{ia}^{\mathrm{eff}}(t))},$$

$${\tilde{s}}_a(q_x,q_y)=\sum _{x=0}^{n-1}\sum _{y=0}^{n-1}{e}^{-\frac{2i\pi }{n}(x{q}_x+y{q}_y)}{s}_a(x,y)$$

$$h_{ia}^{\mathrm{eff}}(t)=\sum _b\left(\frac{J_{ba}}{V_{{\partial }_{ib}}}\sum _{j\in {\partial }_{ib}}{s}_{jb}(t)\right)+k_a{h}_i-h_0,$$

$$P(s_{ia}(t+1)=1)=\frac{e^{\beta h_{ia}^{\mathrm{eff}}(t)}}{2\cosh (\beta h_{ia}^{\mathrm{eff}}(t))}$$

$$T_{SA}(t)=T_0\times a^t$$

$$⟨s_{ia}⟩=\frac{1}{1+e^{\beta ⟨h_{ia}^{\text{eff}}⟩}}$$

$$\langle s_{ia}\rangle=\frac{1}{1+e^{\beta\big{[}\frac{J_{aa}}{V_{\partial_{ia}}}\sum_{j\in\partial_{ia}}\langle s_{ja}\rangle+k_{a}h_{i}-h_{0}\big{]}}}$$

$$\langle s_{ia}\rangle=\frac{1}{1+e^{\beta\big{[}J_{aa}\langle s_{ia}\rangle+k_{a}h_{i}-h_{0}\big{]}}}$$

$$c_{ia}(t^{\prime })=c_{ia}(t)+\frac{1}{1+e^{-u_a({\mathbf{c}}_i,g_i)}}.$$

$$c_{ia}(t+1)=\varepsilon\bigg{[}\sum_{|j-i|<2\sigma_{a}}G_{\sigma_{a}}(|j-i|)[c_{ja}(t)+\frac{1}{1+e^{-u_{a}(\mathbf{c}_{j})}}]\bigg{]}$$

$$\sum _{|j-i|<2{\sigma }_a}{G}_{{\sigma }_a}(|j-i|)[c_{ja}(t)+\frac{1}{1+e^{-u_a({\mathbf{c}}_j)}}]=c_{ia}(t^{\prime })$$

$$c_{ia}(t+1)=\varepsilon\bigg{[}c_{ia}(t)+\frac{1}{1+e^{-u_{a}(\mathbf{c}_{i},g_{i})}}\bigg{]}$$

$$c_{ia}^{\ast }=\frac{1}{1+e^{-u_a({\mathbf{C}}_i)}}$$

$$c^{*}_{ia}=\frac{1}{1+e^{\big{[}\phi_{aa}\frac{c^{*}_{ia}(1-\varepsilon)}{\varepsilon}+\kappa_{a}g_{i}-\theta_{0}\big{]}}}$$

$$\displaystyle\beta\big{[}J_{aa}\langle s_{ja}\rangle+k_{a}h_{i}-h_{0}\big{]}=\frac{\phi_{aa}\varepsilon}{1-\varepsilon}c^{*}_{ia}+\kappa_{a}g_{i}-\theta_{0}.$$

$$⟨s_{ia}⟩=\frac{\epsilon }{1-\epsilon }{c}_{ia}^{\ast }.$$

$$\epsilon =0.5$$

Full text

Pattern Formation, Self organization, Turing, Morphogenesis, Spin, Ising

I Introduction

The mechanisms underlying global pattern formation from local interactions are studied in many fields of physics, chemistry and biology Castets et al. (1990); Koch and Meinhardt (1994); Diego et al. (2018). Reaction diffusion is among possibly the most widespread patterning model. It uses continuous variables such as the concentrations of chemical species that can form complex patterns: e.g. stripes, propagating fronts and oscillatory waves Turing (1952). Another class of pattern formation models are Ising models. Ising variables are two-states (up and down) discrete spins interacting on a lattice. One can distinguish two types of interactions: infinite and finite-range. Infinite-range are more commonly referred to as long-range interactions, i.e. decreasing with the distance as a power law $1/r^{\alpha}$ with $\alpha$ larger or equal to the spatial dimension. Infinite and finite-range interactions show differences in terms of behavior and are investigated using specific computational methods Angelini et al. (2014). In terms of patterning, finite-range interactions mostly lead to ferromagnetic (all spins up or down) or anti-ferromagnetic (alternating up and down spins) states. In the presence of frustration, one-site width stripes can also occur Jin et al. (2013). On the other hand, infinite-range interactions can lead to the formation of more complex patterns. The combination of ferromagnetic nearest neighbor interactions with anti-ferromagnetic infinite-range interactions found for instance in ultra-thin magnetic films Fonseca et al. (2012); De’Bell et al. (2000) is known to produce patterns such as mazes or bubbles. Up to now, infinite-range interactions were considered as a necessary condition for the formation of patterns with a typical width of several lattice constants. Our goal in the present paper is to derive a variant of the Ising model which produces such patterns without infinite-range interactions.

Reaction-diffusion and Ising models are both well-adapted to describe emergent collective behaviors from individual evolution rules. Yet, they have been confronted in few restricted cases. An example is the rigorous proof of the existence of travelling fronts in a ferromagnetic Ising model under Glauber dynamics using a mapping of the magnetization evolution equation with the reaction diffusion Allan-Cahn equation De Masi et al. (1994, 1995).

In this paper, we build upon the similarities between reaction-diffusion and Ising models to construct a dynamic finite-range Ising model giving rise to complex patterns. We introduce in Sec. III a Stochastic Reaction Diffusion Automaton (SRDA) on a lattice, inspired by embryogenesis (Sec. II). Then we construct in Sec. IV our Asymmetric Dynamic Ising Model (ADIM) and establish a mapping between the ADIM and the SRDA in Sec. V. From this mapping we show in Sec. VI that the ADIM can reproduce the variety of patterns observed in reaction-diffusion systems even if it only relies on finite-range interactions.

II An aside on embryogenesis

One of the most emblematic example of pattern formation is embryogenesis, which motivated Turing to develop his well known reaction-diffusion model in 1952 Turing (1952). The major question in this field is to understand how initially identical cells can express different genes to form various tissues and organs at the proper time and location during the development of the embryo.

Turing proposed that interactions between diffusing gene products (called here after species) could lead to the spatial organization of gene expression. In the case of the so-called Turing patterns, two species (an activator $a$ and an inhibitor $b$) can self-organize in stripes, dots, or maze-like patterns. This is driven by three main ingredients: an auto-catalysis of $a$, an asymmetry in reciprocal interactions (i.e. $a$ enhances the formation of $b$ while $b$ penalizes the formation of $a$), and a quantitative difference in their diffusion coefficients. Besides these self-organized patterns Cooke and Zeeman (1976); Corson et al. (2017); Cotterell et al. (2015) suggested by Turing, embryo-genesis provides examples of externally-driven patterns that are pre-formed at large scales, for instance by chemical gradients in the case of the Drosophila’s fertilized egg Wolpert (1969); Gregor et al. (2007). These gradients activate specific genes in different locations of the embryo and the products of these genes subsequently act on other genes in a combinatorial manner. In the current view, the spatial regulation of gene expression in developing embryos is explained by combinations of these two different and complementary classes of scenarios Green and Sharpe (2015).

Various models implementing these two classes have been proposed, either based on differential equations or on cellular automata Sharpe (2017). Recently, an Ising model Hillenbrand et al. (2016) has also been developed to model gene patterning during embryo-genesis. Each spin corresponds to a gene that can be in one of the two states: active or inactive. Each site represent a nucleus and can contain several spins/genes. The spatial proximity between nuclei defines the interaction lattice (non zero couplings), while coupling values defines the interaction network between spin/genes. These couplings are a simplification of the molecular mechanisms at work: gene transcription in RNA, RNA translation in proteins and protein diffusion to neighboring nuclei where they can modify the transcription rate of other genes. While being different from more detailed reaction-diffusion models, this Ising model features the most important ingredient: short range interactions between spins/genes Hillenbrand et al. (2016). Based on this similarity, we also placed ourselves in this context to build an extended variant of the Ising model, implementing asymmetric interactions between spin types, that could lead to Turing patterns. To achieve this goal, we start by constructing a reference reaction diffusion model in the next section.

III Stochastic Reaction Diffusion Automaton (SRDA)

Reaction diffusion models describe the evolution in time and space of the concentrations $\mathbf{c}=(c_{a},c_{b},\dots)$ of $n_{s}$ chemical species $a$, $b$,…undergoing two processes: diffusion, associated to diffusion coefficients $D_{a},D_{b},\dots$, and local reactions $R_{a}(\mathbf{c},g),R_{b}(\mathbf{c},g),\dots$ resulting in the creation of units of $a,b,\dots$ per time unit, where $g$ is a possible external gradient forcing the system. The continuous space and time equations are for each specie $a$ of the form:

[TABLE]

As a discrete space and time model of reaction-diffusion, we use here a probabilistic automaton on a 1D lattice. The vector $\mathbf{c}_{i}(t)=(c_{ia}(t),c_{ib}(t),\dots)$ contains the concentrations of the $n_{s}$ considered species on each lattice site $i$ at discrete time $t$. The external gradient $g$ is a prescribed time-independent additional specie of concentration $g_{i}$ on site $i$. Each time step is divided in three sub-steps. The first one ($t\to t^{\prime}$) is the production event: each concentration is incremented by 1 with probability:

[TABLE]

It is a common usage to model gene regulation by sigmoidal functions. This exact form has been used for Drosophila development modelingIlsley et al. (2013). We choose $u_{a}(\mathbf{c}_{i},g_{i})$ of the form:

[TABLE]

where the $\phi_{ba}$ are the reaction constants between species (the sum is over all species, including $a$ itself), $\kappa_{a}$ represents the effects of the gradient $g$ on specie $a$ and $\theta_{0}$ is an activation threshold common to all species.

The second event ($t^{\prime}\to t^{\prime\prime}$) is diffusion, approximated by a Gaussian kernel $G_{\sigma_{a}}$ of standard deviation $\sigma_{a}$ and mean value 0:

[TABLE]

The third event ($t^{\prime\prime}\to t+1$) is degradation, assumed to occur at the same rate $\varepsilon$ for all species:

[TABLE]

where $0<\varepsilon<1$.

To represent the steady states, concentrations are normalized by $c_{\textrm{max}}=\frac{\varepsilon}{1-\varepsilon}$, so $c^{*}=\frac{c}{c_{\textrm{max}}}$ is comprised between 0 and 1.

As described in Sc. II, two classes of patterns can emerge from our model : gradient-induced patterning (Fig 1-i) and Turing instabilities (Fig 1-ii). The two types of mechanisms can be combined (Fig 1-iii) leading to the appearance of hybrid patterns Green and Sharpe (2015).

IV Asymmetric Dynamic Ising Model (ADIM)

The classical Ising model has been created in 1920 Brush (1967) as a toy model describing ferromagnetism. The spins are bivalued (usually $S_{i}=\pm 1$, but here, we equivalently choose $S_{i}=0,1$). They are placed on a lattice and interact with their nearest neighbors with an interaction constant $J>0$ that tends to align them. In the presence of a space-dependent magnetic field $h_{i}$, the energy of the system is:

[TABLE]

where $\partial_{i}$ contains the nearest neighbors of site $i$.

The first step to map our reaction-diffusion automaton is to design a multiple-spin Ising model, similarly to Hillenbrand et al. (2016). To each site $i$ of the lattice, we associate $n_{s}$ spin types corresponding to the different species $a$, $b$, …. $S_{i}$ becomes a vector $\mathbf{S}_{\rm{i}}=(s_{ia},s_{ib},\dots)$ of components [math] or $1$. The concentration of the SRDA now corresponds to the average value of the spin, named magnetization. Each spin type interacts with all other spin types including itself. We thus rewrite $J$ as a matrix $\bar{J}$ whose coefficient $J_{ab}$ represents the interactions between $a$-spins and $b$-spins. To mimic the gradient effect on the different genes, $h_{i}$ is supposed to be coupled differently to each spin-type and we introduce a coupling vector $\mathbf{K}=(k_{a},k_{b},\dots)$, where $k_{a}$ represents the effect of $h_{i}$ on the $a$-spins. Eq (6) becomes:

[TABLE]

To get patterns typical of reaction-diffusion, we also introduce spin type dependent interaction range $r_{a}$. An $a$-spin now acts on its neighboring sites up to a distance $r_{a}$, including itself. We denote by $\partial_{ia}$ this set of spins (e.g. $\partial_{ia}=[i-r_{a}:i+r_{a}]$ on a 1D chain). The interaction parameters are re-scaled using the volume $V_{\partial_{ia}}$, which is the number of sites in $\partial_{ia}$ (e.g. $V_{\partial_{ia}}=2r_{a}+1$ in 1D).

We notice that only the average values $\frac{J_{ab}+J_{ba}}{2}$ are important in the determination of the equilibrium state of Eq (7), obtained at finite temperature $T\neq 0$ using e.g. the Metropolis algorithm N. Metropolis and S. Ulam (1949). However, the matrix $\bar{J}$ needs to be effectively asymmetric to reproduce the asymmetry of the reaction constants $\phi_{ab}$ leading to the formation of interesting patterns. We thus take a step further from Hillenbrand et al. (2016) by implementing parallel dynamics in our model, similarly to kinetic asymmetric Ising models Mézard and Sakellariou (2011); Roudi and Hertz (2011), making it a non-equilibrium model.

At each time step $t$, we calculate an effective field $h^{\rm{eff}}_{ia}(t)$ influencing the $s_{ia}$ spin:

[TABLE]

where we introduce $h_{0}$ an homogeneous external gradient.

All spins are then updated to give the configuration at time $t+1$ according to the following probability distribution:

[TABLE]

where $\beta=1/T$. The temperature $T$ accounts for the noisiness of the system.

V Model comparison

The SRDA and ADIM are both governed by Markovian dynamics and share important features: asymmetry of the coupling matrices and locality of the interactions (see App. A for the detailed numerical implementation). In the ADIM however, local reactions and transport are merged in nearest neighbors interactions. Still, many parameters appear to have similar roles: $\bar{J}$ and $\Phi$, $k$ and $\kappa$, $h_{0}$ and $\theta_{0}$, and $r$ and $\sigma$, as presented Table 1. Two parameters are nevertheless model specific: the degradation rate $\varepsilon$ and the temperature $T$.

We wish to investigate the similarity in patterns that can be obtained by confronting the two models. We study the case $n_{s}=1$, which corresponds to a classical Ising model under a space-dependent external field. We fix $T=1$ in the ADIM since this temperature gives the most direct equivalence with Eq (2) (see App. B.). For the sake of simplicity, we fix $r_{a}=\sigma_{a}=1$ so that only between nearest neighbors interact and we choose the space-dependent external field to be linear. In this conditions, we show by a calculation using the homogeneous solutions in mean-field approximation (see App. B) that the degradation rate that gives the best mapping between interaction parameters in both models is $\varepsilon_{\textrm{opt}}=0.5$.

Fig 2-a-d presents the mean concentration $\langle c^{*}\rangle$ and the magnetization $M$ of the patterns obtained by the SRDA and ADIM in the planes ($\phi_{aa}$, $\kappa_{a}$) and ($J_{aa}$, $k_{a}$) with $h_{0}=\theta_{0}=1$. They represent the fraction of time during which $a$ is being produced or spins $s_{a}$ are up, averaged over each positions. For all probed regions of the parameter space, very similar patterns are observed for both the ADIM and the SRDA. Two broad zones of $\langle c^{*}\rangle=M=1$ and $\langle c^{*}\rangle=M=0$ correspond to full activation (Fig 2-$\blacktriangle$) and full inhibition (Fig 2-$\blacktriangledown$). Transition between theses two states can be done either through the formation and displacement of a sharp boundary (Fig 2-$\blacksquare$), corresponding to strong inhibition by the external gradient and auto-activation, or by the appearance and intensification of a smooth gradient-like pattern (Fig 2-$\blacklozenge$), corresponding to a gradient activation balanced by a null or very low auto-inhibition.

Yet, the sharpness of the $\blacksquare$-pattern is different in both models. Intuitively this sharpness depends on the correlation length of the system. The difference in patterns thus reflects a quantitative difference between the roles of the diffusion coefficient and the interaction range. Fig 2-e represents the evolution of the $\blacksquare$-pattern in SRDA and ADIM as a function of respectively on $\sigma_{a}$ and $r_{a}$. The maps confirm that $\sigma_{a}$ and $r_{a}$ are the spatial parameters that control the boundary sharpness. In both models, increasing these parameters causes an increase in the correlation lengths and thus a decrease of the border sharpness. It is therefore possible to choose a value for $r_{a}$ for which the slope exactly matches the slope obtained for any given value of $\sigma_{a}$ (see App. C for more discussion on this matter).

VI Turing patterns in the ADIM

In light of these similarities between SRDA and ADIM in the case $n_{s}=1$, we now look at the case $n_{s}=2$ with no external gradient and try to reproduce Turing patterns. To investigate a richer variety of Turing patterns, we also switch to 2D with periodic boundary conditions. A main feature enabling Turing instabilities in reaction diffusion models is the difference of diffusion coefficients for both species. Likewise, the ratio ${r_{b}}/{r_{a}}$, where $a$ is the activator spin-type and $b$ is the inhibitor spin-type, needs to be large enough to form patterns of different spin orientations. 2D patterns as a function of the interactions $J_{ab}>0$ and $J_{ba}<0$ are presented in Fig 3-a for $r_{a}=1$ and $r_{b}=5$, on a $128\times 128$ square lattice. A subtle trade-off between the values of the interactions gives rise to different types of patterns, including ”maze” $\bigstar$ and ”bubbles” $\bullet$. These two patterns are used for further analysis in the next part.

To quantitatively study the 2D patterns obtained with our ADIM, we compute the discrete 2D Fourier Transform (FT) of $s_{a}$ on our periodic square lattice of linear size $n$:

[TABLE]

These patterns have one characteristic length resulting in a circular distribution in the reciprocal space (Fig 4a). The FT gives two quantities: the pseudo order parameter $\Lambda$ and the characteristic length $L_{\Lambda}$. We average $\tilde{s}_{a}$ over circles of radius $q$ to get the mean intensity $\lambda(q)$ and define the order parameter $\Lambda=\lambda(q_{\textrm{peak}})$ as the peak value of these curves. The corresponding characteristic length in the direct space is $L_{\Lambda}=\frac{1}{q_{\textrm{peak}}}$. For illustration, Fig S2 shows the $\lambda(q)$ curves at different temperatures for both ”maze” and ”bubbles” patterns.

We first investigate the role of $r_{b}$ in the emergence of patterns for a fixed $r_{a}=1$. Fig 4a shows the ”maze” pattern and its FT for different values of $r_{b}$. While there is no pattern for small values of $r_{b}$, maze patterns appear for higher values. Computing $L_{\Lambda}/r_{a}$ for $r_{b}$ ranging from 1 to 25 we conclude that the maze pattern appears for $r_{b}\geq 3$ (Fig 4b). The characteristic length of the pattern increases linearly with $r_{b}$ with a coefficient close to 1. For higher values of $r_{a}$, $L_{\Lambda}/r_{a}$ increases linearly with the ratio $r_{b}/r_{a}$ with the same coefficient. We conclude that if one takes $r_{a}$ as a unit length, the behavior of adimensional quantities $L_{\Lambda}/r_{a}$ and $r_{b}/r_{a}$ is independent of the lattice spacing. When the lattice spacing becomes small compared to the characteristic length $r_{a}$ we can refine the value of the ratio $r_{b}/r_{a}$ for which patterns appear. In this case we find that this value tends towards 2.

We next investigate the effect of temperature on the stability of these patterns using the pseudo order parameter $\Lambda$ defined above. In infinite-range interacting Ising models at equilibrium, phase transitions between smectic, nematic or liquid phases can occur when $T$ varies Abanov et al. (1995). Here, the various patterns do not break any symmetry and no such phase transition can occur. Fig 5 shows $\Lambda$ for different system sizes and different temperatures. A dynamic crossover occurs between a patterned and a disordered phase for both patterns, albeit at different temperatures. The crossovers do not depend on the lattice size, thus excluding potential finite size effects.

In the Ising model, the temperature reflects thermal noise, whereas in a reaction diffusion system, noise is related to the number of molecules: a lower number leads to larger relative fluctuations in the local concentrations. Molecular dynamics studies point out the importance of fluctuations for the emergence of Turing patterns Lemarchand and Nowakowski (2011). Our results on the stability of patterns with temperature suggest on the other hand that if the number of molecules becomes too small (corresponding to high temperatures, hence high fluctuations) the patterns could disappear.

VII Conclusion

To summarize, using a embryo-genesis-inspired mapping between an Ising model and a reaction-diffusion automata, we have constructed the ADIM, a finite-range out-of-equilibrium Ising model that can give rise to Turing patterns with a typical length scale of several lattice spacings. Such Turing patterns had previously only been observed in infinite-range Ising models. It is worth noticing that our variant of the Ising model can be seen as a cellular automata featuring rules similar to the ones which can be derived from reaction-diffusion equations (see for instance Weimar and Boon (1994)). Our approach shows that a discrete toy model with different ranges of interaction and an asymmetry in the interaction coupling between species is sufficient to produce Turing-like patterns.

Appendix A Numerical methods

ADIM— The system is composed of $N_{S}$ sites. On each site is associated $n_{s}$ spin types. Parallel dynamics is implemented as follows.

We initialize the system at $t=0$ by randomly attributing a value [math] or $1$ at each of the $n_{s}\times N_{S}$ spins. At each time step $t$ we calculate the effective field $h^{\rm{eff}}_{ia}(t)$ for each $a$-type spin on each site $i$:

[TABLE]

where $J_{ba}$ is the action of $b$-spins on $a$-spins, $h_{i}$ the value of the space dependent external magnetic field at site $i$, $k_{a}$ the coupling between $\mathbf{h}$ and $a$-type spins and $h_{0}$ a global external field. $\partial_{ib}$ designates the nearest neighbors of spin $s_{ib}$, defined by the interaction range $r_{b}$ for $b$-spin type such as $\partial_{ib}=\{j,|j-i|<r_{b}\}$ and $V_{\partial{ib}}=Card(\partial_{ib})$.

We calculate the corresponding Boltzmann probabilities:

[TABLE]

where $\beta=T^{-1}$ is the inverse temperature.

In 1D, the simulation runs at a fixed temperature $T$. In 2D, the stationary state is calculated using a simulated annealing (SA) to avoid freezing in configurations with long relaxation times.Cooling strategy is exponential:

[TABLE]

with $0<a<1$.

The spins at time $t+1$ are simultaneously chosen according to these probabilities. The operation is repeated during $n_{t}$ time steps, and in case of SA, $a$ is chosen so that $T_{SA}(n_{t})=T$, ie a=\big{(}\frac{T}{T_{0}}\big{)}^{n_{t}^{-1}}. Then the simulation is pursued and recorded for $n_{\textrm{rec}}$ time steps. Finally $\langle s_{ia}\rangle$ is obtained by averaging over the $n_{\textrm{rec}}$ recorded configurations. Typically, $n_{t}=250$ and $n_{\textrm{rec}}=50$.

SRDA— The system is composed of $N_{S}$ sites at which is associated $n_{s}$ specie concentrations $c_{a}$. Each concentration is randomly chosen at $t=0$ between [math] and $1$. Each time step is subdivided in 3 events as presented in the main text.

This is repeated for $n_{t}$ time steps to reach the stationary state. Then each concentration is normalized by $c_{\textrm{max}}$ to obtain $0<c^{*}_{ia}<1$. Typically, $n_{t}=200$.

In both models, the boundary conditions are set to reflecting in the case 1D to have gradient-dependent patterns. In the 2D ADIM, we choose periodic boundary conditions.

Appendix B Determination of $\varepsilon_{\textrm{opt}}$ using homogeneous stationary solutions

We use here the homogeneous steady state solutions in mean-field approximation for both models to establish a relationship between parameters in the case of one specie $a$ ($n_{s}=1$), and more specifically to find the optimal value $\varepsilon_{\textrm{opt}}$ for which both model lead to the same patterns for equal network parameters.

ADIM— The mean-field approximation for the ADIM is (see Eq.(9)):

[TABLE]

For only one spin type $a$ ($n_{s}=1$):

[TABLE]

Since we consider a homogeneous state and we are in 1D, $\sum_{j\in\partial_{ia}}\langle s_{ja}\rangle=(2r_{a}+1)\langle s_{ia}\rangle$ and $V_{\partial_{ia}}=2r_{a}+1$, leading to:

[TABLE]

SRDA— In the mean-field approximation for the SRDA model, we average over several realizations of time evolution for each time step. Thus, we replace the probabilistic increment of the production step (Eq (2)) by a deterministic one:

[TABLE]

By combining this new equation to the two other events (diffusion (Eq (4)) and degradation (Eq (5)) ), we get:

[TABLE]

In the classical PDE reaction-diffusion approach as defined in Eq (1), searching for an homogeneous steady state equates to look for the solution of $R_{a}(\mathbf{c},g)=0$ since $\partial_{t}c_{a}=0$ (steady state) and $\nabla^{2}c_{a}=0$ (homogeneous). In our SRDA, the homogeneity translates in $\forall j$ $c_{ja}(t)=c_{ia}(t)$, and using the fact that $\sum_{r<2\sigma}G_{\sigma}(r)=1$, we obtain:

[TABLE]

which leads to a simplification of Eq (S8) :

[TABLE]

Being at steady-state $c_{ia}(t+1)=c_{ia}(t)$, and using $c_{max}=\frac{\varepsilon}{1-\varepsilon}$, we retrieve $c^{*}$ :

[TABLE]

For $n_{s}=1$:

[TABLE]

Optimal mapping value of $\varepsilon$— Eq (S12) and Eq (S4) are combined to search for the solution of $\langle s_{ia}\rangle=c^{*}_{ia}$, leading to:

[TABLE]

For $\beta=T^{-1}=1$, identical network parameters ($J_{aa}=\phi_{aa}$, $k_{a}=\kappa_{a}$ and $h_{0}=\theta_{0}$) and under the same external gradient ($h_{i}=g_{i}$), we get:

[TABLE]

Finally, we obtain the optimal degradation parameter $\varepsilon$ for which both models gives equivalent homogeneous steady-state for equivalent parameters, as presented in the results of the main content of this letter Fig 2 :

[TABLE]

Appendix C Equivalence between $\sigma_{a}$ and $r_{a}$

To evaluate the correspondence between $\sigma_{a}$ in the SRDA and $r_{a}$ in the ADIM we use the sharp boundary pattern presented in the main text for $n_{s}=1$ under a linear external gradient. This pattern appear when the spin/gene activates itself but is repressed by the gradient ($J_{aa}=7$ and $k_{a}=-5$). As seen in the main text and shown again in Fig S1-a, as $\sigma_{a}$ or $r_{a}$ increases, the sharpness of the boundary decreases. We first measure this sharpness by fitting the boundary interface ($0.2<c_{a}^{*},\langle s_{a}\rangle<0.8$ ) by a affine function of the position (Fig S1-b). The slope defines the sharpness of the boundary. Both the diffusion constant $\sigma_{a}$ and the interaction range $r_{a}$ as a function of sharpness are fitted by power laws $y=\frac{a}{x^{b}}$ (Fig S1-c). From these fits we calculated $\sigma_{a}^{\textrm{fit}}$ and $r_{a}^{\textrm{fit}}$ for any given sharpness and plotting $r_{a}^{\textrm{fit}}$ as a function of $\sigma_{a}^{\textrm{fit}}$ gives us a equivalence between $r_{a}$ and $\sigma_{a}$ (Fig S1-d).