Stochastic fluctuations and quasi-pattern formation in reaction-diffusion systems with anomalous transport

TL;DR

This paper develops a Gaussian approximation for intrinsic noise in reaction-subdiffusion systems, revealing how stochastic fluctuations and anomalous transport can induce and enhance pattern formation.

Contribution

It introduces a novel analytical framework to incorporate intrinsic noise into fractional reaction-diffusion equations, advancing understanding of noise-driven pattern formation in subdiffusive systems.

Findings

Subdiffusion can promote both deterministic and stochastic pattern formation.

Intrinsic noise reduces the diffusion ratio threshold for pattern emergence.

The combined effect of subdiffusion and stochasticity enhances pattern formation more than either factor alone.

Abstract

Many approaches to modelling reaction-diffusion systems with anomalous transport rely on deterministic equations and ignore fluctuations arising due to finite particle numbers. Starting from an individual-based model we use a generating-functional approach to derive a Gaussian approximation for this intrinsic noise in subdiffusive systems. This results in corrections to the deterministic fractional reaction-diffusion equations. Using this analytical approach, we study the onset of noise-driven quasi-patterns in reaction-subdiffusion systems. We find that subdiffusion can be conducive to the formation of both deterministic and stochastic patterns. Our analysis shows that the combination of subdiffusion and intrinsic stochasticity can reduce the threshold ratio of the effective diffusion coefficients required for pattern formation to a greater degree than either effect on its own.

Equations

$$\frac{\partial P_{i,t}}{\partial t}=\sum _{i^{\prime }}\left({\varphi }_{i,i^{\prime }}-{\delta }_{i,i^{\prime }}\right)t_0^{-\gamma }{}_0{D}_t^{1-\gamma }{P}_{i,t},$$

$${}_0{D}_t^{1-\gamma }f\left(t\right)=\frac{1}{\Gamma \left(\gamma \right)}\frac{\partial }{\partial t}{\int }_0^t\frac{f\left(t^{\prime }\right)}{{\left(t-t^{\prime }\right)}^{1-\gamma }}d{t}^{\prime },$$

$${\mathcal{L}}_t\left\{{}_0{D}_t^{1-\gamma }f\left(t\right)\right\}\left(u\right)=u^{1-\gamma }\hat{f}\left(u\right),$$

$${\varphi }_{i,i^{\prime }}=\{\begin{array}{ll}\frac{1}{3},& \mathrm{if}i=i^{\prime }+1\\ \frac{1}{3},& \mathrm{if}i=i^{\prime }\\ \frac{1}{3},& \mathrm{if}i=i^{\prime }-1\end{array}.$$

$${\psi }^{\alpha }\left(t\right)$$

$$\varnothing$$

$$2X_A+X_B$$

$$X_A$$

$$\varnothing$$

$$X_A$$

$$4X_A+X_B$$

$$n_{i,\tau +\Delta ,t+\Delta }^{\alpha }-n_{i,\tau ,t}^{\alpha }=$$

$$-\sum _{i^{\prime }}\frac{k_{i^{\prime },i,\tau ,t}^{\left(H\right)\alpha }}{\Delta },$$

$$n_{i,0,t+\Delta }^{\alpha }=\sum _r\frac{{\ell }_{i,r,t}^{\left(R\right)}{\nu }_r^{\alpha }\theta \left({\nu }_r^{\alpha }\right)}{\Delta }+\sum _{i^{\prime }}\frac{k_{i,i^{\prime },\tau ,t}^{\left(H\right)\alpha }}{\Delta },$$

$$Z\Big{[}\{\Xi_{i,\tau,t}^{\alpha}\}\Big{]}=\left\langle{\exp\left(i\sum_{i,\alpha,\tau,t}\Xi^{\alpha}_{i,\tau,t}x^{\alpha}_{i,\tau,t}\right)}\right\rangle_{\{x^{\alpha}_{i,\tau,t}\}},$$

$$\frac{\partial n_{i,\tau ,t}^{\alpha }}{\partial t}+\frac{\partial n_{i,\tau ,t}^{\alpha }}{\partial \tau }$$

$$n_{i,0,t}^{\alpha }$$

$$p_{i,t}^{\alpha }{n}_{i,\tau ,t}^{\alpha }$$

$${\gamma }_{i,t}^{\alpha }$$

$$\frac{\partial n_{i,t}^{\alpha }}{\partial t}=\sum _{i^{\prime }}\left\{\left({\varphi }_{i,i^{\prime }}-{\delta }_{i,i^{\prime }}\right){\left(t_0^{\alpha }\right)}^{-{\gamma }^{\alpha }}{e}^{-{\int }_0^t{p}_{i^{\prime },T^{\prime }}^{\alpha }d{T}^{\prime }}{}_0{D}_t^{1-{\gamma }^{\alpha }}\left[n_{i^{\prime },t}^{\alpha }{e}^{{\int }_0^t{p}_{i^{\prime },T^{\prime }}^{\alpha }d{T}^{\prime }}\right]\right\}+f_{i,t}^{\alpha }+{\eta }_{i,t}^{\alpha },$$

$${\eta }_{i,t}^{\alpha }$$

$$\displaystyle-e^{\int^{t}_{0}p_{i^{\prime},T^{\prime}}^{\alpha}dT^{\prime}}\left(t_{0}^{\alpha}\right)^{-\gamma^{\alpha}}{}_{0}D_{t}^{1-\gamma^{\alpha}}\Bigg{(}\int_{0}^{t}\Psi^{\alpha}\left(\tau\right)e^{\int^{t-\tau}_{0}p_{i^{\prime},T^{\prime}}^{\alpha}dT^{\prime}}\int^{\tau}_{0}\frac{\xi^{\alpha}_{i^{\prime},T,T+t-\tau}}{\Psi^{\alpha}\left(T\right)e^{-\int_{0}^{T}p_{i^{\prime},T^{\prime}+t-\tau}^{\alpha}dT^{\prime}}}dTd\tau\Bigg{)}\Bigg{]}\Bigg{\}}$$

$$+{\int }_0^t{\xi }_{i,\tau ,t}^{\alpha }d\tau +{\xi }_{i,0,t}^{\alpha }.$$

$$\frac{\partial n_{i,t}^{\alpha }}{\partial t}=\sum _{i^{\prime }}\left[\left({\varphi }_{i,i^{\prime }}-{\delta }_{i,i^{\prime }}\right){\left(t_0^{\alpha }\right)}^{-1}{n}_{i^{\prime },t}^{\alpha }\right]+f_{i,t}^{\alpha }+{\eta }_{i,t}^{\alpha }.$$

$${\eta }_{i,t}^{\alpha }={\int }_0^t{\xi }_{i,\tau ,t}^{\alpha }d\tau +{\xi }_{i,0,t}^{\alpha }.$$

$$n_{i,t}^{\alpha }={\bar{n}}_{i,t}^{\alpha }+{\delta }_{i,t}^{\alpha }.$$

$${\underline{\hat{\tilde{\delta }}}}_{u,q}={\underline{\underline{\hat{\tilde{m}}}}}_{u,q}{\underline{\widehat{\tilde{ϵ}}}}_{u,q},$$

$${\theta }_{\gamma }^{\mathrm{act}}=\frac{{\left(t_0^{\left(\mathrm{act}\right)}\right)}^{\gamma }}{t_0^{\left(\mathrm{inh}\right)}{\left({\bar{p}}^{\left(\mathrm{act}\right)}\right)}^{1-\gamma }},$$

$${\theta }_{\gamma }^{\mathrm{inh}}=\frac{t_0^{\left(\mathrm{act}\right)}{\left({\bar{p}}^{\left(\mathrm{inh}\right)}\right)}^{1-\gamma }}{{\left(t_0^{\left(\mathrm{inh}\right)}\right)}^{\gamma }},$$

$$x_{i,\tau +\Delta ,t+\Delta }^{\alpha }-x_{i,\tau ,t}^{\alpha }-\frac{\Delta }{N}{f}_{i,\tau ,t}^{\alpha }\left(\{x_{i,\tau ,t}^{\alpha }\}\right)-\frac{\sqrt{\Delta }}{N}{\chi }_{i,\tau ,t}^{\alpha }$$

$$x_{i,0,t+\Delta }^{\alpha }-\frac{1}{N}{f}_{i,0,t}^{\alpha }\left(\{x_{i,\tau ,t}^{\alpha }\}\right)-\frac{{\chi }_{i,0,t}^{\alpha }}{N\sqrt{\Delta }}$$

$$C^{\alpha,\alpha^{\prime}}_{i,i^{\prime},\tau,\tau^{\prime},t,t^{\prime}}\equiv\big{<}\chi^{\alpha}_{i,\tau,t}\chi^{\alpha^{\prime}}_{i^{\prime},\tau^{\prime},t^{\prime}}\big{>}.$$

$$\frac{\partial n_{i,\tau ,t}^{\alpha }}{\partial \tau }+\frac{\partial n_{i,\tau ,t}^{\alpha }}{\partial t}$$

$$n_{i,0,t}^{\alpha }$$

$$Z\Big{[}\{\Xi_{i,\tau,t}^{\alpha}\}\Big{]}=\int\prod_{i,\alpha,\tau,t}\Big{\{}dx^{\alpha}_{i,\tau,t}\Big{\}}\left[\exp\left(i\sum_{i,\alpha,\tau,t}\Xi^{\alpha}_{i,\tau,t}x^{\alpha}_{i,\tau,t}\right)\mathcal{P}\left(\{x^{\alpha}_{i,\tau,t}\}\right)\right]$$

Full text

I Introduction

Reaction-diffusion schemes are a well-established tool for modelling non-linear pattern-forming phenomena in a wide variety of systems, ranging from developmental biology Gierer and Meinhardt (1972); Murray (1981a, b) to chemical reactions Petrov et al. (1997), from the self-organisation of slime mold Falcke and Levine (1998) and the distribution of mussels Van de Koppel et al. (2008) to Rayleigh-Bénard convection in fluids Rayleigh (1916). Most notably, Turing showed in his seminal work Turing (1952), that systems which have a stable fixed point when diffusion is not present could exhibit an instability in the presence of diffusion. He showed that diffusion, normally considered a stabilising, equilibrating influence, combined with activator-inhibitor reaction dynamics, could provide a simple mechanism for pattern formation.

A difficulty that arises in the Turing interpretation of pattern formation is that the diffusion constants of the activator and inhibitor species are often required to be quite disparate in order for patterns to be found Murray (2004). These large ratios of diffusion constants are unphysical.

