Effective diffusion coefficients in reaction-diffusion systems with anomalous transport

We show that the Turing patterns in reaction systems with subdiffusion can be replicated in an effective system with Markovian cross-diffusion. The effective system has the same Turing instability as the original system, and the same patterns. If particles are short-lived, the transient dynamics are captured as well. We use the cross-diffusive system to define effective diffusion coefficients for the system with anomalous transport, and we show how they can be used to efficiently describe the Turing instability. We also demonstrate that the mean squared displacement of a suitably defined ensemble of subdiffusing particles grows linearly with time, with a diffusion coefficient which agrees with our earlier calculations. We verify these deductions by numerically integrating both the fractional reaction-diffusion equation and its normally diffusing counterpart. Our findings suggest that cross-diffusive behaviour can come about as a result of anomalous transport.


I. INTRODUCTION
In his seminal paper [1] Turing provided a simple mechanism for spontaneous pattern formation in chemical reaction systems. Since this work, the reaction-diffusion paradigm [2] has been used ubiquitously as a tool for the description of pattern-forming phenomena, ranging from predator-prey interactions [3] to developmental biology [4,5].
In this context, one typically assumes that constituents diffuse in a manner such that the mean squared displacement of individuals grows linearly with time. Mathematically, r 2 ∼ t, where r denotes the position vector of a typical particle, and where the angular brackets define an average over trajectories. The rate of diffusion can then be quantified by means of the diffusion constant D = r 2 /(2t) [6,7]. In one spatial dimension this simplifies to D = x 2 /(2t). However, the situation may not always be this simple. Due to the many complicated interactions between the reactants and their environment, the diffusion may be of an 'anomalous' nature. For example, the particles or individuals involved may be subject to 'trapping' effects, which can result in subdiffusive behaviour: x 2 ∼ t α where 0 < α < 1 [8][9][10]. Reactiondiffusion equations with anomalous (or 'non-Fickian') diffusion have been used to describe systems ranging from the neolithic transition in Europe [11,12] to photoluminescence in semi-conductors [13] to morphogen gradient formation [14,15].
How, then, might one go about defining a quantity which describes the speed of anomalous diffusion? What is the analogue of the diffusion constant D? In order to be useful, such a quantity must not only inform us as to how quickly different reactant types diffuse, but it must also be of physical relevance and have predictive power. For example, a criterion for pattern formation in normally diffusive reaction-diffusion systems can be written in terms of the diffusion constants and the other system parameters [16]. So, by introducing an 'effective diffusion coefficient', we might aim to produce a similar criterion for pattern formation in systems with anomalous diffusion.
The notion of an effective diffusion coefficient has been utilised previously in the study of transport in confined geometries [17] and for tilted periodic potentials [18,19]. The possibility of defining an effective diffusion constant in a reaction-subdiffusion system was mentioned in [20]. In this paper we present a more detailed argument for the definition of effective diffusion coefficients, and seek to interpret the results.
We approach the problem as follows. We start from a reaction-diffusion system in which the reactants can exhibit subdiffusive behaviour. We refer to this as the 'original' system. We then derive an 'effective' normally diffusing reaction-diffusion system, which exhibits pattern formation for the same sets of model parameters as the original anomalously diffusing system. Furthermore, we argue that the normally diffusing system gives rise to the same patterns as the original system in the long-term and, under certain conditions, replicates the transient dynamics of the original system to a good approximation. It is the diffusion coefficients in this normally diffusing system which we use to characterise the pattern formation in the original subdiffusive system. Although we use subdiffusion as the primary example in this work, our method could be applied to other varieties of non-Fickian diffusion than just the subdiffusive kind.
We find that the effective diffusion in reaction-diffusion systems with anomalous transport is in general crossdiffusive. That is, the diffusion of one substance may be affected by a concentration gradient in another. As such, we suggest that anomalous transport can be a possible mechanism for the production of cross-diffusive effects in experiments, in addition to those that are already known.
We also show that the collective behaviour of constituents in the stationary state can give rise to a mean squared displacement which scales linearly with time, x 2 ∼ t, despite the underlying subdiffusive motion of the individuals. To demonstrate this, we consider an en-semble in which particles are subject to diffusion and removal events but are replaced upon removal, thus keeping the ensemble size the same, as would be the case in the homogeneous steady state. The effective diffusion coefficient found from this reasoning agrees with our earlier definitions.
The remainder of this paper is set out as follows: In Section II, we present the necessary background information. In particular, we summarise the main features of reaction-subdiffusion systems, which we use in the rest of the work. In Section III, we construct the corresponding 'effective' normally diffusing system. We go on to show that the condition for Turing pattern formation in the subdiffusive system can be written in terms of the effective diffusion coefficients and that the normally diffusive system produces the same patterns as the original system. Additionally, we show that, under certain circumstances, the normally diffusing system approximates the time-dependent behaviour of the original system. The interpretation of the effective diffusion coefficients in terms of the underlying microscopic processes is also discussed. In Section IV, we confirm the theoretical results by numerical integration of the reaction-subdiffusion equation, and of its normally diffusing counterpart. Finally, we summarise and discuss our findings in Section V.

