Coloured Noise from Stochastic Inflows in Reaction-Diffusion Systems

In this paper we present a framework for investigating coloured noise in reaction-diffusion systems. We start by considering a deterministic reaction-diffusion equation and show how external forcing can cause temporally correlated or coloured noise. Here, the main source of external noise is considered to be fluctuations in the parameter values representing the inflow of particles to the system. First, we determine which reaction systems, driven by extrinsic noise, can admit only one steady state, so that effects, such as stochastic switching, are precluded from our analysis. To analyse the steady state behaviour of reaction systems, even if the parameter values are changing, necessitates a parameter-free approach, which has been central to algebraic analysis in chemical reaction network theory. To identify suitable models we use tools from real algebraic geometry that link the network structure to its dynamical properties. We then make a connection to internal noise models and show how power spectral methods can be used to predict stochastically driven patterns in systems with coloured noise. In simple cases we show that the power spectrum of the coloured noise process and the power spectrum of the reaction-diffusion system modelled with white noise multiply to give the power spectrum of the coloured noise reaction-diffusion system.


I. INTRODUCTION
One of the central challenges in mathematical biology is understanding mechanisms involved in development processes.Within the context of developmental biology, the emergence of large scale spatial structure, has often been theoretically explored through a common framework of deterministic partial differential equations defining reaction-diffusion systems [1][2][3].While current frameworks may explain a variety of phenomena in development, they can also suffer from over-simplification [4], with the additional caveat that finding theoretical models both consistent with mechanism-based knowledge and capable of predicting observed patterns is a highly complex task, suffering from both model and parameter fine tuning [5,6].
Many models that describe pattern formation assume parameters are constant; however, this deterministic assumption is not suitable for certain conditions.Some systems are better suited towards a stochastic approach.When a system is coupled to external and stochastic drivers, then the parameter values can change.The stochastic driving is often represented by stochastic parameters, that is a parameter which is drawn from a certain distribution at each time step or spatial point, and referred to as extrinsic noise below.This contrasts with most previous work on stochastic pattern formation, referred to as intrinsic noise, which assumes low copy number and it does not assume external drivers as the source of noise [4,[7][8][9][10][11].
The structure of intrinsic noise is often taken to be be highly constrained and in particular uncorrelated in time, leading to white noise representations.For instance, if the noise is due to low copy number induced dynamics, Gaussian white noise forcing emerges from the chemical Langevin equation approximation to the chemical master equation [12].Even with such constraints there is a rich diversity of noise-induced phenomena, such as spatio-temporal pattern formation [7,10], stochastic oscillations [8], metastability [4], waves [9] or enhanced oscillation amplitude [13].
However, when the source of the noise is extrinsic, other forms of noise are permissible.
In particular, the defining special properties of white noise may in general be relaxed and hence stochastic forcing can emerge with non-trivial temporal correlations, often termed coloured noise.Our objective is to show how extrinsic noise influences spatio-temporal reaction-diffusion patterns, in particular by developing power spectral methods to analyse the impact of coloured noise.
To proceed, we first note deterministic pattern formation reaction diffusion systems in biological applications with n interacting biochemical species take the form [1,3] in one spatial dimension, where x ∈ R n ≥0 is a vector of n chemical concentrations.The models we consider require the term f (x) = (f 1 (x), . . ., f n (x)) T to be a vector of rational functions.The rational functions describe the underlying reaction network between the species and D(∂ 2 x/∂s 2 ) with D = diag(D 1 , . . ., D n ) describes the diffusion of x.Let there be M chemical reactions in the spatially homogeneous system whose dynamics are described by ∂x/∂t = f (x).Each reaction is parametrised by a reaction rate k i > 0 and therefore we have a vector of reaction rates k = (k 1 , . . ., k M ) ∈ R M >0 .Throughout this paper we restrict our analysis to mass action kinetics [14], rendering f (x) a vector of polynomials and we use homogeneous Neumann boundary conditions where L is the domain length.Note that, although our work is only applied to systems of one spatial dimension, the theory can easily be extended to an arbitrary number of spatial dimensions.
Counter-intuitively, it has been shown that under certain conditions [1] diffusion can drive an otherwise spatially uniform stable state to instability.Such unstable systems, which can form stable patterns, such as stripes or spots, are called Turing Systems [3,15].To avoid the complications of bistable dynamics, and later stochastic analogues such as switching between steady states, we impose the constraint that the spatially homogeneous solution has a unique stable steady state, x * , such that f (x * ) = 0 [1].To find models which have this property we use techniques from real algebraic geometry [16] which provides simple tools, which ensure there exists only a single steady state in a model.Since Turing's work in 1952 many biological patterning systems have been suggested to be Turing systems [17][18][19].
Fundamentally, the application of a set of partial differential equation (PDE) models for a biochemical reaction system assumes that the species of interest are in high enough concentration to allow continuum modelling.By contrast, when the number of particles in the biological system is low, intrinsic stochasticity of the system must be included in the model [4,[7][8][9][10], which typically yields studies that investigate the impact of white noise.
However, as mentioned, stochasticity in biological systems can also arise from extrinsic noise and hence temporal variations in parameter values [20].Thus, instead we generalise the deterministic system to include stochastic parameter values.Specifically, we focus on the effect of stochastic fluctuations in the constant term, the "inflow" term, of the chemical reaction network f (x).As the extrinsic noise can arise from a vast number of different mechanisms, such as varying experimental conditions [21], it is largely free of microscopic constraints, especially the absence of temporal correlation.However, the impact of correlated external noise has, to the authors' knowledge, received little theoretical modelling attention.
Thus we proceed to develop a framework to study external coloured noise forcing of the above deterministic system for pattern formation in biologically motivated scenarios, analysing how temporal correlation in noise, as described by colour, impacts self-organisation properties.
The paper is organised as follows.In section II we introduce required notions of chemical reaction network theory to select a class of models relevant to our framework, i.e. those whose number of steady states is unaltered by changes in parameter values, temporal or otherwise.
We then introduce stochastic inflow parameters and describe their connection coloured noise.
In section III we highlight the impact of noise colouring on the spatio-temporal patterns formed by example of the Schnakenberg system.In particular, we discuss noise arising from stochastic subnetworks and varying experimental conditions.Our numerical results are discussed in section IV where we summarise the distinguishable differences between the various noise colours.

II. THEORETICAL RESULTS
In this section we introduce the theoretical background of this paper.In subsection II A we introduce chemical reaction networks and use results from real algebraic geometry to select models where changing a parameter does not effect the uniqueness of the steady state of the chemical reaction network.Networks that satisfy a unique steady state are required for our framework.In subsection II B we define stochastic inflows and describe their potential sources.In subsection II C we first show how the generic law of mass action reaction-diffusion system can be approximated to highlight the mathematical equivalence of external and internal noise in our framework.Next, we introduce power spectra as an analytic tool to study spatial and temporal pattern formation and, finally, we calculate a general formula for power spectra which includes the effects coloured noise.

A. Chemical Reaction Network Theory
A central aim of chemical reaction network theory (CRNT) is to describe the properties of a chemical reaction network from its reaction graph alone [22,23].One such property is the capacity for multiple steady states.Define a chemical reaction network by the multiset N = {S, C, R}, where S, C, R are defined below.We begin by embedding the network into n dimensional Euclidean space R n by associating a basis vector e i to each chemical species X i such that X 1 → e 1 = (1, 0, 0, . . . ) T , X 2 → e 2 = (0, 1, 0, . . . ) T and so on.Let S = {X 1 , . . ., X n } be the set of all chemical species in the network, then a generic reaction can be expressed as The constants α i and β i are stoichiometric coefficients which give information about how many molecules of X i are consumed and produced in each reaction.By letting X i → e i we can formulate a reaction vector describing the consumption, or production, of a species X i in a reaction If an entry of r is negative, then a species is consumed, whilst if an entry of r is positive, then a species is produced.
Denote the set of all reaction vectors in a network by R = {r 1 , . . ., r M }.To complete the description of a chemical reaction network in terms of sets and their embedding into Euclidean space we introduce the notion of complexes.Complexes are linear combinations of species which appear on the left or right hand sides of reaction vectors.In equation , where C 1 is the reactant complex and C 2 is the product complex.Let the set of all complexes be C = {C 1 , . . ., C l }.
The reaction network, N , is therefore a directed graph with vertex set C and a directed edge between vertices C i and C j if and only if C j − C i ∈ R, with the same embedding, X i → e i .Note that the description of the reaction network does not include any notion of rate constants k, however, many results in chemical reaction network theory include the reaction rates explicitly as a positive vector k In this paper we study the influence of noise on reaction-diffusion systems that are able to produce patterns in a well-defined parameter region.Critically, we will see that coloured noise is able to have both a constructive influence outside of this previously defined parameter region (i.e. the noise stabilises patterns where we would not expect them deterministically), as well as a destructive influence on patterns which normally would arise.In particular, we would like to exclude the intrinsically stochastic effect of stochastic switching, which occurs for a system that has multiple steady states in some parameter region.Stochastic switching describes the phenomenon that a chemical reaction network can jump from one steady state to another when subject to finite stochastic perturbations.To avoid this scenario we use a network structure based tool described in [16] which a priori excludes multistationarity.
Take a chemical reaction network N and embed its complexes and reactions into R n .
Then define the matrices where {b i } is the set of reactant complex vectors of N .Let diag(k 1 , . . ., k M ) be the diagonal matrix of reaction constants.Using the law of mass action the dynamics of the species concentrations can be described by where Note that the reactant complexes and the reactions are structural properties of the reaction graph.Further, let a ∈ R n and define its sign vector σ(a) = {−, 0, +} n by applying the sign function to each component of a.Therefore, the i th component of the sign vector σ(a) i = sign(a i ) ∈ {+, 0, −}.
The link from network structure to multistationarity is outlined in [16,Theorem 1.4] and it concerns the injectivity of the map f : x → dx/dt.If f (x) is injective then there exists a unique vector x * ∈ R n such that f (x * ) = 0.In other words the network is monostationary.
The conditions for an injective f (x) hold for all parameter values k = (k 1 , . . ., k M ).A corollary of [16,Theorem 1.4] Note, that this condition depends on the network structure, specifically, the spaces spanned by the kernel of the reaction matrix, m 1 , and the image of the source complex matrix, m 2 , but not on the parameter values.

B. Stochastic Inflow in Chemical Reaction Networks
We will now show how chemical reaction networks (CRNs) that satisfy injectivity (preclusion of multistationarity) can be applied to parameter stochasticity, in particular, to the case of stochastic inflows.
Much effort in the CRNT literature investigates the effect of inflows into the chemical system, especially when the network is changed to a "fully open system" in which every species has an inflow reaction [24][25][26][27].By contrast, here we change the nature of the model by generalising the inflow reaction rates to stochastic processes, rather than simply adding inflow processes as deterministic reactions.Sometimes these inflow reactions are the called"inputs" in engineering applications and the concentrations are the "outputs" and, therefore the study of chemical reaction networks with varying inflows can be rephrased into a study of the "input-output" relation as has been the approach in previous work [28].
Consider a chemical reaction network N .An inflow reaction for species X i is manifested in the reaction graph as Usually, the work in CRNT assumes the reaction rate k in to be constant in time (and space).
However, the zero complex, ∅, symbolises a process not included in the model such as the the production of X i by another network, which we call an "auxiliary network", or a mechanical addition of X i to an experimental setup.Both mechanisms are often subsumed into ∅ and usually approximated as constant or "perfect" inflow, however, they can exhibit intricate dynamics.We model the dynamics of "non-perfect" inflows by a stochastic process whose origin can be two-fold.
(a) Stochastic sub-systems: We assume that the inflow into the (deterministic) reaction diffusion system is provided by the output of another chemical reaction network (with a unique fixed point).When the number of particles in the auxiliary network is large the system will be in a steady state and the inflow rate k in will simply be proportional to the steady state value of a species of the auxiliary system.However, when the particle number in the auxiliary system is low, stochastic fluctuations cannot be ignored and, while still proportional to the concentration of a species in the subsystem, the actual influx parameter k in will be a stochastic process K in (t).This again renders the influx parameter k in into a stochastic process K in (t).
Whereas the two sources of parameter noise are conceptually different, their mathematical description is the same.Consider a stochastic process K(t) with an underlying distribution D(t).Therefore, at every time t we have K(t) ∼ D(t).To be biologically relevant (i.e. to ensure all chemical concentrations are non-negative at every t) we require D = 0 for all K(t) < 0. To simplify the following mathematical analysis we approximate the distribution of K(t) as a truncated Gaussian, since we would like to connect our analysis to the internal noise case.Therefore, we can describe K(t) as such that: where ξ ∼ N (0, 1), Ω is positive constant and g(t, t ) is a positive function.We assume that the standard deviation of the Gaussian perturbation to the mean parameter value is small such that K(t) < 0 is a rare event.Due to the negligible truncation, we approximate the moments of the truncated Gaussian as the moments of the full Gaussian such that In the remainder of this paper we will show how various functions g(t, t ) can arise in mathematical modelling, especially when stochastic auxiliary networks are considered.
We can further connect the parameter σ to physical quantities of the underlying chemical reaction network by again considering the sources of stochastic inflow.We assume that the origin of the stochastic inflow is stochastic auxiliary networks or experimental imperfections.
Therefore, in the limit of either a perfect experiment or a deterministic auxiliary network the inflow should be deterministic and ( For stochastic auxiliary networks the quantity Ω which parametrises the size of the stochastic fluctuations can be related to the system size [29].Similarly to the system size expansion, we require Ω to be sufficiently large.Large Ω helps to maintain a positive solution of the stochastic partial differential equation, although, positivity cannot be guaranteed when the noise is correlated. Next, consider a chemical reaction network N in which a subset of species has an inflow reaction.Denote the set of species with (stochastic) inflow by S in ⊆ S. If there is more than one (stochastic) inflow |S in | > 1 the stochastic process will be multi-dimensional Brownian motion.The description of the stochastic process outlined in this section still applies in this case, if we have K in (t) as the vector of stochastic inflow such that and B representing the potential covariance between the stochastic inflows.
Since our assumptions made the inflow process additive and uncoupled to the species, the stochastic reaction network can be described by the system of stochastic partial differential equations (SPDEs) where η is a vector of stochastic processes such that its support is supp(η) = S in .
Note that in our limit randomising inflows leaves N as well as f (x) invariant, and hence, the steady state structure of the system; i.e. if we start out with an injective function f (x) then adding stochasticity will not change this injectivity.

C. Power Spectra for Stochastic Inflows
To fully classify the patterns arising from the addition of stochastic inflows we calculate the power spectra of the patterns.Power spectra are an analytic tool showing which spatial and temporal frequencies are present in a pattern [7,30,31].Peaks in power spectra give information about dominant frequencies and, hence, about oscillatory behaviour of a system in space and time.
First, we linearise equation ( 9) about the fixed point x ∼ x * + δx.Where δx represent small perturbations; these should decay in the deterministic limit as the system will converge to the steady state outside the Turing parameter regime.Inside the Turing regime the system is assumed to converge to a stable pattern rather than a stable state which is not considered in this paper.Hence, our analysis is valid only outside the Turing pattern regime.
Substituting the linearisation ansatz into equation ( 9) and keeping lowest order terms only where J| x=x * is the Jacobian of f (x) evaluated at the fixed point x * , which for notational convenience we will denote as J.
We make one further assumption, namely the scaling of the perturbation δx with the parameter Ω controlling the magnitude of the stochastic input.Let δx = Ω α δx.Then, if α < −1/2 and in the limit of Ω → ∞ the leading order term is the stochastic process η only, similarly, if α > −1/2, then η could be neglected to leading order in Ω.Hence, we let α = 1/2 and equation (10) simplifies to Note, that equation ( 11) is mathematically equivalent to a chemical Langevin equation of compartmentalised diffusion in the limit of the compartment size, ∆ s , going to zero [32].
However, in our derivation the source of the noise is external.To emphasise the mathematical connection with internal noise we will in fact discretise equation (11) into a finite number of compartments, such that for an n species network with K compartments we have and the matrices D and J turn into block matrices such that The ij th block of J is J ij I K×K to give Compartmentalisation results in a system of nK coupled SDEs with ).To calculate the power spectra of the system we introduce the discrete cosine transform [33] The cosine transform incorporates the von Neumann (no flux) boundary conditions, which are commonly used for reaction-diffusion systems, however, other boundary conditions can easily be implemented [33].Due to the boundary conditions we require κ = mπ/l with m ∈ {0, 1, 2, . . .}.We refer to m as the spatial mode.Note that the use of (j − 1) instead of j which is due to the fact that the compartment labelling starts at j = 1.Hence, by applying the spatial Fourier transform we reduce the system of nK SDEs (15) to a system of n coupled SDEs Finally, applying the temporal Fourier transform we get an expression for δ xκ (ω).Note that the Fourier transform always exists if we consider the system to have a fixed point [34].
Therefore, the power spectra of the pattern, P δx (κ, ω), can be expressed as the diagonal elements of where we defined and † denotes the hermitian conjugate of a matrix.In order to be guaranteed real power spectra we need N = N † .In the case of Gaussian noise we have where F cos denotes the cosine transform and F the temporal Fourier transform.Note that when all temporal correlations are equal such that G ij (t, t ) = g(t, t ) the transform factor becomes where P Correlations is the power spectrum of the correlation function g(t, t ) and B κ is the n × n matrix of covariances.For white noise we take g(t, t ) = δ(t − t )/dt to reduce (15) to the Itô form with η i (t)η j (t ) = B ij δ(t − t ) and, therefore, P Correlations = 1.Hence, we see that for white noise the power spectra are For general temporal correlations g(t, t ) we have such that the spectra of the noise colour appear as a multiplicative factor modulating the white noise spectra.

D. Application to the Schnakenberg System
We illustrate our analysis via the example of the Schnakenberg kinetics.The Schnakenberg system is an n = 2 species reaction-diffusion system which, despite its apparent simplicity, exhibits a wealth of different behaviours in the deterministic [35][36][37][38] as well as stochastic [7,31] regimes.The Schnakenberg system is a fully open system and, hence, both species can be subject to stochastic inflow.The system follows the reaction scheme The dynamics of the Schnakenberg system is described by the system of PDEs First we check for a potential multistationarity in the Schnakenberg system by formulating the matrices m 1 and m 2 .
Next, we consider stochastic inflows which result in k i → k i + 1/ √ Ωη i for i ∈ {1, 2} and the equations Linearising equations ( 31) about the steady state and discretising space we get a system of linear stochastic differential equations similar to the ones arising from the study of internal noise [7] with where a, d are tridiagonal matrices with diagonal elements and off-diagonal (sub-and super-diagonal) elements The matrices b and c are diagonal matrices with entries Note that the vector x has 2K components.There are K components for species x 1 and K components for species x 2 .
Hence, Fourier transforming equation (32) in space and time we can compute the power spectra, with In the following section we will discuss the effect of various noise correlations on the power spectra, and hence, the patterns generated.

III. COMPUTATIONAL RESULTS
In this section we highlight the computational patterns generated by a Turing system with coloured noise.In subsection III A we start with a brief discussion of the numerical methods used before we discuss the differences between parameter noise as considered in this paper and internal noise in subsection III B. We then review systems with white noise which were mainly derived in [7] and in the following subsections we proceed to relax the assumption of uncorrelated noise to see the effect of noise colouring.This section is divided into two parts, the first one treating stochastic processes which may arise from experimental set-ups, such as exponential broadening of the noise correlation (Ornstein-Uhlenbeck noise) or excitation of the zero frequency modes.The second part will treat sources of noise due to stochastic subnetworks and noise mixtures.An overview of our main results is given in

A. Numerical Methods
The system of SDEs ( 15) is simulated by an Euler-Maruyama scheme [39] with time step ∆t such that The important step of the integration is to find the correct vector η.
The white noise and Ornstein-Uhlenbeck noise are generated by an auxiliary noise process.
In particular, at teach time step the white noise is sampled from a multivariate Gaussian distribution η∆t ∼ N (0, B∆t).
The Ornstein-Uhlenbeck process is a generated by the stochastic differential equation where ξ is a standard white noise vector.Therefore, at time t the vector η(t) is added to the to the system of SDEs.The Ornstein-Uhlenbeck "auxiliary equation" is integrated by an Euler-Maruyama scheme as described in [40].
To simulate power law noise for which, in general, no auxialiary SDE exists we use inverse transforms [41].To generate a vector η(t) with correlations η(t)η T (t ) = Bg(t − t ), we first use the fact that η = √ Bξ where ξ may be correlated in time but not in space, i.e. ξ(t)ξ T (t ) = δ ij g(t − t ).Then, let the power spectrum of ξ i be 1/ω α as described in subsection III D and use the algorithm of [41] to create a time series.Multiplication of ξ(t) with √ B gives the desired noise process, η(t), which can be added to the SDE system.
Auxiliary networks as in III E can be simulated by either method and in this paper the auxiliary network input is generated by using the inverse transforms technique.

B. White Noise
First, we consider external white noise.In this special case the noise vector η is just a Wiener process with correlation matrix B. This case is mathematically identical to the case of internal noise in the weak noise limit as studied in [7].The main differences between internal noise and the parameter noise considered in this paper are the amplitude of the noise and the exact forms of the covariances of the stochastic processes η.
Both approximations (the weak noise limit, as well as our "truncated Gaussian" approximation) assume that the stochastic effect is a perturbation to the deterministic limit.
However, as derived in [7], the covariance matrix of the stochastic processes η is determined by the microscopic processes whereas for external noise, whose origin can be manifold, the covariance matrix is be arbitrary.For mathematical simplicity we choose the covariance matrix Throughout the remainder of this paper we fix the parameter values to with diffusion coefficients which results in steady state concentrations of x * 1 = 1, x * 2 = 1000.The parameters chosen are not generic, but they represent a particularly interesting point in parameter space.The parameters are in the (stable) oscillatory regime of the Schnakenberg system and outside the Turing space.The main function of the noise will be to temporarily move the system into the Turing regime.Note that the system has large inflows (and outflows) compared to the chemical reaction parametrised by k 3 which allows for larger variations in the noise translating into larger pattern amplitudes.The modifications to the power spectra due to coloured noise are generic and apply all points in parameter space.However, the visibility of these modifications depends on the point in parameter space and the correlation matrix.
Similarly, the results on amplitude of the patterns are parameter-dependent.To simulate the system, we discretise the space into 40 compartments of width ∆ s = 0.0025 (which gives a total domain length of L = 0.1).
Simulating equation (32) under the influence of white noise and with the given parameter values we obtain Figure 2 and its corresponding averaged power spectrum.The power spectrum shows the temporal frequencies, ω, and spatial modes, m, present in the pattern.FIG.2: An overview of the deterministic and white noise behaviour of the system with parameters as in (41).The peak in 2c at zero spatial mode, m = 0, indicates temporal oscillations which result from deterministic limit cycle oscillations.Due to the choice of reaction kinetics the temporal oscillations of species x 1 and x 2 are in phase.
The amplitude of the oscillations is about 3% of the steady state value for species x 1 whereas for species x 2 it is much lower.Therefore, we assume that any patterning of x 2 is not generally measurable and, therefore, we focus our attention on x 1 .The prominent temporal oscillations with no spatial dependency manifest themselves as a sharp peak in the power spectrum.This is due to the fact that the chosen parameter point is in the oscillatory regime of the Schnakenberg model.Upon close inspection slight spatial variations in the pattern of x 1 can be seen, but with white noise these are not very pronounced.In the following sections we show how noise colour can enhance spatial modes, create additional oscillations or even create a stable pattern.

C. Ornstein-Uhlenbeck Noise
Next, we investigate the effect of exponentially correlated noise also known as Ornstein-Uhlenbeck noise [42].This stochastic process has a finite correlation time τ which we interpret as a response time.Consider an experiment in which the inflow rate is highly sensitive to the ambient temperature.Further, suppose this temperature undergoes random fluctuations about a regulated mean value.Therefore, the correlation time τ could represent the average response time of the temperature regulator.The correlation time τ could also represent an intrinsic relaxation time of the system such as found in motility induced phase separation of bacteria [43].
We make the simplification that all the temporal correlations in equation ( 8) originate from Ornstein-Uhlenbeck processes with the same correlation time, τ , such that Hence, it follows that where B κ is the same matrix as in the white noise case and therefore P (k, ω) Ornstein-Uhlenbeck = P (k, ω) white /(ω 2 τ 2 +1).Hence, the noise colouring will dampen temporal frequencies of ω = 0 and, therefore, stabilise the pattern.Consequentially, in systems with Ornstein-Uhlenbeck noise we would expect any present stationary patterns that may be obscured by transient noise effects to be more prominent.We simulated the Schnakenberg system with Ornstein-Uhlenbeck noise and plotted the resulting patterns and power spectra in Figure 3.The pattern of x 2 has the same characteristic as in the white noise case, except for a smaller amplitude.Interestingly, x 1 shows a very different behaviour to the white noise case.Clear spatial structures can be seen, in particular the phenomenon of polarity switching [7,31].
A Turing pattern of mode m = 1 is generated with a given polarity, namely, a minimum at s = 0 and a maximum at s = L or vice versa.Polarity switching describes the inherently stochastic phenomenon of a sudden change in polarity as can be observed in Figure 3.The temporal oscillations, although still present, are reduced to a practically unobservable level.
We attribute this to the fact that the exponential noise correlations dampen the oscillations present in the white noise system.

D. Power Law Noise
By power law noise we mean a stochastic process whose frequency distribution follows a power law in Fourier space Various types of power law noise are common in engineering and are defined by colours summarised in Table II.In this paper we study red noise, white noise (subsection III B) and violet noise to illustrate the effects of a positive, zero and negative exponent, α.Red noise amplifies small temporal frequencies with 1/ω 2 → ∞ as ω → 0. Therefore, the zero frequency behaviour dominates the pattern.As can be seen in Figure 4, due to the amplification of small ω, the pattern becomes stable, even though the amplitude grows FIG.3: The patterns and power spectra for the species x 1 (left) and x 2 (right) with noise generated by an Ornstein-Uhlenbeck process and τ = 100.Spatial patterns become visible in the Ornstein-Uhlenbeck case which can be attributed to its dampening effect on temporal oscillations.This can be observed as a visible excitation of the m = 1 mode in the power spectrum of x 1 with the dashed line representing the analytical prediction.The well-known phenomenon of polarity switching [7] is observed in the pattern of x 1 .The spectra were averaged over 50 repetitions with T = 1000 and subsystem size Ω = 100.beyond the biologically viable (non-negative) bound indicating that non-linear effects need to be considered in the model.To investigate this further, we look at the noise vector η.
For the simulation times chosen, T = 200, the noise vector is independent of the time t such that each η i has a constant, but random value.Therefore, the effect of red noise is a deterministic one and to get an accurate description of the system behaviour a full nonlinear, deterministic compartment model with generic (random) inflow parameters needs to be analysed.The red noise case is similar to models which assume a cell-to-cell variability where a parameter is perturbed by a constant, sampled from a given distribution [21,44].In contrast to conventional Turing theory in which the inflow parameter is constant throughout the domain, our analysis indicates that randomly varying inflows across the domain might actually facilitate pattern formation.
When investigating the power spectra we see that for species x 1 the simulations match the mathematical prediction.For x 2 , however, we see that the peaks at the deterministic oscillation frequency have different heights.This indicates that either the oscillation amplitude is underestimated by the analytical prediction or, vice versa, the simulations exaggerate the oscillation amplitudes due to numerical limitations.Due to the agreement of prediction and simulations at all other points we disregard the different peak heights as numerical artefacts.
To illustrate the behaviour of the system on the other end of the colour spectrum we simulated violet noise which stabilised the temporal oscillations further and due to large power at large ω drove the oscillation amplitude in the linear treatment to non biologically viable values.Again, this indicates that violet noise cannot be treated in biological applications with a simple linear theory, but a full non-linear theory must be used.For violet noise the effect of the truncation of the Gaussian became too large to get an accurate prediction of the power spectra as can be seen in Figure 5. Further, while noise colours with positive exponent are actively studied in biology [45], the potential emergence of blue or violet noise in biological systems is not clear.

E. Stochastic Auxiliary Networks
We now proceed to the second major source of random inflows, namely the dependence of a stochastic auxiliary network.An auxiliary network is a chemical reaction network which provides an input into the Turing system (see Figure 6).The auxiliary network is connected to the main network only via an inflow reaction and, if the auxiliary network reaches a steady state, it can be subsumed into the zero complex for practical modelling.If, however, the auxiliary network exhibits more complex dynamics such as deterministic or stochastic FIG.4: The patterns and power spectra for the species x 1 (left) and x 2 (right) with noise generated by a red noise process.As predicted, the spectra are dominated by the behaviour the ω = 0, which amounts to the stabilisation of a particular spatial mode.The species concentrations, however, go negative indicating that a full non-linear model needs to be used.The peak in the power spectrum for x 2 is a numerical artefact.The spectra were averaged over 50 repetitions with simulation time T = 200.The subsystem size was Ω = 5000.oscillations, then it needs to be treated as a part of the main network.We focus on auxiliary networks which are deterministically stable, but show stochastic quasi-cycles.The dynamics of such systems are modelled by using a correlated stochastic inflow parameter as outlined in section II B. As an illustrative example we use the predator-prey system described in [8], In [8] it was shown that in the deterministic regime the system (45) has exactly one attractor for any choice of rate constants.Hence, as the inflow parameters k 1 and k 2 are proportional to the concentration of Y 1 , they will have a constant value.However, when the copy numbers of Y 1 and Y 2 are small and stochastic fluctuations are important the predator-prey model ( 45) can exhibit so-called "quasi-cycles" which are stochastic analogue of limit-cycles and manifest themselves as peaks in the power spectra.
Consider this "predator-prey noise" in a Schnakenberg system.Suppose the inflows k 1 and k 2 in the Schnakenberg system depend on the presence of the chemical species Y 1 , such that k 1 ∝ y 1 and k 2 ∝ y 1 , where y 1 denotes the concentration of Y 1 .Again, we assume g ij (t, t ) = g(t, t ) for mathematical simplification.In Fourier space the power spectrum of where we use the parameter values from [8] α = 0.000384, β = Γ = 0.04, Ω 2 0 = 0.016.Equation (46) has a peak at ω ∼ 0.06, so we expect peaks at a similar frequency for each non-zero spatial mode in the power spectrum of the Schnakenberg system.The resulting patterns and the corresponding average power spectra are presented in Figure 7.The power spectra gained a second peak in the m = 0 mode and peak in all modes m > 0 at ω ∼ 0.06.Therefore, it can be seen that a deterministic parent system can inherit the dynamics of a stochastic auxiliary network and mix it with its own intrinsic dynamics.

F. Mixed Noise
Next, we investigated the case of mixing noise processes, in particular, we consider an Ornstein-Uhlenbeck process which has a different correlation time for x 1 and x 2 , τ 1 and τ 2 respectively.Then we let τ 2 → 0 in order to recover the white noise case for species two.
We introduce the auxiliary normal stochastic processes ξ and hence, the spatial modes of the noise η κ are described by [42] where b κ is a 2 × 4 matrix satisfying b κ b T κ = B κ .We Fourier transform (47) and define the matrix, to give Letting τ 2 → 0 and computing the covariance matrix gives where B k,ij are the ij th elements of the B κ matrix.Note that N is hermitian and therefore we expect real power spectra for the patterns of x 1 and x 2 .
Simulating equation ( 11) with the noise process described by (47) gives rise to the mixed patterns seen in Figure 8. Note, that computationally the limit τ 2 → 0 is equal to setting τ 2 to the time step dt.The patterns of both species appear similar to the pure Ornstein-Uhlenbeck patterns, however, with reduced amplitude.

G. Reduced Stochastic Inflows
In the final subsection we instigate the case of only species x 1 being subjected to stochastic inflows.The species x 2 is assumed to only have deterministic inflow at a rate k 2 and η 2 = 0.
Hence, the correlation matrices N will be reduced to In the noise processes studied in this paper the behaviour of the species is virtually unchanged and the patterning is robust with respect to disregarding cross-correlations.

IV. CONCLUSION
In this paper we investigated the effect of stochastic inflows on a deterministic reactiondiffusion system.We restricted the class of networks considered to monostationary systems and used results from real algebraic geometry to show how monostationarity is related to The power spectrum of x 2 is a mixture of White noise and Ornstein-Uhlenbeck noise.The power spectra were averaged over 50 simulations with simulation time T = 1000 and subsystem size Ω = 5000.network structure.We then introduced a stochastic perturbation to an inflow reaction as a truncated Gaussian process with the expectation value at the deterministic inflow parameter.
After linearising we computed the power spectra for arbitrary noise colours and showed that in simple cases the power spectra can be derived as a multiplicative factor.The remainder of the paper consisted of applying our analysis to the Schnakenberg system and highlighting the effects various noise colouring can have on a reaction-diffusion system.
We briefly restated the results from the white noise analysis, proceeding to add coloured noise and then demonstrating its effect by computing the power spectra.In particular, for the simple case when all species experience the same temporal correlations, the power spectra of the correlations appear as a multiplicative factor in the total power spectra.Hence, depending on the nature of the correlations they will suppress or excite temporal frequencies at all spatial modes.Ornstein-Uhlenbeck noise has a Lorentzian frequency distribution and, therefore, suppresses positive frequencies.Power law noise can either completely stabilise or destabilise the pattern depending on the sign of the exponent, α.This is due to the fact that for positive α temporal frequencies at = 0 go to infinity and, therefore, all oscillations are suppressed which results in a stable pattern which resembles a pattern arising from simulations inside the Turing regime.
For the simulation times used in this paper, the noise vector was actually constant, and, therefore, deterministic methods could be used to study the stabilising effect of pink or red noise.The opposite is the case for negative α where the frequency contribution in the power spectrum increases as ω increases and oscillations are enhanced.However, certain aspects of the red, blue, or violet, noise cases cannot be explained by our simple analytic prediction.
The biological significance and the correct analytic description of blue/violet noise could form a part of further work.
When turning to stochastic auxiliary networks we see that the reaction-diffusion system can inherit traits of the dynamics of the auxiliary network.In particular, we showed that when the auxiliary network exhibits a predator-prey dynamic with stochastic oscillations these oscillations can still be observed in the deterministic main network.We concluded the range of applications by considering mixed noise with Ornstein-Uhlenbeck input for species x 1 and white noise for species x 2 .In this case, the system behaves similar to a system with pure Ornstein-Uhlenbeck dynamics and differences can only be found by a power spectral analysis.Finally, we observed in the cases considered that when only one species experiences stochastic inflow the patterns created are, except for potential special cases, similar to the ones when the inflow to both species is randomised.However, the extent to which such results may hold in generality is for further work.
Further directions could also include attempting to relate these studies to potential mechanisms of left-right symmetry breaking amplification in developmental biology, in particular the impact of induced Nodal production on one side of the embryonic node, which is hypothesised to be driven by ciliary fluid flows and also highly error-prone [46].It is generally asserted that the resulting interactions of the gene products Nodal and Lefty, which are major contenders as Turing morphogens [47][48][49], amplify this initial error-prone signal to generate robust patterning, driving downstream developmental left-right asymmetry.However, the ability of Turing systems to amplify the spatially localised, error-prone, and thus stochastic, influx of activator morphogen [46,Figure 2] to induce a robust symmetry breaking self-organisation, as well as any additional constraints required to do so, is theoretically untested.Thus examining the mechanistic basis of these postulates in this critical developmental biological process provides a fundamental application for the theoretical foundations developed here.

FIG. 1 :
FIG. 1: The two mechanisms contributing to stochastic inflows.The boxes on the left are usually treated as black boxes resulting in some constant inflow modelled by the parameter k in .We differentiate between experimental fluctuations, symbolised by a correlation τ in (a) and auxiliary networks described by f aux (y) in case (b).
The deterministic solution of the Turing system at the given parameter point The power spectra of the analytical calculations (dashed lines) and their simulated curves averaged over 50 realisations of length T = 1000 and subsystem size Ω = 5000.We attribute the differences in peak height to the finite time steps used.

6 FIG. 5 :
FIG.5:The patterns and power spectra for x 1 (left) and x 2 (right) with noise generated by a violet noise process.The discrepancy in the power spectra originates from the truncation of the Gaussian process.The negative species concentrations indicate that a nonlinear theory needs to be used to fully model a violet noise system.The spectra were averaged over 50 repetitions with simulation time T = 200.The subsystem size was Ω = 5000.

8 FIG. 7 :
FIG.7:The power spectra for the species x 1 (left) and x 2 (right) with noise generated by a stochastic predator-prey network[8].The power spectra inherited a second peak from the subsystem which activates oscillatory modes on top of the deterministic oscillation frequency.The spectra were averaged over 50 repetitions and simulation time T = 200 with m indexing the spatial mode.The subsystem size was Ω = 100.

8 FIG. 8 :
FIG.8:The patterns generated for the species x 1 (left) and x 2 (right) with Ornstein-Uhlenbeck noise and τ = 100 for x 1 and x 2 with white noise.The patterns are phenomenologically similar to the Ornstein-Uhlenbeck patterns, however, with reduced amplitude.

TABLE I
: An overview of our main results on various noise colours.

TABLE II :
The various colours of power law noise.In this paper we focus on red and vio-