One potential solution to this is that a system may not have to undergo a deterministic Turing instability in order to exhibit pattern formation. Stochastic quasi-patterns, which arise due to the intrinsic noise of individual-based systems, may be sufficient to explain the spatial ordering observed in some natural systems Butler and Goldenfeld (2009, 2011). The emergence of such quasi-patterns requires a lower threshold ratio of diffusion coefficients than the emergence of deterministic patterns Biancalani et al. (2017, 2010). Noise-induced quasi-patterns have been observed recently in experimental biological systems Karig et al. (2018); Di Patti et al. (2018).

Another proposed solution to this problem is modelling one or more of the components as subdiffusing Metzler and Klafter (2000, 2004); that is, diffusing with a mean squared displacement which is sub-linear in time. This has also been shown to reduce the ratio of the effective diffusion coefficients required for pattern formation Yadav et al. (2008); Baron and Galla (2019). Furthermore, subdiffusion has been observed in many experimental systems with low particle numbers, such as proteins moving in the cell membrane Ghosh and Webb (1994); Schwille et al. (1999) and in other transport phenomena involving obstacles or binding sites Saxton (1994, 1996). It has been argued that morphogens could be subject to similar binding and trapping effects and move subdiffusively as a result Hornung et al. (2005); Yuste et al. (2010); Fedotov and Falconer (2014).

However, as of yet, the treatment of reaction-subdiffusion systems has been restricted to deterministic equations which ignore the stochastic effects due to the finite numbers of particles involved. Examples of the deterministic treatment include Fedotov and Falconer (2014, 2013); Yuste et al. (2010). Stochastic effects due to demographic noise are important particularly in biological systems, where subdiffusion has been observed, and where the particle numbers involved can be sufficiently small for noise to be non-negligible. From a more theoretical perspective, subdiffusion is a non-Markovian phenomenon. Because of this, reaction-diffusion systems with subdiffusing components provide a good opportunity to study the combination of memory effects and intrinsic noise.

In this work, we use an individual-based approach, which explicitly takes into account intrinsic stochasticity, to study systems with reactions and anomalous diffusion. We carry out a path-integral (or generating-functional) analysis of the stochastic dynamics, combining transport and reactions, and perform an expansion in the strength of the intrinsic noise. To lowest order, effectively neglecting fluctuations, we recover the familiar deterministic reaction-subdiffusion equations. Taking into account sub-leading orders of fluctuations, we find additional coloured Gaussian noise terms in these equations, encapsulating the intrinsic stochasticity of the individual-based system. We show that the analytical expressions for this noise can be used to characterise the emergence of stochastic quasi-patterns in reaction systems with subdiffusing components. Specifically, we use both the Brusselator Glansdorff and Prigogine (1971) and Lengyel-Epstein Lengyel et al. (1990); Lengyel and Epstein (1991) systems as examples. We conclude that the combination of the two effects, noise-driven quasi-patterns and subdiffusion, can lower the critical ratio of diffusion coefficients required for pattern formation significantly, and to a greater degree than either effect on its own.

The remainder of this paper is set out as follows: In Section II, we describe the individual-based model in detail and briefly present some background material related to subdiffusion and to the Brusselator and Lengyel-Epstein models. Section III contains an outline of our analysis of the individual-based system and the generating-functional approach to approximating the noise; further details of this calculation can be found in the Supplemental Material. In Section IV, we show that this approach can be used to derive the familiar reaction-diffusion equations with fractional diffusion, now with additional noise terms. In Section V, we describe how one can then compute the fluctuations about the deterministic solution, using the linear-noise approximation. We verify our approach by comparing our theory predictions to stochastic simulations. In Section VI, we go on to use our theoretical approach to find the parameter regions where stochastic patterns are present and where deterministic patterns are present. Finally, in Section VII we discuss our results and conclude.

II Model construction and background

In this section we introduce our model and the general notation, and we briefly summarise some background material. We also define the Brusselator and Lengyel-Epstein models, which we later use to illustrate the formation of quasi-patterns in stochastic subdiffusive systems.

II.1 Individual-based model

We consider a general class of individual-based reaction-diffusion models. We refer to the individuals as ‘particles’ from here on in, but these particles may represent biological entities or molecules of chemical reactants.

Several species of particle react with one another and hop around on a discrete lattice. We use $\alpha$ to index the different species. Each lattice site can be occupied by multiple particles simultaneously, and we denote the number of particles of type $\alpha$ at site $i$ at time $t$ by $n^{\alpha}_{i,t}$. Reactions between particles occur locally within a given site of the lattice, and result in the annihilation and creation of particles, or the conversion of one type of particle into another. For the purposes of simplicity, we will confine ourselves to a discrete one-dimensional lattice with periodic boundary conditions. Most of our analysis could be generalised to other spatial arrangements or to higher dimensions, and indeed the continuum limit may be taken.

Particles are assigned a random waiting-time at birth (creation), drawn from a, as yet unspecified, waiting-time probability density function $\psi^{\alpha}\left(\tau\right)$. These may be different for the different species of particle, as indicated by the superscript. The particle then hops to a new location once it has waited for the assigned time, assuming that it has not been eliminated in a reaction. A new waiting-time is drawn from $\psi^{\alpha}\left(\tau\right)$ once the particle has hopped and the process begins again. This means that a particle which has been at its current site for an amount of time $\tau$ hops with a rate $h^{\alpha}_{\tau}=\frac{\psi^{\alpha}\left(\tau\right)}{\Psi^{\alpha}\left(\tau\right)}$, where $\Psi^{\alpha}\left(\tau\right)=1-\int^{\tau}_{0}\psi^{\alpha}\left(\tau^{\prime}\right)d\tau^{\prime}$ is the survival probability. The quantity $h^{\alpha}_{\tau}$ is often referred to as the hazard rate Metzler and Klafter (2000); we use subscript notation for $\tau$ for later convenience. When the particle does hop, its new location is drawn from the hopping kernel $\phi_{i,i^{\prime}}$. Here $\phi_{i,i^{\prime}}$ is the probability that the particle hops from location $i^{\prime}$ to $i$, given that a hopping event occurs. For our purposes, $\phi_{i,i^{\prime}}$ will be a function of $\lvert i^{\prime}-i\rvert$ only, so as to ensure the translational invariance of the problem. In principle the hopping kernel could also be different for the different species. In order to keep the model simple we use the same hopping kernel for all species; it is straightforward to extend the model to the more general case.

During the particles’ sojourn periods at a given site, they may undergo reactions. We index the various reaction types with $r$. The rate $\lambda_{i,r,t}$ at which reactions of type $r$ occur at site $i$ is in general dependent on the number of particles of the various types currently at site $i$. We denote the number of particles of type $\alpha$ that are produced or annihilated in a reaction of type $r$ by $\nu^{\alpha}_{r}$, which can be positive or negative. The constants $\nu^{\alpha}_{r}$ are the so-called stoichiometric coefficients McNaught and Wilkinson (1997). If a particle is annihilated in a reaction then the hop which was scheduled to occur for that particle no longer occurs. As a result of this and of the fact that the reaction rates depend on local concentrations, the hopping process and the reactions are interdependent and cannot easily be separated as in conventional reaction-diffusion equations (such as those used by Turing Turing (1952)).

II.2 Subdiffusion

We first discuss the phenomenon of subdiffusion in the case where particles undergo hopping events but no reactions. Particles are said to undergo subdiffusion if the mean-squared displacement for a single particle behaves as follows in the long term: $\langle x(t)^{2}\rangle\sim t^{\gamma}$, where $0<\gamma<1$. The behaviour approaches normal diffusion as $\gamma\to 1$. Subdiffusive transport may be produced by choosing a long-tailed waiting-time distribution $\psi\left(t\right)$ for particles which hop around as described in the previous section (for the time being we suppress the dependence on the particle species $\alpha$). The distribution we will use in this paper to model subdiffusion is that of Mittag-Leffler, which produces the desired behaviour and is particularly convenient for the theoretical analysis. That is, we choose $\psi\left(t\right)=-\frac{d}{dt}E_{\gamma}\left[-\left(\frac{t}{t_{0}}\right)^{\gamma}\right]$ where $E_{\gamma}\left[x\right]$ is the Mittag-Leffler function Mainardi and Gorenflo (2000). The parameter $t_{0}$ sets the overall scale of the hopping process. This distribution has the convenient property that its Laplace transform is given by $\hat{\psi}\left(u\right)=\frac{1}{1+\left(t_{0}u\right)^{\gamma}}$. Notably, for $\gamma=1$ the Mittag-Leffler function reduces to an exponential such that $\psi\left(t\right)=\frac{1}{t_{0}}e^{-\frac{t}{t_{0}}}$.

For a single species of diffusing particle, which undergoes no reactions, one may use the Mittag-Leffler waiting-time distribution, along with the Montroll-Weiss formula Michael F. Shlesinger (2017); Montroll and Weiss (1965), to recover the fractional diffusion equation Klages et al. (2008); Metzler and Klafter (2000); Mendez et al. (2010)

[TABLE]

where $P_{i,t}$ is the probability of finding a particular such particle at position $i$ at time $t$. The Riemann-Liouville fractional derivative ${}_{0}D_{t}^{1-\gamma}$ is defined as

[TABLE]

where $\Gamma(\cdot)$ denotes the gamma function. The fractional derivative has the property

[TABLE]

where we use $\mathcal{L}_{t}$ to denote the Laplace transform with respect to $t$. One can show from Eq. (1) that $\langle i^{2}\rangle=\frac{2}{3\Gamma\left(1+\gamma\right)}\left(\frac{t}{t_{0}}\right)^{\gamma}$, if one uses the symmetrical hopping kernel

[TABLE]

In our system, different species of particles may hop with different typical waiting times $t_{0}^{\alpha}$ and different anomalous exponents $\gamma^{\alpha}$. Precisely, the waiting-time distributions for each species are given by

[TABLE]

II.3 Reaction schemes

As example systems, we will use the Brusselator Glansdorff and Prigogine (1971) and the Lengyel-Epstein Lengyel et al. (1990); Lengyel and Epstein (1991) models. The Brusselator is an activator-inhibitor model, first conceived of to describe oscillatory chemical reactions such as the Belousov-Zhabotinsky reaction Tyson (1976). The reactions between the two species of particle involved, $A$ and $B$, are given by

[TABLE]

Without hopping the system has a homogeneous stable deterministic fixed point at $\bar{n}^{\left(A\right)}=aN,\bar{n}^{\left(B\right)}=\frac{b}{a}N$ so long as $b<1+a^{2}$, where $\bar{n}^{\left(\alpha\right)}$ is the number of particles of type $\alpha$ at the fixed point.