II. MODEL DEFINITION AND BACKGROUND
In this section, we discuss the formulation of reactiondiffusion systems in terms of a continuous-time random walk (CTRW) model. The CTRW approach is a standard method by which to derive reaction-diffusion equations with anomalous diffusive behaviour [8][9][10]. We also discuss the stability of the homogeneous steady states in such systems.

A. General reaction-diffusion model
We consider a system in one spatial dimension which involves different types of particles. These particles may hop from location to location and be involved in reactions. When a reaction occurs, particles of each flavour can be created or destroyed. In principle this defines a stochastic system if the overall particle number is finite. In this paper, however, we ignore the stochastic fluctuations and examine only the deterministic behaviour. In doing so, we are assuming that the total number of particles is sufficiently large, so that the fluctuations are comparatively negligible. A further convenient consequence of dealing with large particle numbers is that the discrete quantity of particles can be approximated as a continuous concentration.
Within this approximation, the local particle number density (or concentration) of a particular species i at a time t and at location x is denoted by ρ i (x, t). This means that x2 x1 ρ i (x, t) dx is the number of particles of type i in the region x 1 < x < x 2 . We write ρ for the vector (ρ 1 , ρ 2 , . . . ), i.e., ρ(x, t) indicates the concentrations of all particles types at a given location and time. Given the initial location of a particle of type i, a waiting time is drawn from a distribution ψ i (t). Waiting times for each particle are drawn independently. We do not specify these distributions at this point for the purposes of generality. Each particle then remains where it is until it has waited the allotted time, at which point it hops by a distance x, drawn from a hopping kernel φ i (x). The spatial coordinate x can be continuous or can denote a position on a discrete lattice. Once the particle has hopped, a new waiting time is drawn and the process repeats itself.
Given these definitions, one can define the hazard rate of hopping dτ is the probability that a particle has not hopped for a time t since its last hop. The hazard rate characterises the proclivity of particles to hop a time t after their last hop. It is constant in the case of regular diffusion but time-dependent in the general case.
Particles may also be involved in reactions during their sojourn periods at a fixed location. We label these different types of reactions with an index r. They occur with rates which depend on the concentration of particles at that location. The result of a reaction is the local production and/or annihilation of particles. We denote the rate at which a reaction of type r occurs by λ r [ρ (x, t)] and the number of particles of type i which are produced or annihilated in such a reaction by ν r i . If particles are annihilated, scheduled hopping events involving these particles no longer occur. Particles which are involved in reactions have their waiting-times redrawn; this approach applies to irreversible reactions, but special consideration has to be given in the case where the reactions are reversible [21]. We do not consider such reaction schemes in this paper.
In order to deal with the non-constant hazard rates h i (t) and to derive an equation which describes the time evolution of the concentrations ρ(x, t), one can introduce an additional temporal coordinate τ -the time since the last hop for an individual particle. The concentrations can then be subdivided by introducing the quantities ρ(x, τ, t) where ρ(x, t) = t 0 ρ(x, τ, t)dτ . Using such a coordinate, the problem is recast as Markovian. One is then able to derive a generalised reaction-diffusion equation. A detailed derivation of equations of this type can be found in [11,22]. Carrying out a calculation along these lines, one obtains We have defined the following objects: where H(u) = 1 if u ≥ 0 and H(u) = 0 otherwise. The quantities R ± i are the net production and annihilation rates, respectively, of particles of type i. In Eq. (1) K i (τ ) is the 'memory kernel' defined through its Laplace transform, The Laplace transform is given by (1), the reaction term, R + (ρ) − R − (ρ), and the diffusion term [which contains convolutions with both the hopping kernel φ (x) and the memory kernel K i (t)] are coupled via the exponential factor in the convolution with the memory kernel. This coupling arises from the fact that particles may be annihilated before they hop: If the death rate R − i were identically zero, the contribution of the concentration at a time t − τ to the integral over time in Eq. (1) would be K i (τ ) ρ i (x , t − τ ). To account for a non-zero death rate, this contribution is weighted by the probability that a particle of type i survives the time interval from t − τ to t; this probability is given by e . This coupling term is absent in the case of Markovian hopping. That is, if we choose the exponential waiting time distribution ψ i (t) = 1 τi e − t τ i we obtain K i (t) = 1 τi δ (t), and the reaction and diffusion terms in Eq. (1) decouple. In general though, the timeevolution of the concentrations depends on the history of the system, and the process is manifestly non-Markovian.