The Lengyel-Epstein model was introduced primarily as a way of modelling the $\textrm{ClO}_{2}^{-}$-$\textrm{I}^{-}$-MA reaction, which exhibits Turing patterns experimentally Castets et al. (1990). We use a simplified two-species version of the full model, which effectively assumes that the concentrations of the remaining components are constant. The reactions are given by

[TABLE]

where the rate $R=\frac{1}{\left(n^{\left(A\right)}\right)^{3}}\frac{cN}{dN^{2}+\left(n^{\left(A\right)}\right)^{2}}$ is dependent on the concentration of type-A particles. The Lengyel-Epstein system has a homogeneous deterministic fixed point at $\bar{n}^{(A)}=\frac{aN}{5b},\bar{n}^{(B)}=\frac{bdN}{c}\left[1+\left(\frac{\bar{n}^{(A)}}{\sqrt{d}N}\right)^{2}\right]$ so long as $ca>\frac{3}{5}a^{2}-25b^{2}d$.

In both models, species $A$ is the activator and species $B$ is the inhibitor. The parameter $N$ in both systems characterises the typical number of particles per site and will become useful to us later when we perform a system-size expansion in order to analyse noise in the these systems. The reaction rates (the number of reactions which occur per unit time) are calculated according to the usual mass action kinetics Érdi and Tóth (1989). For example, the rate at which the last reaction in Eq. (7) occurs is $R\times\left(n^{\left(A\right)}\right)^{4}\times n^{\left(B\right)}=\frac{cNn^{\left(A\right)}n^{\left(B\right)}}{dN^{2}+\left(n^{\left(A\right)}\right)^{2}}$.

III Approximation of the fluctuations in particle number using Gaussian noise

III.1 Formulation of the problem

In order to capture the noise-driven effects in the models described in section II.1, we carry out an expansion in the inverse system-size. The idea is similar to the principles underlying the Kramers–Moyal expansion or the system-size expansion by van Kampen for Markovian systems Van Kampen (1992). However, there are also conceptual differences due to the non-Markovian nature of subdiffusion.

To carry out the expansion, we use the reciprocal of the parameter $N$ in Eqs. (6) and (7); $N$ characterises the typical particle number per site. More specifically, ensemble-averaged particle numbers are of the order $N$ at each lattice site, whereas the fluctuations in particle numbers are of the order $\sqrt{N}$. Thus, for large $N$, the noise is small in relative terms. The limit $N\to\infty$ reproduces the deterministic behaviour, and when $N$ is large but finite the expansion can be expected to accurately describe stochastic corrections. By performing the expansion, we obtain a set of stochastic differential equations (SDEs) Gardiner (2009) with Gaussian noise terms. These SDEs encapsulate not only the deterministic trajectory of the system but also the next-order stochastic noise corrections.

A complicating feature of systems where the waiting-time distribution for the hops is non-exponential is that the hazard rate for hopping is non-constant, meaning that what happens at a given point in time is dependent on the history of the system. A sensible and common way to deconvolute the problem is to introduce the age coordinate $\tau$, which denotes the length of time that a particle has resided at a particular location since its last hop Cox (2017). That is, we denote the number of particles of type $\alpha$ which have resided at site $i$ for a time between $\tau_{1}$ and $\tau_{2}$ by $\int_{t_{1}}^{t_{2}}n^{\alpha}_{i,\tau,t}d\tau$. This effectively recasts the problem as Markovian. The quantity $n^{\alpha}_{i,\tau,t}$ is a density of particles per time $\tau$.

In order to develop the formalism we discretise time into steps of size $\Delta$. We assume that all reaction rates and hopping rates remain constant during each step, similar to the $\tau$-leaping approach to Gillespie simulations in discrete time Gillespie (2001). Eventually, we will take the limit $\Delta\to 0$ to restore continuous time.

For our system, the quantity $n^{\alpha}_{i,\tau,t}$ may change due to two effects in each time step: reactions and hopping. We represent the changes due to these effects by the following quantities respectively: $k^{\left(R\right)\alpha}_{i,r,\tau,t}$ and $k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}$. That is, $k^{\left(R\right)\alpha}_{i,r,\tau,t}$ is the number of particles of type $\alpha$ and age $\tau$ at position $i$ which are annihilated in the time step from $t$ to $t+\Delta$ due to reactions of type $r$. We note that newly created particles have age zero; this will be dealt with separately below. Similarly, $k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}$ is the number of particles of age $\tau$ hopping away from position $i^{\prime}$ to $i$ at time $t$. The sets of integer variables $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}$ and $\{k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}\}$ are stochastic, that is they take different values for every realisation of the system.

We then have

[TABLE]

where $\theta\left(\cdot\right)$ denotes the Heaviside function. We note that, in the discrete-time setup, the total number of particles at site $i$ at time $t$ is given by $n^{\alpha}_{i,t}=\Delta\sum_{m}n^{\alpha}_{i,\tau=m\Delta,t}$; the quantity $n^{\alpha}_{i,\tau,t}$ therefore has dimensions of inverse time, in-line with the expression on the right-hand side of Eq. (8).

Eq. (8) only includes reactions which annihilate particles, i.e. it describes the outflux of particles from position $i$. When particles are produced, or when they newly arrive at a location after hopping, they have age $\tau=0$. Therefore, the influx of particles to position $i$ is given by

[TABLE]

where $\ell_{i,r,t}^{\left(R\right)}$ is the number of reactions of type $r$ firing in the time window $t$ to $t+\Delta$ at position $i$. We note also that $\theta\left(-\nu_{r}^{\alpha}\right)\ell^{\left(R\right)}_{i,r,t}\lvert\nu_{r}^{\alpha}\rvert=\sum_{\tau}k^{\left(R\right)\alpha}_{i,r,\tau,t}$. Furthermore, particles of a given species which are annihilated due to a reaction at a given site are selected at random from all particles of this species at the site, irrespective of age.

III.2 Generating functional approach to the system-size expansion

To approximate the stochastic fluctuations in particle number with Gaussian noise, one would normally write down a master equation and expand in powers of $N^{-1}$ in order to obtain a Fokker-Planck equation, see for example Risken (1996); Gardiner (2009); McKane and Newman (2005); Biancalani et al. (2010). From this Fokker-Planck equation one could read off the corresponding SDE. Since the hopping of the particles in our system is history-dependent, there is no straightforward way to write down a master equation, even after the introduction of the additional age coordinate $\tau$. The most natural way to analyse a non-Markovian stochastic system is with the Martin-Siggia-Rose-Janssen-De Dominicis (MSRJD) path integral, which takes into account all possible histories of the system Martin et al. (1973); Janssen (1976); De Dominicis (1976). Such an approach allows one to perform the expansion in inverse powers of the system-size without writing down a master equation explicitly Brett and Galla (2013, 2014).

In order to reflect the system-size dependence in our calculation we introduce $x^{\alpha}_{i,\tau,t}=\frac{n^{\alpha}_{i,\tau,t}}{N}$. These quantities are of order $N^{0}$. The MSRJD generating functional then takes the form of a path integral over all possible trajectories of the variables $\{x^{\alpha}_{i,\tau,t}\}$ Altland and Simons (2010). It can be written as

[TABLE]

where $\left\langle{\dots}\right\rangle_{\{x^{\alpha}_{i,\tau,t}\}}$ denotes an average over all trajectories of the system. The $\{\Xi^{\alpha}_{i,\tau,t}\}$ are source variables. The procedure for performing the expansion in $N^{-1}$ is similar to that used in Brett and Galla (2013, 2014): We find the joint probability distribution for the sets of variables $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}$ and $\{k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}\}$ and rewrite the generating functional in term of these quantities. We then average these random numbers against their joint distribution, to obtain the generating functional in terms of only the coordinates $\{x^{\alpha}_{i,\tau,t}\}$ and the model parameters. We then carry out an expansion in $N^{-1}$ up to and including sub-leading order. This approximate generating functional is recognised as that of an effective SDE. The leading-order terms in this SDE correspond to the deterministic dynamics; after further re-arrangement it reproduces the well-known deterministic reaction-subdiffusion equation, as shown in Sec. IV. This deterministic approximation is accurate in the limit $N\to\infty$. The next-order terms correspond to Gaussian noise corrections, with a standard deviation of order $N^{-1/2}$.

The details of this procedure are provided in the Supplementary Material (Sections S1 and S2); here we only quote the final result, in which we have restored continuous time. It is given by the following stochastic equations

[TABLE]

In these expressions, $p_{i,t}^{\alpha}$ is the per capita removal rate for particles of species $\alpha$ at position $i$ and time $t$; $\gamma_{i,t}^{\alpha}$ is the total production rate for the particles of type $\alpha$ at $i$ and $t$, that is

[TABLE]

The quantities $\lambda_{i,r,t}$, $\nu_{r}^{\alpha}$, $\phi_{i,i^{\prime}}$ and $h_{\tau}^{\alpha}$ are defined in Section II.1. The fraction $n_{i,\tau,t}^{\alpha}/n_{i,t}^{\alpha}$ in the first of the relations in Eq. (12) reflects the fact that the particles which are to be removed are selected irrespective of their age $\tau$. The quantities $\{\xi^{\alpha}_{i,\tau,t}\}$ in Eq. (11) represent white Gaussian noise of zero mean. Eqs. (11) are approximations of the full random process given by Eqs. (8) and (9), in which the randomness is discrete. While the noise $\{\xi^{\alpha}_{i,\tau,t}\}$ is white, we note that the components are correlated across species and lattice sites. The expressions for these correlations are somewhat lengthy; we derive and give them in Section S2 of the Supplemental Material. It is also important to note that the noise is multiplicative, i.e. the elements in the covariance matrix of the $\{\xi^{\alpha}_{i,\tau,t}\}$ depend on the variables $\{n^{\alpha}_{i,\tau,t}\}$. This is similar to the outcome of a Kramers–Moyal expansion of the master equation for conventional Markovian systems Gardiner (2009); Risken (1996).

IV Fractional reaction-diffusion equation with noise

Using the waiting-time distributions discussed in Section II.2 in conjunction with Eqs. (11), we are able to obtain the fractional reaction-diffusion equation reported in the literature (see e.g. Yadav et al. (2008); Vlad and Ross (2002); Yuste et al. (2010); Fedotov and Falconer (2014)) if we neglect the noise terms. Including these terms, we are able to capture effects driven by the stochasticity of the original individual-based dynamics.

One proceeds from Eqs. (11) by integrating out the age variables $\tau$, with the aim of finding a time-evolution equation in terms of only $n_{i,t}^{\alpha}$ and the noise. Using the Mittag-Leffler function from Eq. (5) as the waiting-time distribution, one arrives at the following fractional reaction-diffusion equation for the particle number $n_{i,t}^{\alpha}$

[TABLE]

where we have defined the total reaction rate $f^{\alpha}_{i,t}=\gamma^{\alpha}_{i,t}-p^{\alpha}_{i,t}n^{\alpha}_{i,t}$.

In the limit of continuous space, the sum involving the hopping kernel becomes a Laplacian operator (i.e. $\sum_{i^{\prime}}\left(\phi_{i,i^{\prime}}-\delta_{i,i^{\prime}}\right)\to\sigma^{2}\nabla^{2}$, with a suitable constant $\sigma^{2}$, related to the variance of hopping distances). Thus, the fractional reaction-diffusion equation Eq. (13) corresponds to Eq. (3) in Yadav et al. (2008), but with the addition of a noise term. We recover the deterministic reaction-subdiffusion equation if we take the infinite system-size limit, $N\to\infty$, whereupon the noise becomes negligible.

One notes that the reaction and diffusion terms in Eq. (13) are coupled. This is due to the non-Markovian nature of the hopping and the fact that particles may be annihilated in reactions and are thus unable to perform scheduled hops. The exponential factors in Eqs. (13) correspond to the probability of particles surviving without being annihilated in a reaction.

The noise also involves such exponential factors and is given by

[TABLE]

The derivation of Eqs. (13) and (14) from Eqs. (11) broadly follows the method in Vlad and Ross (2002) or Fedotov and Falconer (2012) but also handles the additional noise terms. It is given in the Supplemental Material (Section S3). The noise variables $\{\eta^{\alpha}_{i,t}\}$ involve integrals of the $\{\xi^{\alpha}_{i,\tau,t}\}$, and as a consequence they are correlated in time, in addition to their correlations across components and lattice sites. Given that the statistics of $\{\xi^{\alpha}_{i,\tau,t}\}$ depend on the variables $\{n^{\alpha}_{i,\tau,t}\}$, this noise is multiplicative as well.

Although the expression for the noise $\eta^{\alpha}_{i,t}$ in Eq. (14) appears cumbersome at first glance, a great deal of simplification can be achieved in the regime where one considers small deviations about the homogeneous fixed point of the deterministic system. We discuss this in the next section.

One notes that for $\gamma\to 1$ we recover the normal reaction-diffusion equation where the reactions and diffusion are uncoupled:

[TABLE]

In this case, the expression for the noise also simplifies greatly and contains no exponential factors:

[TABLE]

This can be seen from Eq. (14) by using the fact that the Mittag-Leffler function reduces to an exponential in the limit $\gamma\to 1$ and that the fractional derivative becomes an identity operator. In this special case, $\langle\xi^{\alpha}_{i,\tau,t}\xi^{\alpha}_{i^{\prime},\tau^{\prime},t^{\prime}}\rangle\propto\delta\left(t-t^{\prime}\right)$, and we then also have $\langle\eta^{\alpha}_{i,\tau,t}\eta^{\alpha}_{i^{\prime},\tau^{\prime},t^{\prime}}\rangle\propto\delta\left(t-t^{\prime}\right)$. Therefore, in the Markovian limit, we recover white noise. This is explored further in the Supplemental material (Section S2).

V Linear-noise approximation and comparison with simulations

V.1 Linear-noise approximation

We now derive explicit expressions for the deviations of the stochastic system (i.e., with a finite number of particles per lattice site) from the solution of the deterministic reaction-subdiffusion equation.

We write $\bar{n}^{\alpha}_{i,t}$ for the solution of Eq. (13) with the noise term removed, and define the deviation from this deterministic solution through

[TABLE]

We focus on fluctuations about the steady state, i.e. we assume that $\bar{n}^{\alpha}_{i,t}\equiv\bar{n}^{\alpha}$ is the homogeneous fixed point of the deterministic system. From the results stated so far, we are able to compute the correlation matrix of the fluctuations $\langle\delta^{\alpha}_{i,t}\delta^{\alpha^{\prime}}_{i^{\prime},t^{\prime}}\rangle$, which we then compare with results from stochastic simulations of the individual-based system.

The expansion in powers of $N^{-1}$ in Section III was based on the assumption that the fluctuations about the deterministic trajectory were small in comparison to the total numbers of particles. A more precise formulation of this assumption allows us to make further simplifications. Since the noise term $\eta^{\alpha}_{i,t}$ in Eq. (13) is of order $N^{1/2}$, it is also reasonable that the fluctuations would be such that $\delta^{\alpha}_{i,t}\sim\mathcal{O}(N^{1/2})$. This is similar to the assumption used by van Kampen Van Kampen (1992) to perform the system-size expansion now named after him. Using this, one can expand Eq. (13) about the deterministic fixed point $\bar{n}^{\alpha}$.

We note that the $\delta^{\alpha}_{i,t}$ are small deviations about the homogeneous deterministic fixed point. The correlators $\langle\delta^{\alpha}_{i,t}\delta^{\alpha^{\prime}}_{i^{\prime},t^{\prime}}\rangle$ therefore do not describe any any deterministic pattern-forming features of the system. Instead, they characterise stochastic phenomena induced by fluctuations due to finite system sizes $N$.

To leading order (i.e., sending $N\to\infty$) we obtain the deterministic dynamics already mentioned, and in which all noise terms are removed. The terms to sub-leading order result in a linear expression for the fluctuations $\delta^{\alpha}_{i,t}$ in terms of $\eta^{\alpha}_{i,t}$. Within this order of approximation the quantities dependent upon $n^{\alpha}_{i,t}$ in the noise correlators $\langle\eta^{\alpha}_{i,t}\eta^{\alpha}_{i^{\prime},t^{\prime}}\rangle$ are evaluated at the fixed point $\bar{n}^{\alpha}$. The noise in the dynamics for the $\{\delta^{\alpha}_{i,t}\}$ is therefore now additive rather than multiplicative. The procedure is discussed further in Section S4 of the Supplemental Material.

Once the linearised equations are found, the problem is most naturally handled by taking Fourier and Laplace transforms with respect to the spatial and temporal coordinates respectively. We use $\mathcal{F}_{n}\left\{f_{n}\right\}\left(q\right)=\sum_{n}e^{inq}f_{n}=\tilde{f}_{q}$ to denote the discrete Fourier transform and $\mathcal{L}_{t}\left\{g_{t}\right\}\left(u\right)=\int_{0}^{t}e^{-ut}g\left(t\right)dt=\hat{g}_{u}$ to denote the Laplace transform. We then obtain equations of the form

[TABLE]

with suitable Gaussian noise variables $\{\hat{\tilde{\epsilon}}^{\alpha}_{u,q}\}$. In the time domain this noise is coloured. We use underscores to indicate vectors in the space of species. The square matrix $\underline{\underline{\hat{\tilde{m}}}}_{u,q}=(\hat{\tilde{m}}_{u,q})^{\alpha\beta}$ has dimension equal to the number of particle species. Its precise form is given in the Supplemental Material (Section S4), as are full expressions for the correlators $\langle\tilde{\epsilon}^{\alpha}_{q,t}\tilde{\epsilon}^{\alpha^{\prime}}_{q^{\prime},t^{\prime}}\rangle$.

It then remains to invert the Fourier and Laplace transforms. Expressions for the resulting correlators $\langle\delta^{\alpha}_{i,t}\delta^{\alpha^{\prime}}_{i^{\prime},t^{\prime}}\rangle$ are given in the Supplemental Material (Section S5). In practice, the evaluation of the correlators $\langle\delta^{\alpha}_{i,t}\delta^{\alpha^{\prime}}_{i^{\prime},t^{\prime}}\rangle$ is performed using numerical inverse Laplace transform methods (such as that of Zakian Halsted and Brown (1972)) and fast Fourier transform routines. This is because analytical expressions for the inverse transforms of the elements of $\underline{\underline{\hat{\tilde{m}}}}_{u,q}$ are in general difficult to find. That being said, the equal-time equal-wavenumber correlator, $\langle\tilde{\delta}^{(\alpha)}_{q,t}\tilde{\delta}^{(\alpha^{\prime})\star}_{q,t}\rangle$, can indeed be found without the use of a numerical inverse Laplace transform technique (see Section S5 of the Supplemental Material).

V.2 Comparison of theory to simulation results

The theory hitherto presented can be expected to be accurate when fluctuations about the deterministic homogeneous fixed point are small enough so that the linear approximation is valid. The examples in Figs. 1 and 2 show that one obtains good agreement with individual-based simulations of the full system with typical particle numbers per lattice site as low as $N\approx 4000$.

Fig. 1 demonstrates the accuracy of the theory for the steady-state power-spectrum of fluctuations $\langle\lvert\tilde{\delta}^{\left(A\right)}_{q,t}\rvert^{2}\rangle$ for the case of the Brusselator model. We compare the results of simulations to our theory for the cases where the activator subdiffuses and the inhibitor diffuses normally, and vice versa. Fig. 2 again compares the theory predictions to Monte Carlo simulations but this time in the case of the Lengyel-Epstein model with a subdiffusing activator and a subdiffusing inhibitor. The auto-correlators $\langle\delta^{\left(A\right)}_{i,t}\delta^{\left(A\right)}_{i^{\prime},t^{\prime}}\rangle$ and $\langle\delta^{\left(B\right)}_{i,t}\delta^{\left(B\right)}_{i^{\prime},t^{\prime}}\rangle$ and the cross-correlator $\langle\delta^{\left(A\right)}_{i,t}\delta^{\left(B\right)}_{i^{\prime},t^{\prime}}\rangle$ are shown for different temporal separations $t-t^{\prime}$ and spatial separations $i-i^{\prime}$. The agreement in both Figs. 1 and 2 is good apart from some minor deviations. These arise from a combination of the stochastic nature of the simulations and the limits to the accuracy of the numerical inverse Laplace transform. The discrepancy reduces as the number of trials, over which the average is taken, is increased.

To carry out the simulations, it was necessary to modify the Gillespie algorithm Gillespie (1976) in a similar way to Anderson (2007). This was due to the non-Markovian nature of the hopping processes. Broadly, the simulations involve keeping a list of the scheduled hopping times of particles, carrying these out at appropriate times and in the right sequence, while performing reactions in the intermediate times. The full procedure is given in Appendix A. To the best of our knowledge, individual-based simulations of reaction-diffusion systems with subdiffusive transport have not been performed in this way previously.