B. Reaction-subdiffusion system
In this paper, we will primarily concern ourselves with a specific case of non-Markovian transport: subdiffusion. We choose the waiting-time distribution to be the Mittag-Leffler distribution [23] where 0 < α i < 1 and where E αi (·) is the oneparameter Mittag-Leffler function [24]. This function can be thought of as a generalised exponential function and is defined as where Γ (·) denotes the standard gamma function. The constant τ i has units of time; its role is to fix the overall time scale of the hopping process. It is important to note that the Mittag-Leffler distribution has the 'heavy-tail' That is, in absence of annihilation the first moment of the waiting times between hops diverges. This makes the Mittag-Leffler function suitable for modelling the hopping of particles subject to 'trapping' effects [8].
We allow the hopping kernel to take a general symmetrical form. If one presumes that the hopping kernel is of the form φ i (x) = Φ i (x/ξ), where ξ characterises the typical hopping distance and Φ i (·) is even, the Fourier transform of the hopping kernel has following asymptotics for small ξ, where σ i is a constant related to the variance of the hopping distances such that − ∂ 2 ∂k 2φ (k) | k=0 = 2σ 2 i ξ 2 . If one were to track the movement of an ensemble of particles which hopped with waiting-times drawn from the distribution in Eq. (4) and hopping distances drawn from a hopping kernel described by Eq. (6), one would obtain subdiffusive behaviour, characterised by the mean squared displacement x 2 ∼ t αi [8]. This can be shown using the Montroll-Weiss formula [25]. The parameter α i characterises the degree to which the diffusion is subdiffusive, with α i → 0 describing extremely subdiffusive behaviour and α i = 1 leading to normal diffusion. We refer to α i as the anomalous exponent.
The Laplace transform of the memory kernel [see Eq.(3)] for the distribution in Eq. (4) has the following simple form [26], However, the Laplace transform for this function cannot easily be inverted. This makes it difficult to find the memory kernel K i (t) in the time domain in explicit form. As a result, it is not convenient to write the reactiondiffusion equation in the form of Eq. (1). Instead, we use the Riemann-Liouville fractional derivative, defined as [27] This has the following property prodivded that lim t→0 [27]. It is reasonable to assume that this condition applies to the functions on which the fractional derivative acts in this work, as these functions are related to particle concentrations.
Using these definitions, we define the re-scaled time η i = τ i ξ −2/αi . This re-scaling hints at the underlying subdiffusive nature of the transport. Taking the limit ξ → 0, Eq. (1) can be written as where is the total reaction rate, and where we define The history-dependent quantity M i (x, t) characterises the current number of particles of type i hopping away from location x per unit time. References [11,14,22] contain a more detailed discussion of the origins of Eq. (10). In a similar way to Eq. (1), the reaction and diffusion terms in Eq. (10) are coupled through exponential terms in the fractional derivative. We note that the evolution of the concentrations described by Eq. (10) is dependent on the history of the system and is therefore non-Markovian. However, for α i = 1 the Mittag-Leffler distribution Eq. (4) reduces to an exponential distribution, and one obtains the normal reaction-diffusion equation In this case the length of time since the last hop is irrelevant to the current hopping rate, and the process is Markovian. We note that Eq. (10) can also be used as an approximation to the dynamics when a general class of heavy-tailed waiting-time distributions is used, not just the Mittag-Leffler function [8,9,14].
C. Linear stability and pattern formation in the reaction-subdiffusion system Now that we have discussed the general formulation of reaction-subdiffusion systems, we consider the stability of their fixed points, and discuss conditions for pattern formation. From here on in, we will focus on systems which have two species of reactants.
A stable fixed point in a non-spatial system can become unstable with the introduction of diffusive proceses. As a result of the instability, a limited band of Fourier modes with non-zero wavenumbers can become excited, resulting in the formation of stationary patterns [2,16]. Turing [1] demonstrated this effect in reaction-diffusion systems with normal diffusion, such as the one described by Eq. (12), and established a condition for such an instability to occur.
A similar instability has been found in the reactionsubdiffusion system whose evolution is described by Eq. (10), see [20,22]. We denote the homogeneous fixed point of this system by ρ 0 , which is defined by f i (ρ 0 ) = 0 for all i.
Unless indicated otherwise the f ij are always to be evaluated at the homogeneous fixed point. We denote deviations from the homogeneous fixed point by for the Fourier transform with respect to the spatial coordinate, one finds the following dynamics for the Fourier mode with wavenumber k, to linear order in the perturbations. We have written p i = R − i ρ 0 /ρ 0 i for the per capita removal rate of species i at the fixed point. We have also introduced To keep the notation compact we have omitted the argument ρ 0 in Eq. (14); it is to be understood that the A ij are evaluated at the homogeneous fixed point.
In [22] a condition for the Fourier mode k to be unstable to perturbations was derived for a system with one subdiffusing species and one normally diffusing species. A generalised version of this result for the case where both species subdiffuse is A derivation of this formula is detailed in Appendix A 1; we use an alternative route to that of [22]. We note that a sign error was made in the calculation in [22]. This mistake has been corrected in Eq. (15), see also Appendix A 2.

III. CONSTRUCTION OF THE EFFECTIVE NORMALLY DIFFUSIVE SYSTEM
A. Definition of the system and stationary patterns In the following, we construct a Markovian system, with normal diffusion, which replicates the stationary properties of the reaction-subdiffusion system. Finding a Markovian system such as this allows us to define effective diffusion coefficients and to interpret the behaviour in the original system more easily. We also discuss the conditions under which the transient dynamics of the original system can be approximated by the timedependent behaviour of this effective Markovian system.
At first it may seem surprising that one should be able to replicate the features of a manifestly non-Markovian system with a Markovian system. However, suppose that, in the long-term, the reaction-subdiffusion system reaches a stationary state. This state may be patterned or homogeneous. Since the non-Markovian 'memory' effects in the system are only detectable in dynamic quantities, one might suspect that the stationary behaviour can be described without the use of a non-Markovian memory kernel or a fractional derivative. In this case, one ought to be able to write down a Markovian system which has the same stationary properties as the original non-Markovian system.
We propose a Markovian system of the form where we define The specific choice ofD i (ρ), which replicates the stationary properties of the original system, will be described below. The system in Eq. (16) has the same reaction terms as that in Eq. (10), but the anomalous diffusion term has been replaced by a cross-diffusion term. The (effective) diffusion coefficientD i depends on the local concentration vector ρ(x, t). That is, the transport of one substance can be affected by the presence or absence of particles of all types. It is important to note that the diffusion term in Eq. (16) is not dependent on the history of the system.
To motivate the proposed effective cross-diffusive dynamics and to define the effective diffusion coefficientŝ D i , we first focus on the condition for the instability of a perturbation with wavenumber k. In the system described by Eq. (16) this mode is unstable if A detailed derivation of Eq. (18) is given in [28,29]. We define the effective diffusion coefficients for the subdiffusive system by requiring that Eq. (18) is the same condition as Eq. (15). The following choice satisfies this requirementD Making this choice we obtain By choosing the effective diffusion coefficients in this way, we ensure that the Markovian system and the original subdiffusive system experience the Turing instability for the same sets of system parameters. Furthermore, the same sets of modes are unstable in the Markovian system as in the original subdiffusive system. However, having the same sets of unstable modes does not necessarily mean that the two systems converge to the same patterned states in the long-time limit. This is because in order for the steady-state patterns to be produced, the exponential growth of the unstable modes must be curtailed by the non-linearities in the reaction rates. This is not accounted for in the linear stability analysis. That being said, one can indeed show that the stationary patterns produced in the both systems obey the same stationary equation. We first note the following property of the Riemann-Liouville derivative where Γ (α, x) = ∞ x s α−1 e −s ds is the upper incomplete gamma function [30], which has the property lim x→∞ Using Eq. (22) in combination with Eq. (10) and the fact that the concentrations do not vary in time in the stationary state, one can deduce that the stationary state in the subdiffusive system obeys which, if one notes the definition ofD i (ρ) in Eq. (19), is the relation defining the stationary state of Eq. (16).
B. Critical ratio of the effective coefficients required for pattern formation So far, we have discussed conditions for the instability of perturbations with specific wavenumbers k. We now proceed to establish a condition for the onset of Turing patterns in the subdiffusive system. Such a criterion, a version of which has been written down before in [20], is made far simpler with the use of the effective diffusion coefficients in Eq. (17). Since the effective Markovian system has been constructed so that the criterion for the instability of a particular mode k is the same as in the original subdiffusive system, the criterion for Turing pattern formation in general will also be the same in both systems. We therefore work with the effective system.
Turing patterns are formed when a finite range of Fourier modes with non-zero wavenumbers is unstable but the k = 0 mode is stable, such that Eq. (13) is satisfied. In order to find a condition for the formation of Turing patterns, one identifies the value of k 2 for which the left-hand side of Eq. (18) is minimal. We refer to this value as the 'critical' Fourier mode. In order for any other mode to satisfy Eq. (18), this critical mode must also do so. One finds where we keep in mind that D ij = D ij ρ 0 . Substituting this result into Eq. (18), one obtains the following condition for the presence of Turing patterns [28,29] When the effective cross-diffusion terms are zero (D 12 = D 21 = 0) this reduces to the condition for Turing instability in regular reaction-diffusion systems [16], This reduction to a known result for diffusive systems underlines the fact that the effective diffusion coefficients encapsulate the net diffusivity in the reactionsubdiffusion system.
For the cross-terms to be zero, the death rates of both species must be independent of the concentration of the other, or the diffusion of both species must be Markovian. The choice of whether the criteria Eq. (25) and Eq. (26) contain a plus or a minus sign is determined by which of species 1 or 2 is the activator/inhibitor. This can be seen realising that the critical ratio of the inhibitor diffusion constant to that of the activator for Turing pattern formation to occur should reduce to a number which is greater than unity in the normally diffusing case [5,20].
If species 1 is the activator, then Eq. (25) and Eq. (26) should contain a plus sign.