The full expressions for the correlators $\langle\delta^{\alpha}_{i,t}\delta^{\alpha^{\prime}}_{i^{\prime},t^{\prime}}\rangle$, from which the theory lines in Figs. 1 and 2 are derived, can be found in Section S5 of the Supplemental Material. In the Supplement we also explain how the equal-time correlators in Fig. 1 can be found without the use of a numerical inverse Laplace transform.

VI Onset of stochastic and deterministic Turing patterns

VI.1 Description of the phenomena

In a deterministic system which exhibits a Turing instability, a finite range of Fourier modes with non-zero wavenumbers is unstable. That is to say, if the system is perturbed from its uniform fixed point, the amplitudes of these modes will grow with time. A criterion for the Turing instability in systems with one subdiffusing component and one normally diffusing component is given in Yadav et al. (2008); Baron and Galla (2019). In practice, particle numbers do not deviate infinitely far from the fixed point. Instead, the growth of the unstable modes is curtailed by the non-linearity of the reaction equations. If one simulates systems with such an instability and looks at the Fourier transform of deviations from the uniform state one finds the dominant peak at a non-zero wavenumber. It is this wavenumber which characterises the periodicity of the observed Turing patterns. It is important to note that noise is not required to sustain these patterns once the initial perturbation about the uniform state has been applied. In particular, the amplitude of these deterministic patterns is set by the nonlinearities of the reaction-(sub)diffusion system, and not by any source of stochasticity.

However, one may also observe patterns in systems with finite particle number (i.e., systems with intrinsic noise), even when parameters are such that the purely deterministic system would not show any patterns Butler and Goldenfeld (2009, 2011). These patterns are noise-driven, and their amplitude is set by the strength of the noise; in individual-based systems the variance of the intrinsic noise in turn is determined by the inverse typical particle number $N^{-1}$. We refer to patterns formed in this way as quasi-patterns or stochastic Turing patterns. They can be detected in the Fourier spectra of fluctuations about the uniform deterministic fixed point; these spectra show a peak at a characteristic nonzero wavenumber. In contrast to deterministic Turing patterns, this peak is not observed in the absence of noise. Such patterns have been analysed in systems with regular diffusion, see e.g. Butler and Goldenfeld (2009, 2011); Biancalani et al. (2010); Karig et al. (2018).

In Fig. 3 we show examples of Turing patterns in the Lenyel-Epstein system with subdiffusion. The data in the upper panel is from individual-based simulations in the parameter regime in which the deterministic model shows Turing patterns. The pattern is seen in the stochastic simulation as well, but it is important to stress that this is not a noise-driven pattern. The stochasticity in the individual-based model modulates the deterministic pattern, but the amplitude of the structure is independent of the noise, i.e., increasing the number of particles, $N$, per site does not change the relative magnitude of the pattern. In the central panel we show an example of a quasi-pattern. This data is taken in the regime in which the deterministic Lengyel-Epstein system does not have any instability; the pattern is purely noise-driven and its amplitude decreases with increasing particle number per lattice site. One further main difference between stochastic and deterministic patterns is that the stochastic patterns are not stationary. Instead, regions of high concentration for one chemical species continually assemble, shift and dissipate in such a manner that a typical level of periodicity is maintained. Deterministic patterns on the other hand are stationary in time. In the lower panel of Fig. 3 finally, the parameters are such that neither deterministic nor stochastic Turing patterns are seen. The system fluctuates about the deterministic uniform fixed point, but no particular spatial structure emerges in these fluctuations. That being said, the fluctuations do appear to be quite large in magnitude, given the system size. Similarly, the magnitude of the noise-induced pattern in Fig. 3(b) is similar to that of the pattern in panel (a). This is not an unusual observation and is due to a phenomenon known as ‘stochastic amplification’, which has been reported in the context of noise-driven cycles and patterns in normally diffusing systems (see e.g. McKane and Newman (2005); Alonso et al. (2007); Butler and Goldenfeld (2009, 2011); Biancalani et al. (2010)).

VI.2 Phase diagrams

The theory we developed in the earlier sections allows us to find the equal-time correlator in Fourier space (i.e., the structure factor $\langle\tilde{\delta}^{(\alpha)}_{q,t}\tilde{\delta}^{(\alpha^{\prime})\star}_{q,t}\rangle$) of fluctuations without having to carry out time-consuming simulations. We can use these expressions to search parameter space and to identify the regimes in which stochastic patterns emerge, and the parameter ranges where they do not. This gives a more complete picture of the effect of subdiffusion on stochastic pattern formation. Our results are summarised in the phase diagrams in Figs. 4 and 5, which we discuss in more detail below. While the expressions from the theory have to be evaluated numerically, we note that determining the phase behaviour from the analytical expressions is very efficient; accurately establishing the phase diagrams from individual-based simulations on the other hand would require unrealistic computing time.

It is well-known that the emergence of deterministic Turing patterns relies on the slow diffusion of the activator and the comparatively fast diffusion of the inhibitor. In previous work Baron and Galla (2019) we showed that it is possible to define effective diffusion coefficients for systems in which one reactant undergoes subdiffusion and the other normal diffusion; see also Yadav et al. (2008). To characterise the degree to which the rates of (sub)diffusion of the two particle species differ, we define the ratio of these effective diffusion constants, $\theta_{\gamma}$, see again Baron and Galla (2019); Yadav et al. (2008) for details. The detailed definition depends on which one of the components (activator or inhibitor) is subdiffusing. We have

[TABLE]

when the activator subdiffuses, and

[TABLE]

when the inhibitor subdiffuses. The quantity $\bar{p}$ in these expressions denotes the per capita removal rate of the relevant substance at the uniform deterministic fixed point.

Figs. 4 and 5 show the critical ratios of the effective diffusion constants of the activator and the inhibitor for the onset of stochastic and deterministic Turing patterns in the Lengyel-Epstein and Brusselator models respectively. For any fixed value of the anomalous exponent $\gamma$ the following behaviour is observed in Figs. 4 and 5. For low values of $\theta_{\gamma}$ neither deterministic nor stochastic patterns emerge. As the ratio of effective diffusion coefficients is increased and crosses a first critical threshold, $\theta_{\gamma}^{s}$, a phase is entered in which stochastic patterns are found, but where the deterministic system shows no patterns. As $\theta_{\gamma}$ is increased further the system undergoes a full deterministic Turing instability at $\theta_{\gamma}^{d}$.

Additionally, we find that in general, if the subdiffusion of a particular reactant is conducive to the formation of deterministic patterns, then it is also conducive to the formation of stochastic patterns. Similarly if subdiffusion of a reactant makes deterministic pattern formation more difficult then it also reduces the parameter range in which noise-driven patterns are found. In other words, if the critical value $\theta_{\gamma}^{d}$ for the onset deterministic patterns increases (decreases) with $\gamma$, then $\theta_{\gamma}^{s}$ also increases (decreases) with $\gamma$.

One notes that due to the concentration-independent decay rate of the activator in the Brusselator model, both $\theta_{\gamma}^{d}$ and $\theta_{\gamma}^{s}$ do not vary with the anomalous exponent $\gamma$ when the activator is subdiffusing (left-hand panel in Fig. 5). However, if the death rate of the subdiffusing component is not constant, one finds that, in our examples, as $\gamma$ decreases, so do both $\theta_{\gamma}^{d}$ and $\theta_{\gamma}^{s}$ if the activator is subdiffusing (left-hand panel of Fig. 4). Conversely, if the inhibitor is subdiffusing and the inhibitor death rate is concentration-dependent, a reduced value of $\gamma$ makes it more difficult it is for Turing patterns to form (right-hand panels of Figs. 4 and 5).

One caveat to the approach we have presented is that analytical results are only valid for large system sizes $N$. That being said, the qualitative conclusions we have reached (that the combination of stochasticity and subdiffusion lowers the threshold for pattern formation) can be extended to the regime of moderate $N$, where the theory is not as accurate. We have verified this in simulations, and discuss this point briefly in Section 6 of the Supplemental Material.

VII Summary and conclusions

We have successfully extended the description of reaction-diffusion systems with anomalous transport to include the intrinsic noise which comes about when the number of particles in the system is finite. To do this, we carried out a Gaussian approximation of the individual-based system, using a generating-functional approach. We also provided a prescription for finding the fluctuations about the stable homogeneous fixed point of reaction-diffusion systems with anomalous diffusion. This was done in the example cases of the Brusselator and Lengyel-Epstein systems, and our theoretical predictions were successfully tested against computer simulations. Finally, we used this theory to determine the parameter regimes where one could expect to observe stochastic Turing patterns. We found examples for which a subdiffusing activator encouraged both stochastic and deterministic pattern formation, whereas a subdiffusing inhibitor hindered both stochastic and deterministic pattern formation. This is exemplified in Figs. 1, 4 and 5.

The theory developed in this paper is entirely general and would be equally applicable to any reaction scheme with anomalously diffusing reactants; the two examples discussed here were mainly chosen for illustration. Systems such as the Lotka-Volterra Lotka (1909) dynamics, the Oregonator Field and Noyes (1974); Tyson and Fife (1980); Jahnke et al. (1989) or the Schnakenberg model Schnakenberg (1979) with subdiffusion would have been equally valid candidates for the study of noise-driven patterns. We could also have studied systems with different waiting-time distributions or hopping kernels.

A well-known problem with using Turing’s mechanism as an explanation for pattern formation in nature is the often unphysically large ratio of the diffusion constants required for the instability to occur. Both subdiffusion and stochastic pattern formation have been proposed separately as potential remedies for this problem. In this paper, we have shown that a combination of the two can lower the threshold for pattern formation to a greater degree than either mechanism on its own. This is demonstrated by the fact that the critical ratio of the diffusion constants is always lower in Figs. 4 and 5 for the onset of stochastic patterns than it is for deterministic patterns, and by the fact that this critical ratio can reduce as the activator becomes more subdiffusive.

Recently, the existence of noise-induced Turing patterns in bacterial cultures has been examined Karig et al. (2018). However, these experiments were concerned with a mechanism for pattern formation for which the diffusion was presumed to be Markovian. A possibility for further experiment would be to analyse stochastic pattern formation in a system in which the diffusion of the reactants was known to be subdiffusive. Our approach makes predictions for the way in which noise-driven pattern formation is affected by subdiffusion.

Another question for future study would be to ask how other stochastic phenomena in reaction-diffusion systems, such as waves Biancalani et al. (2011), would be affected by anomalous diffusion. Conversely, one might also wonder about how such phenomena are affected by superdiffusion. The theoretical approach we have developed is sufficiently general, and can be applied to the study of such systems.

Acknowledgements. JWB thanks the Engineering and Physical Sciences Research Council (EPSRC) for funding (PhD studentship, EP/N509565/1).