C. Limit of large removal rates
Naturally, the Markovian and subdiffusive systems described by Eqs. (10) and (16) do not, in general, have identical dynamics. However, one can show that the dynamics of the subdiffusive system are replicated to a good approximation in the corresponding effective Markovian system, when the particle removal rates are large. In this limit, particles survive for only short amounts of time on average, so the time horizon for memory effects is also small and, consequently, the dynamics can be approximated as memoryless.
To show this, let be the per capita removal rate for particles of type i. When the death rate is sufficiently large, one obtains the following We refer the reader to Appendix B for further details of this approximation. Noting that , the right-hand sides of Eqs. (10) and (16) are found to be approximately equal in the regime in which the approximation Eq. (27) is valid.

D. Interpretation
We have shown that an appropriately constructed normally diffusive system can replicate some of the most important features of the original subdiffusive system, in particular the location of any Turing instability in parameter space, and the resulting patterns. This has allowed us to define effective diffusion coefficients for the subdiffusive system. In this section, we provide further insight as to why this is possible.
Before we do this however, let us first contrast the effective diffusion coefficientsD i (ρ) with the quantities D ij (ρ) [see Eq. (17)]. The coefficientsD i (ρ) appear inside the second-order spatial derivative in Eq. (16); they represent the proclivity of particles to hop away from a given location. The coefficient D ij (ρ) on the other hand describes how the transport of particles of type i is affected by a local gradient of the concentration of particles of type j. Which set of coefficients one decides to use is a matter of taste. The coefficients D ij (ρ) have the advantage of highlighting the cross-diffusive nature of the system while theD i (ρ) are more succinct. The D ij (ρ) are also useful for writing down the criterion for Turing instability, as was demonstrated in Section III B.
So far, we have introduced effective diffusion coefficients for subdiffusive systems in the context of Eq. (16), but we have not related them to the statistics of particle trajectories. We would now like to provide further interpretation, and establish how the well-known diffusive law x 2 = 2Dt can be obtained in the context of subdiffusive systems. As we will see this can be achieved by defining a suitable ensemble of particles, and the connection with the effective diffusion coefficients can be made.
To do this, we consider the following equation, in which θ is a parameter which may take values 0 or 1. Eq. (28) encapsulates the subdiffusive terms of Eq. (10), and, for θ = 1, the annihilation reactions when the removal rate p i is constant and uniform (as it would be the case in the homogeneous steady state). In the original system, subdiffusive motion and removal are the only processes an individual particle can undergo once it has been created. We assume that the initial condition at time t = 0 is given by P i (x, t = 0) = δ(x). Let us discuss possible microscopic processes described by Eq. (28). For θ = 1, the equation describes a system in which particles subdiffuse by drawing waiting times from the Mittag-Leffler distribution, and in which they are also subject to a constant removal rate p i . The solution P i (x, t) of Eq. (28) describes the density of particles of type i present at time t at position x. The total number of particles in the system, P i (x, t) dx = e −pit , decreases with time. We define C i (x, t) ≡ P i (x, t) /[ P i (x , t) dx ]; this is the probability density for the position of surviving particles, that is particles still present in the system at time t.
The subdiffusive law is obtained by examining the mean squared displacement of these surviving particles. To see this we define x(t) 2 surv = C i (x, t)x 2 dx. One then finds Eq. (29) is derived in Appendix C. We now turn to the the case θ = 0, which can be obtained from the scenario for θ = 1 by adding a term p i P i (x, t) on the right-hand side of Eq. (28). This indicates additional particle production. More precisely, one representation of the case θ = 0 is a dynamics in which any particle that is removed is immediately replaced by another identical particle, which draws a new waiting time. This reflects the behaviour in the stationary state of the full system in Eq. (1). In the stationary state the particle number at each location is constant, and R + i (ρ 0 ) = R − i (ρ 0 ). The per capita removal rate of particles of type i is given by p i = R − i (ρ 0 )/ρ 0 i , and removal and production balance. The density P i (x, t) remains normalised at all times for θ = 0, and we define x(t) 2 replace = P i (x, t)x 2 dx. It is this ensemble which leads to diffusive behaviour, This relation is again derived in Appendix C. From Eq. (30), one can read off the effective diffusion . This coincides with the earlier definition in Eq. (19).
In summary, if one includes the removal term in Eq. (28) (θ = 1), one finds that the ensemble of surviving particles generates subdiffusive statistics. If this term is not present (i.e. removed particles are replaced), standard diffusive statistics x(t) 2 ∝ t ensue, despite the non-standard transport term in Eq. (28)

A. Lengyel-Epstein model
To test the conclusions of Section III, we numerically integrate the reaction-subdiffusion-equation (10) and the corresponding normally diffusive system in Eq. (16) and then compare the numerical solutions. We do this for the example of the Lengyel-Epstein model [31], which was introduced primarily as a way of modelling the ClO − 2 -I − -MA reaction. This chemical system exhibits Turing patterns experimentally [32,33]. We use a simplified twospecies version of the full model, and focus on the case of one spatial dimension. The system involves two types of particles, labelled A and B, which each undergo a (sub)diffusion process, potentially with different anomalous exponents. Particles at the same location can also react. The reactions are described by the terms where ρ A and ρ B are the particle concentrations, and where a, b, c, d and N are positive model parameters. The first term in the definition of f A describes the production of particles of type A, and the first term in f B represents the production of B-particles. The remaining terms describe removal processes; the corresponding per capita death rates for the two species are given by The Lengyel-Epstein system has a homogeneous deterministic fixed point at ρ 0 Using Eq. (13), the homogeneous fixed point in the well-mixed system is stable so long as ca > 3 5 a 2 − 25b 2 d.

B. Numerical integration method
The numerical integration of Eq. (10) is made difficult by the fact that the evolution of the system depends not only on the concentrations at the present time, but also on the concentrations at all earlier times. This is due to the nature of the fractional derivative 0 D 1−α t , defined in Eq. (8). Further complications arise from the coupling of reaction and subdiffusion terms, viz. the presence of the exponential term inside the fractional derivative in Eq. (10). To our knowledge, most existing numerical methods for integrating reaction-subdiffusion equations have not taken this coupling into account [34][35][36].
The key to integrating Eq. (10) is in recognising that the fractional derivative can be expressed in discretised form using the equivalence of the Grünwald-Letnikov and the Riemann-Liouville definitions of the fractional derivative, which is valid when f (t) is continuous and f (t) is integrable [27,37], where [x] denotes the integer part of x and where α j is the generalised binomial coefficient The numerical approach is then based on first-order Euler-forward integration. We also discretise space and operate on a one-dimensional lattice with spacing ∆x and periodic boundary conditions. The Laplacian in one dimension, At each time-step the fractional deriva- ρi(x,t ) dt ρ i (x, t) must be evaluated for every site on the lattice. This means that it is necessary to keep a history of the concentrations for every point on the lattice. For efficiency, it is also convenient to keep the history of the integral that appears in the exponential. The 'short-memory' principle [27] allows us to simplify matters and speed up the computation. The definition of the fractional derivative in Eq. (8) indicates that times closer to the present contribute more than times further in the past. Using a 'cut-off' t cut in how far back one goes in time to calculate the fractional derivative ensures that the numerical evaluation does not slow to halt as t becomes large. One must be careful however to choose The numerical integration is carried out on a lattice with 101 nodes and with spacing ∆x = 1, which is much smaller than the wavelength of the patterns produced in the examples. The time step used in the Euler-forward integration is ∆t = 10 −4 . Further model parameters are given in the figure captions.
C. Numerical comparison of subdiffusive and corresponding normally diffusive systems

Stationary patterns
As discussed in Section III A, for fixed model parameters we expect the same stationary patterns to result from the reaction-subdiffusion system and the corresponding effective Markovian system in the long-term. The patterns in the Lengyel-Epstein model shown in the lefthand panel of Fig. 1 confirm this; the amplitude, periodicities and general shapes of the patterns in the original and the effective systems match. The similarity of the patterns is further evidenced in the Fourier spectra shown in the right-hand panel. We attribute the remaining differences to inaccuracies in the numerical methods. Fig. 2 shows that the transient approach to the patterned stationary state in the subdiffusive system differs from that in the effective cross-diffusive system.