Appendix A Individual-based simulations

The non-Markovian nature of subdiffusive hopping requires us to modify the traditional Gillespie algorithm Gillespie (1976), in a similar way to Anderson (2007). The statistics of the processes laid out in Section II.1 are preserved by using the following algorithm:

Set $t=0$ and initialise the system. (In our simulations, we start from $n^{\left(1\right)}_{i}=\bar{n}^{\left(1\right)}$ and $n^{\left(2\right)}_{i}=\bar{n}^{\left(2\right)}$.) For each particle draw a waiting time until its next hop from the appropriate distribution, and create a list of all scheduled hopping events.

• 2.

Take note of the earliest hopping time and the position of the particle associated with this time; we label the time at which the next hop occurs $\tau^{\left(H\right)}$.

• 3.

Based on the numbers of particles at each position, calculate the reaction rates for each type of reaction and each point on the lattice, $\lambda_{i,r,t}$. Calculate the sum of the reaction rates $\lambda_{t}=\sum_{r,i}\lambda_{i,r,t}$ and draw the time at which the next reaction is due to occur: $\tau^{\left(R\right)}=t-\frac{1}{\lambda_{t}}\ln\left(s\right)$, where $s$ is an independent random number from the uniform distribution over $(0,1]$. This reaction can occur at any position on the lattice.

• 4.

(a) If $\tau^{\left(R\right)}<\tau^{\left(H\right)}$, perform one single reaction. Choose this to be a reaction of type $r$ at site $i$ with probability $\lambda_{i,r,t}/\lambda_{t}$. If particles are removed during this reaction delete the scheduled hopping events for these particles from the list. If particles are created, draw hopping times from the appropriate distribution, and add these to the list of scheduled hops. Set $t=\tau^{\left(R\right)}$. Go to 2.

(b) If $\tau^{\left(R\right)}>\tau^{\left(H\right)}$, perform the hop associated with $\tau^{\left(H\right)}$ in step 2, and discard the reaction event from step 3. To do this draw a random hopping distance from the hopping kernel $\phi_{i,i^{\prime}}$, and re-locate the particle. Draw a new waiting time until the next hop of this particle, and add this time to the list. Set $t=\tau^{\left(H\right)}$. Go to 2.

To draw waiting times from the Mittag-Leffler distribution in Eq. (5), one can use the following formula $t=-t_{0}^{\alpha}\ln\left(u\right)\left(\frac{\sin\left(\gamma^{\alpha}\pi\right)}{\tan\left(\gamma^{\alpha}\pi v\right)}-\cos\left(\gamma^{\alpha}\pi\right)\right)^{\frac{1}{\gamma^{\alpha}}}$ Fulger et al. (2008); Kozubowski (2000), where $u$ and $v$ are independent random numbers from a uniform distribution on $(0,1]$.

**Stochastic fluctuations and quasi-pattern formation in reaction-diffusion systems with anomalous transport

 

Supplemental Material

**

One main result of our work is a Gaussian approximation of the fluctuations in a reaction-diffusion system with anomalous transport (i.e. generalised waiting-time distributions between hopping events). To derive this approximation, we find the generating functional for the reaction-diffusion system, and then expand this functional to sub-leading order in $N^{-1}$, where the parameter $N$ characterises the typical number of particles per site in the discrete lattice. This approximate path-integral is recognised as the generating functional of a stochastic different equation (SDE), with a special form, to be specified below. One can then compare the two generating functionals and read off the effective SDE for the reaction-diffusion system. This supplement sets out the different steps of the calculation, and subsequent analysis.

Appendix S1 Generating functional for an SDE with age variables

We first present a method of calculating the generating functional for a stochastic process with Gaussian noise. The choice of the form of the equations for which we find the corresponding generating functional is important and is informed by the type of equations we expect for the reaction-subdiffusion systems in the main paper.

To start with, we discretise time using step size $\Delta$. We will eventually take the continuous limit in time. The set of stochastic equations has the discretised form

[TABLE]

where $f^{\alpha}_{i,\tau,t}\left(\{x^{\alpha}_{i,\tau,t}\}\right)$ is an arbitrary function of the coordinates $\{x^{\alpha}_{i,\tau,t}\}$ that informs the deterministic trajectory of the system. The $\{\chi^{\alpha}_{i,\tau,t}\}$ are Gaussian random variables, which encode the fluctuations about the deterministic trajectory. We write their correlations as

[TABLE]

The matrix $C^{\alpha,\alpha^{\prime}}_{i,i^{\prime},\tau,\tau^{\prime},t,t^{\prime}}$ is symmetric with respect to the pairs of indices $(\alpha,\alpha^{\prime}$), $(i,i^{\prime})$, $(t,t^{\prime})$ and $(\tau,\tau^{\prime})$.

In continuous time ($\Delta\to 0$) the expressions in Eqs. (S2,S3) correspond to

[TABLE]

where $\xi^{\alpha}_{i,\tau,t}=\lim_{\Delta\to 0}\frac{\chi^{\alpha}_{i,\tau,t}}{\sqrt{\Delta}}$. We note that the first of these is a stochastic partial differential equation. The second relation is technically not a differential equation (it contains no derivatives with respect to time). We will nevertheless refer to the combination of Eqs. (S5) and (S6) as a stochastic differential equation (SDE).

We begin by writing down the generating functional

[TABLE]

for the process in Eqs. (S2,S3). The quantity $\mathcal{P}\left(\{x^{\alpha}_{i,\tau,t}\}\right)$ is the probability measure in the space of trajectories for the system. Assuming that the correspondence of the set of noise variables $\{\chi^{\alpha}_{i,\tau,t}\}$ to a particular trajectory $\{x^{\alpha}_{i,\tau,t}\}$ is one-to-one, we can write $\mathcal{P}\left(\{x^{\alpha}_{i,\tau,t}\}\right)=\langle\delta\left(\boldsymbol{x}-\boldsymbol{x}\left(\boldsymbol{\chi}\right)\right)\rangle_{\{\chi^{\alpha}_{i,\tau,t}\}}$, where bold face notation indicates entire trajectories.

Then letting $X^{\alpha}_{i,\tau+\Delta,t+\Delta}\equiv x^{\alpha}_{i,\tau+\Delta,t+\Delta}-x^{\alpha}_{i,\tau,t}-\frac{\Delta}{N}f^{\alpha}_{i,\tau,t}\left(\{x^{\alpha}_{i,\tau,t}\}\right)-\frac{\sqrt{\Delta}}{N}\chi^{\alpha}_{i,\tau,t}$ and $X^{\alpha}_{i,0,t+\Delta}\equiv x^{\alpha}_{i,0,t+\Delta}-\frac{1}{N}f^{\alpha}_{i,0,t}\left(\{x^{\alpha}_{i,\tau,t}\}\right)-\frac{\chi^{\alpha}_{i,0,t}}{N\sqrt{\Delta}}$ we have

[TABLE]

where the Jacobian determinant $\lvert\frac{\delta\boldsymbol{X}}{\delta\boldsymbol{x}}\rvert=1$ because the Jacobian itself is a triangular matrix with diagonal elements each equal to unity Altland and Simons (2010).

The generating functional thus takes the following form

[TABLE]

where the $\delta$-functions enforce the stochastic dynamics in Eqs. (S2,S3). For the sake of brevity, we have adopted the Einstein summation convention for sums over repeated indices in the exponential, and have omitted the summation symbols $\sum$. The quantity $P\left(\{\chi^{\alpha}_{i,\tau,t}\}\right)$ is the distribution of the noise variables $\{\chi^{\alpha}_{i,\tau,t}\}$; we keep in mind that these may depend on the variables $\{x^{\alpha}_{i,\tau,t}\}$.

Introducing conjugate variables $\{\hat{x}^{\alpha}_{i,\tau,t}\}$ to rewrite the $\delta$-functions in their exponential representation one obtains

[TABLE]

We can now carry out the Gaussian integrals over the $\{\chi^{\alpha}_{i,\tau,t}\}$, using their joint distribution

[TABLE]

where $\Omega$ is a normalisation constant. We define the inverse of $C^{\alpha,\alpha^{\prime}}_{i,i^{\prime},\tau,\tau^{\prime},t,t^{\prime}}$ such that

[TABLE]

The integral over $\{\chi^{\alpha}_{i,\tau,t}\}$ can be written as follows:

[TABLE]

We now complete the square in the exponent and carry out the integral over the noise variables. To do this let $\chi^{\alpha}_{i,\tau,t}=-i\frac{\sqrt{\Delta}}{N}C^{\alpha,\alpha^{\prime}}_{i,i^{\prime},\tau,\tau^{\prime},t,t^{\prime}}\hat{x}^{\alpha^{\prime}}_{i^{\prime},\tau^{\prime},t^{\prime}}+y^{\alpha}_{i,\tau,t}$. Then,

[TABLE]

Thus, if we change the set of integration variables to \Big{\{}y^{\alpha}_{i,\tau,t}\Big{\}}, the integral becomes

[TABLE]

Therefore, we can rewrite the generating functional as

[TABLE]

In Eq. (S16), sums over repeated indices are implied. Note that the sum over the residence time variables $\tau$ are not to include $\tau=0$; this special case is dealt with separately.

Thus, we have recast the generating functional in terms of only the coordinates and their conjugates, eliminating the noise variables in favour of their correlations. Now that we have the general form for the generating functional for a process of the form in Eq. (S2,S3), we can compute the generating functional for the reaction-subdiffusion system and compare the results. This allows us to write down an effective SDE for the reaction-subdiffusion system.

It should be noted that terms in the exponent in Eq. (S16) which are quadratic in the conjugate variables are associated with the noise, and the coefficients are the correlations between the noise variables. The terms which are linear in the conjugate variables arise due to the deterministic part of the dynamics.

Appendix S2 Generating functional for a reaction-diffusion system with non-constant hazard rate

S2.1 Construction of the generating functional

In order to obtain effective SDEs for the type of (anomalous) reaction-diffusion systems in the main paper, we seek to write down the generating functional for the process and reduce it to a form similar to that in Eq. (S16). We begin by reminding ourselves of Eqs. (8) and (9) from the main paper; recalling the notation $x_{i,\tau,t}^{\alpha}=n_{i,\tau,t}^{\alpha}/N$ these can be written as

[TABLE]

[TABLE]

where $k^{\left(R\right)\alpha}_{i,r,\tau,t}$ is number of particles of type $\alpha$ and age $\tau$ at position $i$ which are annihilated in the time step from $t$ to $t+\Delta$ due to reactions of type $r$. Similarly, $k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}$ is the number of particles of age $\tau$ hopping away from position $i^{\prime}$ to $i$ at time $t$. Finally, $\ell_{i,r,t}^{\left(R\right)}$ is the number of reactions of type $r$ firing in the time window $t$ to $t+\Delta$ at position $i$. It should be noted that $k^{\left(R\right)\alpha}_{i,r,\tau,t}$, $k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}$ and $\ell_{i,r,t}^{\left(R\right)}$ are always non-negative integers and $\theta\left(-\nu_{r}^{\alpha}\right)\ell^{\left(R\right)}_{i,r,t}\lvert\nu_{r}^{\alpha}\rvert=\sum_{\tau}k^{\left(R\right)\alpha}_{i,r,\tau,t}$.

Following a procedure similar to that in Section S1, we can write Eq. (S7) in the following form

[TABLE]

Here, we use $\{k\}$ as a shorthand for the set of all the variables $\left\{\ell^{\left(R\right)}_{i,r,t}\right\}$ ,$\left\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\right\}$ and $\left\{k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}\right\}$ and $\{x\}$ as shorthand for the set $\{x_{i,\tau,t}^{\alpha}\}$. The statistics of these variables is informed by the reaction rates, which in turn depend on the $\{x\}$. The joint probability distribution for the set $\{k\}$, given the set $\{x\}$, is denoted by $\mathcal{P}\left(\{k\}|\{x\}\right)$. The precise form of this distribution is discussed below.

Similarly to the previous section, we write the delta function as the integral of a complex exponential,

[TABLE]

Our aim is to carry out the sum over $\{k\}$ in order to obtain an expression for the generating functional in terms of only the coordinates $\{x^{\alpha}_{i,\tau,t}\}$ and the model parameters. Using this, we can then perform an expansion in $N^{-1}$ to obtain an approximate expression which has the same form as Eq. (S16). This will allow us to read off the set of effective SDEs. It should be noted that the values $\{k\}$ are of order $N^{0}$, as are the $\{x^{\alpha}_{i,\tau,t}\}$.

We must now identify the precise nature of $\mathcal{P}\left(\{k\}|\{x\}\right)$ in order to evaluate the summation over $\{k\}$. First, due to the way that we have deconvolved the problem by introducing the age-variable $\tau$, we observe that the values of $\{k\}$ at different times $t$ are independent. Further, the hopping events are independent of the reactions. However, the values of $k^{\left(R\right)\alpha}_{i,r,\tau,t}$ at a particular time $t$ are dependent upon the number of firings of the different reaction types $\ell^{\left(R\right)}_{i,r,t}$. Introducing the notation $\{\cdot\}_{t}$ to mean the set of all variables for a fixed value of $t$, we write

[TABLE]

We can simplify matters further by noting that the hops of particles of different flavours $\alpha$, ages $\tau$ and starting/ending positions $i$ and $i^{\prime}$ also occur independently. So, we can write

[TABLE]

Similarly, reactions of different types $r$ and at different sites $i$ occur independently of one another. We also presume that reaction rates do not depend on the ages of the particles, only the total local concentration. Thus, we can write

[TABLE]

Finally, the numbers of particles of various ages $\tau$ which are created/annihilated in reactions of type $r$ at location $i$ and time $t$ obey the following relation

[TABLE]

We note again that the condition $\sum_{\tau}k^{\left(R\right)\alpha}_{i,r,\tau,t}=\ell^{\left(R\right)}_{i,r,t}\lvert\nu_{r}^{\alpha}\rvert\theta\left(-\nu_{r}^{\alpha}\right)$ must be met, so the values of $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}_{\alpha,i,r,t}$ are interdependent. Using Eqs. (LABEL:sep1-S23), we can write Eq. (S19) in the form

[TABLE]

where we have gathered the reaction terms

[TABLE]

The hopping terms are given by

[TABLE]

Now we evaluate the sums over $\left\{\ell^{\left(R\right)}_{i,r,t}\right\}$ ,$\left\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\right\}$ and $\left\{k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}\right\}$ in Eqs. (S25) and (S26). In order to do this, we first must specify the distributions $P\left(k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}|x_{i,\tau,t}^{\alpha}\right)$, $P\left(\ell^{\left(R\right)}_{i,r,t}|\{x_{i,t}^{\alpha}\}_{i,t}\right)$ and $P\left(\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}_{\alpha,i,r}|\ell^{\left(R\right)}_{i,r,t},\{x_{i,t}^{\alpha}\}_{i,t}\right)$. After these sums are evaluated, we perform the system-size expansion in powers of $N^{-1}$. The calculation that follows is based on the assumption made for example in the context of Gillespie’s $\tau$-leaping algorithm Gillespie (2001). Reaction rates are assumed to be constant during each time interval $\Delta$. Hence the number of reactions of type $r$ firing in such an internal is a Poissonian random number with parameter $\lambda_{i,r,t}\times\Delta$. The statistics of the number of hopping events are Poissonian as well, with the appropriate rates.

First, we deal with the hopping contribution. The distribution for $k^{\left(H\right)\alpha}_{i,i^{\prime},\tau,t}$ is Poissonian with mean $m^{\alpha}_{i,i^{\prime},\tau,t}=\Delta^{2}\phi_{i,i^{\prime}}h^{\alpha}_{\tau}n_{i,\tau,t}^{\alpha}$ where $h^{\alpha}_{\tau}$ is the rate of hopping for particles which have survived a time $\tau$. That is, we have

[TABLE]

Substituting Eq. (S27) into Eq. (S26) and evaluating the sum over $k_{i,i^{\prime},\tau,t}^{\left(H\right)\alpha}$ , one then obtains

[TABLE]

In the last step we have carried out the expansion in $N^{-1}$, and have retained terms up to and including sub-leading order.

The procedure for the reaction events is more complicated, but follows similar lines. The distribution for $k^{\left(R\right)\alpha}_{i,r,\tau,t}$ is dependent on the number of firings of the different reaction types, $\ell^{\left(R\right)}_{i,r,t}$. The distribution for $\ell^{\left(R\right)}_{i,r,t}$ is Poissonian. The mean number of firings of a reaction of type $r$ in time $\Delta$ is $\Delta\lambda_{i,r,t}$. The rate $\lambda_{i,r,t}$ is dependent on the overall concentrations $\{x_{i,t}^{\alpha}\}$. The distribution of the number of firings of type $r$ at position $i$ and time $t$ is therefore

[TABLE]

The conditional joint distribution for $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}_{\alpha,i,r,t}$ is multinomial if particles of type $\alpha$ are annihilated in the reaction of type $r$. If particles of type $\alpha$ are not annihilated in reactions of type $r$, $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}_{\alpha,i,r,t}$ must all be zero. The distribution is multinomial because the ages $\tau$ of the particles to be annihilated are selected randomly, with each value of $\tau$ coming up with probability proportional to $\frac{x_{i,\tau,t}^{\alpha}}{x_{i,t}^{\alpha}}$. The values of $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}_{\alpha,i,r,t}$ must also be such that $\sum_{\tau}k^{\left(R\right)\alpha}_{i,r,\tau,t}=\ell^{\left(R\right)}_{i,r,t}\lvert\nu_{r}^{\alpha}\rvert\theta\left(-\nu_{r}^{\alpha}\right)$. So the conditional joint distribution for $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}_{\alpha,i,r,t}$ is

[TABLE]

where we use $\delta^{(K)}\left(\cdot\right)$ to denote a Kronecker delta function which is equal to one when the argument is zero and is equal to zero otherwise. We note that the expression $\left(\ell^{\left(R\right)}_{i,r,t}\lvert\nu^{\alpha}_{r}\rvert\right)!\left(\prod_{\tau}\frac{1}{k^{\left(R\right)\alpha}_{i,r,\tau,t}!}\left(\frac{x_{i,\tau,t}^{\alpha}}{x_{i,t}^{\alpha}}\right)^{k^{\left(R\right)\alpha}_{i,r,\tau,t}}\right)$ in Eq. (S30) reduces to the familiar binomial distribution if the number of values that $\tau$ can take is limited to two. As a sanity check we use the multinomial expansion to verify normalisation,

[TABLE]

Using equation (S30), one can evaluate the sum over $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}_{i,r,t}$ in Eq. (S25) to show that

[TABLE]

where again we have used the multinomial expansion to evaluate the sums over $\{k^{\left(R\right)\alpha}_{i,r,\tau,t}\}$. Combining this result with Eq. (S25), one obtains

[TABLE]

If one then uses the expression in Eq. (S29), one can perform a similar deduction to that used to find Eq. (S28) to finally obtain

[TABLE]

This can then be expanded in powers of $\frac{1}{N}$. To the lowest two orders one obtains

[TABLE]

In the above, ‘${\rm pairs}(\alpha,\alpha^{\prime})$’ means each combination of $\alpha$ and $\alpha^{\prime}$ with $\alpha\neq\alpha^{\prime}$ such that $(\alpha,\alpha^{\prime})$ and $(\alpha^{\prime},\alpha)$ are the same ‘pair’. Neglected terms are of order $N^{-3}$. Having found approximate expressions for $R_{i,r,t}$ and $H^{\alpha}_{i,i^{\prime},\tau,t}$ in Eq. (S28) and Eq. (S35), we can write down the final expression for the generating functional

[TABLE]

We have used the fact that $\frac{n^{\alpha}_{i,\tau,t}}{n^{\alpha}_{i,t}}=\frac{x^{\alpha}_{i,\tau,t}}{x^{\alpha}_{i,t}}$. Eq. (S36) contains terms in the exponentials which are at most quadratic in the conjugate variables $\{\hat{x}\}$. By comparing Eq. (S36) with Eq. (S16) we can read off expressions for the effective SDEs for the individual-based system in the Gaussian approximation.

The terms which are linear in the conjugate variables (and with pre-factors of order $N^{-1}$) represent the dynamics in the deterministic limit. The quadratic terms in Eq. (S36) carry pre-factors $N^{-2}$; the coefficients in these terms represent the correlations between the noise variables. These correlators are thus given by

[TABLE]

S2.2 Results of the generating functional approach

Once continuous time is restored by sending $\Delta\to 0$, the effective SDE found from the comparison of Eq. (S36) with Eq. (S16) is given by

[TABLE]

These are the expressions in Eq. (11) in the main paper. When one takes the continuous-time limit one uses the fact that $\frac{\delta_{tt^{\prime}}}{\Delta}\to\delta\left(t-t^{\prime}\right)$ as $\Delta\to 0$, and that in going from Eqs. (S2) and (S3) to Eq. (S38) the Gaussian variables are rescaled in the following manner $\frac{\chi_{i,\tau,t}^{\alpha}}{\sqrt{\Delta}}\to\xi_{i,\tau,t}^{\alpha}$. One finds