Time-evolution for large removal rates
We also deduced in Section III C that the original subdiffusive system and the corresponding effective Markovian system have approximately the same dynamics when the particle removal rates are large. In Figs. 3 and 4, we verify this for the case in which one species subdiffuses and the other diffuses normally. If one component diffuses normally (α i = 1), the evolution equation for this component is identical in both systems. So we only require the subdiffusing component to have a large removal rate in order for the dynamics to be the same in both systems. Fig. 3 demonstrates the accuracy of the approximation for large removal rates. We compare the transient timeevolution of the concentration for the original system and for the effective Markovian system. This is done for two different locations in the spatial domain. Fig.  4 shows the agreement between the systems over longer time scales as patterns begin to emerge.

V. SUMMARY AND DISCUSSION
In a Markovian reaction-diffusion system, the quantification of the motility of a particular species is comparatively straightforward since one can easily define the diffusion coefficient, which can be related to the rate of change with time of the mean squared displacement for an individual particle. The situation in subdiffusive systems is more complicated. In this paper we have defined the effective diffusion coefficient for a reactionsubdiffusion system, and we are now able to compare the motilities of particle species. These effective diffu-  Fig. 1. We show δρA/N as relief. The patterns are fully formed at roughly t = 150 in the effective system but take until around t = 200 to form in the subdiffusive system. The system is initialised at the homogeneous fixed point, with a deviation of magnitude 10 at x = 51 at t = 0. The cut-off time used was tcut = 10.  sion coefficients are dependent not only on the anomalous exponents and the typical waiting times of particles, but also on the reaction rates. This in turn gives rise to cross-diffusive behaviour. Moreover, we have shown that, using the effective diffusion coefficients, the condition for Turing pattern formation in a subdiffusive system can be written in the same form as it would appear in a Markovian system. We have therefore demonstrated that the effective diffusion coefficients are indeed the pertinent quantities with which one can characterise transport effects in a reaction-subdiffusion system.
We emphasize that while the effective Markovian system replicates the behaviour in the subdiffusive system in the stationary state, it reproduces the transient dynamical behaviour only under special circumstances. Interesting time-dependent phenomena which are peculiar to reaction-subdiffusion systems, such as the failure of frontpropagation [38], are not necessarily reproduced fully by the effective system. That being said, we note that the concentration profiles in [39], which were reported to have unique character due to subdiffusion, can also be produced using an effective Markovian system of the type discussed in this paper. This can be seen through the fact that Eq. (11) in [39] is a special case of our Eq. (23) for the particular reaction scheme used in that paper.
In Sec. III D we showed that the mean squared displacement of particles can be thought of as increasing linearly with time even in a subdiffusive system. This view can be taken if one interprets diffusion as a phenomenon defined by an ensemble of particles, rather than by tracking the motion of individual particles. Different choices for the ensemble can lead to subdiffusive or regular diffusive behaviour. In the latter case, the diffusion coefficient, defined from the ensemble view, coincides with the one we obtained from analysing the Turing instability of the system with anomalous diffusion.
There exist many explanations for the observed crossdiffusive effects in real-life systems. In ecological systems for example, predators pursue their prey and prey avoid their predators [40,41]. Mechanisms leading to cross-diffusive behaviour in physical and chemical systems include electrostatic interactions, excluded-volume effects and complexation [42]. Our analysis suggests another possible mechanism: If a system exhibits memory effects due to anomalous transport, cross-diffusion can arise when the removal rate of one species of particle depends on the concentration of another species.
A common way of measuring diffusion coefficients in chemical systems experimentally is the Taylor method, which involves monitoring the spread of the concentration of a drop placed in a laminar flow. A Gaussian profile is then fitted to the data in order to estimate the diffusion coefficients [43,44]. In such an experiment, individual particles are not tracked and so the diffusive properties are inferred based on macroscopic statistical behaviour. Our work would suggest that cross-diffusive coefficients arising in experiments of this kind could possibly originate from, or at least be affected by, non-Markovian transport.
While we restricted most of the discussion to the case of subdiffusion with Mittag-Leffler distributed waiting times, we note that an effective normally diffusive system can be found for general non-Markovian systems described by Eq. (1); the calculations that were carried out starting from Eq. (10) could just as well be performed with Eq. (1). So effective diffusion coefficients can also be found for reaction-diffusion systems with other types of non-standard diffusion. at long times. Eq. (A2) can be rewritten in matrix form In order for the solution to be non-trivial, the matrix M k must be singular, that is Eq. (A4) determines the possible values that λ k may take. The determinant ∆ k (λ k ) can be written as We note that this expression does not reduce exactly to the one given in [22], see Appendix A 2. The mode with wavenumber k is guaranteed to be unstable if Eq. (A4) has a real and positive root.
For 0 < α i < 1, we have ∆ k (λ k ) ≈ λ 2 k for large real λ k , i.e. ∆ k (λ k ) is positive. Additionally, the function ∆ k (λ k ) is continuous. Hence, if ∆ k (0) < 0 there must be at least one real positive zero. In other words, ∆ k (0) < 0 is a sufficient condition for instability. This leads to Eq. (15). We refer to [22] for arguments as to why this is not only a sufficient condition but also necessary.

Correction to the Turing instability criterion in [22] and verification in simulations
In this Appendix, we point out a small (but consequential) error which was made in the calculation of the Turing instability in [22] and built upon in [20]. We also discuss some of the most important consequences.
In going from Eq. (29) to (31) in [22], a sign error was made in one of the terms. Eq. (31) in [22] should read It is the terms involving A 1 and A 2 in Eq. (A6) which have the incorrect sign in [22]. In [20], conditions for Turing pattern formation are given in Eqs. (24) and (31). These are a special case of our more general Eq. (25), in that only one subdiffusing species is considered. Given the sign error in [22], Eqs. (24) and (31) in [20] attract corrections as well.
Following [22], we denote species 1 as the subdiffusing species with per capita removal rate p and total reaction rate f , and species 2 as the normally diffusing species with total reaction rate g. We define γ = 1 − α and θ γ = , the ratio of the effective diffusion coefficientsD evaluated in the homogeneous steady state. The condition for Turing instability can be written θ γ > θ γ,c , where the critical value θ γ,c can take two forms, depending on whether the activator or the inhibitor is subdiffusing.
For a subdiffusing activator one has where as for the case of a subdiffusing inhibitor the condition reads These expressions are different from the ones in [20] in that the signs of A 1 and A 2 have both been inverted.
The sign error has important consequences to the conclusions drawn in [20]. In fact, with the error taken into account the phase diagrams shown in Figs. 2 and 3 in [20] change considerably. We present the corrected phase diagrams in Fig. 5.

Verification in simulations.
To the best of our knowledge the theoretical predictions of [20] have not been tested numerically. To verify Eq. (A8) we use the method described in the main text and in Appendix D. The version of the Lengyel-Epstein model used in [20] can be obtained from the setup described in Eqs. (31) and (32) by the re-scaling ρ A = ρ A N , ρ B = cρ2 N and d = b = 1, c = b , a = a (dashed quantities are the ones used in [20]). This yields An example of the outcome of the numerical integration is shown in Fig. 6. One finds patterned stationary states only when θ γ is greater than the critical value given in Eq. (A8). At the same time patterns can be observed below the threshold given in Fig. 2 of [20], confirming our version of the calculation of the phase diagram.
Importantly, these numerical results demonstrate that, in this model, one can observe Turing patterns for lower values of θ γ when the activator is subdiffusing than when both reactants are diffusing normally or when the inhibitor subdiffuses. The opposite is reported in [20].  (d)): Model parameters a , b and γ are as on the left, but η1 = 0.008, η2 = 0.01, resulting in θγ = 4. This is just below the threshold for pattern formation in Fig. 5. The perturbation decays, and no patterns are found.

Thus we have
pi(x,t )dt dτ .
We now assume that the removal rate for particles of species i at the homogeneous fixed point is large. So, we would like to find a series expansion for integrals of the form where M is a large dimensionless parameter. This can be done using a method analogous to Laplace's method [45]. We presume that the functions f and g can be expanded as Taylor series, which converge for 0 < τ < t, such that f (t − τ ) = f (t) − τ f (t) + 1 2! τ 2 f (t) + · · · , g (t − τ ) = g (t) − τ g (t) + 1 2! τ 2 g (t) + · · · . (B4) In our case, g (t) > 0. Noting that the series for the exponential function is absolutely convergent for any value, one then obtains We now note the following where γ(s, x) = x 0 t s−1 e −t dt is the lower incomplete gamma function. Expanding further in Eq. (B5) and re-organising terms, one obtains the following asymptotic series The lower incomplete gamma function fulfills following relation [30] γ (s, and keep a record of all previous values of these quantities, up to the cut-off.

Go to 2.
With regards to choosing the cut-off time t cut , a simple test as to whether the cut-off is suitably long is to evaluate e −pit 0 D 1−αi t [e pit ] and ensure that this agrees with the expected analytical result of p 1−αi i for long times. Using an infinite t cut would clearly be the most accurate choice of cut-off. That being said, we found that it was always possible to find a finite cut-off time, which greatly increased the efficiency of the calculations and did not interfere with the results. We tested in selected examples that identical results are obtained for all intents and purposes if no cut-off is used.
Additionally, we note that more sophisticated methods of evaluating the Grünwald-Letnikov derivative do exist [46][47][48], but we found that the above method was sufficient for our purposes.