[TABLE]

In the Markovian case the noise variables of concern are $\eta^{\alpha}_{i,t}=\int_{0}^{t}\xi^{\alpha}_{i,\tau,t}d\tau+\xi^{\alpha}_{i,0,t}$ (as stated in Eq. (16) of the main paper) and the hazard rate $h^{\alpha}_{\tau}=h^{\alpha}$ is a constant with respect to $\tau$. As a check, we calculate the correlators of these noise variables and find

[TABLE]

This is the same result one would obtain by using a standard master equation approach and the Kramers-Moyal expansion or by the direct application of Kurtz’ theorem Brett and Galla (2013); kurtz; Gardiner (2009).

To summarise, we have found the approximate generating functional for a non-Markovian reaction-diffusion system. This was done by identifying the probability distributions for the numbers of the various events that might occur, carrying out the summations over $\{k\}$ and, finally, by using a system-size expansion to obtain the same form as in Eq. (S16). This enabled us to read off the effective SDEs for the system, given in Eq. (11) of the main paper, and the correlation coefficients for the stochastic variables by equating terms of the same order in $\frac{1}{N}$.

Appendix S3 Derivation of the anomalous reaction-diffusion equation with noise

Beginning with Eq. (S38), we wish to derive the general reaction-diffusion equation

[TABLE]

This is Eq. (13) in the main paper. This is done using a similar method to those used in Vlad and Ross (2002), Fedotov and Falconer (2012) or yadav. Here, however, we also handle the noise variables. Using the method of characteristics, one can show from Eq. (S38) that

[TABLE]

where we have used $\Psi^{\alpha}\left(\tau\right)=e^{-\int_{0}^{\tau}h^{\alpha}_{s}ds}$. Also, one can integrate the first of Eqs. (S38) with respect to $\tau$ in order to obtain

[TABLE]

All that remains is to eliminate the integral over $\tau$ in Eq. (S43) in favour of an expression involving only $n^{\alpha}_{i,t}$. We do this using the convolution theorem for Laplace transforms and Eq. (S42). Integrating Eq. (S42) with respect to $\tau$ and then taking the Laplace transform we obtain

[TABLE]

One the other hand we can also multiply both sides of Eq. (S42) by $h_{\tau}^{\alpha}$, then integrate over $\tau$, and then take the Laplace transform. We then arrive at

[TABLE]

One can use the combination of Eqs. (S44) and (S45) to eliminate $\mathcal{L}_{t}\left[n^{\alpha}_{i,0,t}e^{\int_{0}^{t}p^{\alpha}_{i,s}ds}\right]$, and to find an expression for $\int_{0}^{t}h^{\alpha}_{\tau}n^{\alpha}_{i,\tau,t}d\tau$ in terms of $n^{\alpha}_{i,t}$. Defining the memory kernel $K^{\alpha}_{t}=\mathcal{L}^{-1}_{u}\left(t\right)\left\{\frac{\psi^{\alpha}\left(u\right)}{\Psi^{\alpha}\left(u\right)}\right\}$ and inverting the Laplace transforms, this finally leads one to

[TABLE]

with the noise term

[TABLE]

If we now choose the waiting-time distribution to be the Mittag-Leffler function,

[TABLE]

then we obtain Eq. (S41) since in this case $\int_{0}^{t}K^{\alpha}_{t-T}f\left(T\right)dT=\left(t_{0}^{\alpha}\right)^{-\gamma^{\alpha}}{}_{0}D_{t}^{1-\gamma^{\alpha}}\left\{f\left(t\right)\right\}$.

Appendix S4 Linear-noise approximation

We define the deviation from the homogeneous fixed point as $\delta^{\alpha}_{i,\tau,t}=n^{\alpha}_{i,\tau,t}-\bar{n}^{\alpha}$. We presume that $\delta^{\alpha}_{i,\tau,t}\sim\mathcal{O}\left(\sqrt{N}\right)$. One can then expand Eq. (S46) as a series in $\frac{1}{\sqrt{N}}$. We discard terms $\mathcal{O}\left(N^{0}\right)$. The leading terms are of order $N$, and correspond to the deterministic trajectory. The sub-leading terms are $\mathcal{O}\left(\sqrt{N}\right)$, and encapsulate the behaviour of the fluctuations to a linear approximation. Carrying out this expansion, we have

[TABLE]

where the concentration-dependent terms in the expression for the noise are now evaluated at the homogeneous fixed point such that

[TABLE]

The correlators of $\{\xi^{\alpha}_{i,\tau,t}\}$ are also evaluated at the fixed point. In order to simplify the above expression for the noise, one notes that

[TABLE]

Using the Laplace transform and the corresponding convolution theorem, this enables us to define the following simplified noise variables

[TABLE]

The expression in Eq. (S52) is preferable that in Eq. (S50) because it involves integrals over at most two variables, as opposed to three. This simplification makes the calculation of the correlators of the noise variables much more straightforward. Taking the Laplace transform of Eq. (S49) and multiplying through by $\hat{\Psi}^{\alpha}\left(u+\bar{p}^{\alpha}\right)$, one obtains

[TABLE]

This can then be inverted to find an expression for the fluctuations in terms of the noise variables. Observing that Eq. (S53) represents a matrix equation of the form $\underline{\underline{\hat{\tilde{m}}}}^{-1}_{q,u}\underline{\hat{\tilde{\delta}}}_{q,u}=\underline{\hat{\tilde{\epsilon}}}_{q,u}$, we obtain the solution

[TABLE]

The elements of the matrix $\underline{\underline{\tilde{m}}}_{q,t}$ are the Green functions of Eq. (S49) without the noise term.

Appendix S5 Correlators for the noise variables and fluctuations

The correlators of the noise variables in Eq. (S52) can be derived using Eqs. (S39) and are given by

[TABLE]

and

[TABLE]

We have defined the function $\chi^{\left(\alpha\right)}\left(t\right)=\int_{0}^{\infty}\Psi_{\star}^{\left(\alpha\right)}\left(t+T\right)dT$ and used the shorthand $\Psi_{\star}^{\left(\alpha\right)}\left(T\right)=\Psi^{\left(\alpha\right)}\left(T\right)e^{-\bar{p}^{\left(\alpha\right)}T}$. We note also that $\hat{\chi}^{\left(\alpha\right)}\left(u\right)=\frac{1}{u}\left(\hat{\Psi}\left(\bar{p}^{\left(\alpha\right)}\right)-\hat{\Psi}\left(u+\bar{p}^{\left(\alpha\right)}\right)\right)$. The remaining coefficients in Eqs. (S55) and (S56) are given by

[TABLE]

We have used the long-time limit to approximate the barred quantities as constant in time. Note that

[TABLE]

The following identities for the double Laplace transform allow us to simplify our expressions considerably

[TABLE]

Using Eq. (S54) one can then use the expressions for the noise correlators to calculate the correlators and power spectra of the fluctuations $\delta^{\alpha}_{i,t}$. These are given by

[TABLE]

[TABLE]

In order to evaluate the integrals in Eq. (S61), the inverse Laplace transforms of the functions in the integrands are taken numerically using the Zakian method Halsted and Brown (1972) and the integrals are performed numerically.

We note that in the long-time limit the correlators for both the noise variables and the fluctuations can be shown to be time-translation invariant, i.e., they are a function of only $t-t^{\prime}$. This is due to the fact that, since the functions in the integrands decay as their arguments are increased, the integrals can be approximated in the following way in the long-time limit

[TABLE]

The evaluation of the equal-time correlator can be performed without the use of the numerical inverse Laplace transform. The equal-time correlator involves integrals of the form

[TABLE]

In the long-term, one can use the following theorem to evaluate the equal-time correlator

[TABLE]

where $c$ is a real number to the right of all the poles of either function in the integrand and we have used the dual convolution theorem for Laplace transforms debnath. This result is an analogue of Parseval’s theorem, which applies to Fourier transforms. The result Eq. (S64) allows one to calculate the equal-time power-spectra of the fluctuations without using a numerical Laplace transform. Instead, one merely has to evaluate a contour integral numerically. Let $A^{\beta\beta^{\prime}}_{q}\left(u,u^{\prime}\right)=\left(u+u^{\prime}\right)\langle\hat{\tilde{\epsilon}}^{\left(\beta\right)}_{q,u}\hat{\tilde{\epsilon}}^{\left(\beta^{\prime}\right)\star}_{q,u^{\prime}}\rangle$. The equal-time power spectrum is then given by

[TABLE]

Appendix S6 Accuracy of the system size expansion and low-$N$ behaviour

The system-size expansion performed in Section S2 assumes a large value of $N$. In this section, we briefly discuss how large the system-size is required to be for our theory to be accurate and whether the conclusions drawn in the main paper can be extended to lower system sizes.

Fig. S1 demonstrates the convergence to the analytical result for the power spectrum as the system size is increased, in the case of the Brusselator model. One notices that the inaccuracy is greatest for lower $q$ and that, as a result, the stochastic patterns are more prominent in the system for lower $N$. That is, the peak of the power spectrum surpasses the value at $q=0$ to a greater extent for lower $N$ than for higher $N$.

For reference, we also include examples of how the deterministicly-driven and noise-driven patterns manifest in a system which has fewer particles. Fig. S2 demonstrates that noise plays a more prominent role for smaller $N$, resulting in the deterministicly-driven patterns being far less orderly and uniform than for the large-$N$ case presented in the main text. The noise-driven patterns however are slightly less affected than the deterministicly-driven patterns. This could be attributed to the fact that the peak of the power spectrum is slightly more prominent for low $N$, as discussed above.

In Fig. S3, we show the power spectra for the Lengyel-Epstein model in four cases: $N=4000$, $\gamma=0.5$; $N=100$, $\gamma=0.5$; $N=4000$, $\gamma=0.7$; $N=100$, $\gamma=0.7$. All other system parameters, including the ratio of effective diffusion coefficients $\theta$, are held constant. We find that the peak in the power spectrum rises (relative to the value at $q=0$) as $N$ is reduced (regardless of the value of $\gamma$) and as $\gamma$ is reduced (regardless of the value of $N$). This is a common feature for other sets of parameters. This verifies that the threshold value of $\theta$ required for stochastic pattern formation still reduces as the activator is made more subdiffusive for smaller values of $N$, as is the case in Fig. 4(a) of the main text. Furthermore, the critical value of $\theta$ actually reduces as $N$ is made smaller, meaning that the region of parameter space for which stochastic patterns occur actually appears to be larger for smaller $N$.