Gaussian approximations for stochastic systems with delay: chemical Langevin equation and application to a Brusselator system

We present a heuristic derivation of Gaussian approximations for stochastic chemical reaction systems with distributed delay. In particular we derive the corresponding chemical Langevin equation. Due to the non-Markovian character of the underlying dynamics these equations are integro-differential equations, and the noise in the Gaussian approximation is coloured. Following on from the chemical Langevin equation a further reduction leads to the linear-noise approximation. We apply the formalism to a delay variant of the celebrated Brusselator model, and show how it can be used to characterise noise-driven quasi-cycles, as well as noise-triggered spiking. We find surprisingly intricate dependence of the typical frequency of quasi-cycles on the delay period.


I. INTRODUCTION
Traditionally, chemical reaction systems are modelled by sets of differential equations, also known as rate equations.These equations describe the time evolution of the continuous real-valued concentrations of the different particle types, and they are derived from the microscopic reactions using mass-action principles [1].This framework is mathematically convenient: the theory of ordinary and partial differential equations is well developed, and the tools for their analysis are readily available.Descriptions based on deterministic differential equations have one major drawback though, they systematically neglect all effects of stochasticity.Mathematically, deterministic descriptions are only appropriate for large (formally infinite) systems.
It is now commonly accepted that intrinsic stochasticity, arising due to the finite numbers of particles in reaction systems, can have significant effects on the dynamics [2].This includes noise-induced cycles [3], patterns and waves [4,5], and phenomena such as extinction and fixation [6].In order to fully capture these effects models must describe the reaction dynamics at the micro-scale and keep track of integervalued particle numbers.Assuming a well-mixed reactor, the state of the system is fully characterised by the number of molecules of each type present in the system at a given time.The time dependence of the probability distribution over the space of states is then governed by a chemical master equation (CME) [7,8].
The CME is in general difficult to solve, exact solutions only exist for a limited set of examples, such as simple one-step birth-death processes and those satisfying detailed balance [8,9].For the majority of cases approximative schemes present the best opportunity for analytical progress.Deterministic rate equations, as described above, are the simplest such approximation.They neglect all stochasticity, and formulate an effective dynamics for the first moment of the probability distribution over microstates.These equations can either be written down based on intuition, or they can formally be derived from the CME by means of the so-called system-size expansion.This requires the presence of a large parameter, typically the volume of the reactor or the total particle number.Its inverse is then a small quantity, and serves as an expansion parameter.The deterministic approximation is obtained from the lowest order of this expansion, corresponding to an infinite volume or particle number [8].Retaining the sub-leading terms in the expansion on the other hand ultimately leads to a stochastic differential equation with Gaussian noise.There are multiple methods by which the expansion of the CME can be carried out, e.g. the method by van Kampen [8] or the Kramers-Moyal expansion [10].Depending on the details one obtains either additive noise, this is referred to as the linear-noise approximation (LNA), or multiplicative noise.The resulting stochastic differential equation in the latter case is frequently referred to as the chemical Langevin equation (CLE).We will refer to any of these methods as Gaussian approximations, as the noise in the resulting effective dynamics is Gaussian, but it should be noted that the distribution of the quantities described by those effective dynamics itself may not be Gaussian in the case of multiplicative noise.The various different approximations are all closely related, see e.g.[11] for further details.The starting point for most of the expansion techniques mentioned above is the CME, which can generally only be formulated for Markovian systems.These are systems without memory, i.e. systems for which the dynamics at any one time only depends on the state of the system at that time, but not on the previous history by which the system has arrived at this state [9].This implies that all effects of chemical reactions must occur instantaneously, no reaction event has any further effects on the system at later times.
The times between reaction events follow exponential statistics, or equivalently the number of events occurring in a fixed time interval is Poissonian.These assumptions do not hold in a variety of applications.One of these is gene expression dynamics, in which there is a characteristic time delay associated with transcription and translation events.These may be triggered at a given time, but their products are only generated at a later time [12,13].Another example is in the modelling of epidemics, where an individual may be infected at a given time, and where recovery is at a later time, drawn from a distribution peaked around a typical infectious period.Recovery is then not an exponential process.We will refer to models of this type as delay models.The delay description can be considered an effective description at a coarsegrained level.On a much finer level a Markovian description of the underlying reactions may be appropriate.A key distinction is between models with a fixed delay and models with so-called distributed delay [14].In models with distributed delay the delay times are stochastic variables themselves, whenever a delay reaction is initiated a delay time is drawn from an underlying distribution, and delayed effects materialise at this later time.Fixed-delay models are the special case in which the delay kernel is a δ-function.
Non-Markovian systems can in principle be mapped onto high-dimensional Markovian systems [8].This procedure is not applicable to general delay kernels though, or it leads to Markovian dynamics with an infinite number of degrees of freedom.This is often unsatisfactory, and as a consequence many existing analyses of delay systems have focused on the deterministic limit [15][16][17].It is only more recently that a more systematic stochastic description of delay reactions has been attempted [18][19][20][21][22].This existing work has predominantly focused on extensions to the CME to account for delay reactions [18,20,21], but these equations typically do not close and are limited to specific classes of model systems.Systems with constant delay have been studied using the LNA in [18,20].
In a previous piece of work [22] we proposed the use of a method known from condensed matter physics [23] to describe the time evolution of discrete reaction systems with distributed delays.This technique, the socalled Martin-Siggia-Rose-Janssen-de-Dominicis generating functional [24] takes a path-based view.It considers the space of all possible time courses of the system and formulates the probability for a given path to occur as a dynamic generating functional or path integral, effectively representing the Fourier transform of the probability measure in the space of all dynamic paths.This is a powerful formulation applicable to a wide class of delay systems, and it can be used to derive effective Gaussian approximations, in particular an equivalent of the CLE for delay systems.The approach involves relatively complex mathematics though, and the purpose of the present paper is to show how this machinery can be bypassed.We present a heuristic, more intuitive procedure to derive the CLE for delay systems.To demonstrate its utility this method is applied to a variant of the celebrated Brusselator system with delay.In particular we focus on effects of stochasticity in parameter regimes in which the deterministic delay system approaches a fixed point.We are able to characterise the stationary distribution of the system, and the noise-induced quasi-cycles the dynamics generate.From our analytical calculations we find a surprisingly intricate dependence of the typical frequency of the cycles on the delay.We also study choices of parameters in which the deterministic delay dynamics constitute an excitable system.In the presence of intrinsic noise stochastically triggered spikes are found, we show how these can be simulated efficiently using the CLE.The remainder of the paper is organised as follows: In order to keep our paper self-contained and pedagogical we briefly summarise earlier work by Gillespie [25] in Sec.II.In particular we describe how the Gaussian approximation is obtained for Markovian systems, and how the corresponding CLE is derived.In Sec.III we then carry out a similar analysis for delay systems, and derive both the delay-CLE and subsequently the corresponding LNA.These results are then applied to the specific example of the Brusselator with delay dynamics in Sec.IV.We present a summary our work in Sec.V, and outline possible future work.In the Appendix we briefly describe how the well-known modified next-reaction method is modified to accommodate distributed delays, and we provide further supplementary details of our analytical calculations.

II. MARKOVIAN DYNAMICS
Before considering delay reactions it is useful to first summarise the main results for systems without delay obtained from the above expansion methods.We consider a well-mixed system composed of particles of different chemical species X α , α = 1, . . ., S. The corresponding particle numbers are written as n α ∈ N 0 .Interactions occur via a set of reactions, i = 1, . . ., R. Each chemical reaction, i, is written in the form The stoichiometric coefficient s i,α represents the number of X α particles entering the reaction, and q i,α is the number of such particles exiting.It is useful to define the quantities v i,α = q i,α − s i,α , so that v i,α indicates the change in the number of particles of type α when one reaction of type i occurs.The above notation indicates that the rate with which reaction i occurs is given by when the system is in state n = (n 1 , . . ., n S ).These rates are scaled with an overall parameter Ω, representing for example the volume of the system or a scale for the total number of particles.Each rate is of order Ω, such that the number of reactions occurring per unit time scales as Ω as well, in-line with standard conventions [8].The quantity r i (n/Ω) can be understood as the 'intensive' rate.The CME is then given by where P n (t) is the probability that the system is in state n at time t.The quantity T nn is the total transition rate from n to n , and similar for T n n .Next, it is convenient to introduce concentration variables x = n/Ω.Carrying out a Kramers-Moyal expansion of the CME in powers of Ω −1 up to and including sub-leading terms the following CLE is found (via a Fokker-Planck equation) [8] Here, η α is a white Gaussian noise term with mean η α (t) = 0 and with correlations across species given by The above result for the CLE is derived from a controlled expansion procedure, we do not re-iterate the full details here.Alternatively, the same result can be obtained from a heuristic argument following the lines of [25], which we summarise here.As a first step one discretises time into intervals of duration ∆.Assuming for the time being that reaction rates remain constant in each such interval, the number of reactions of type i firing in the time interval [t, t + ∆) is an integer random variable k i,t drawn from a Poissonian distribution with parameter R i (x t )∆.One then has The {k i,t } are statistically independent for different i and t, and, given their Poissonian statistics, we have and where we have used R i (x) = Ωr i (x).One finds (8) Eq. ( 8) together with Eq. ( 7) is recognised as the Euler-Maruyama discretisation of the continuoustime CLE: with where the noise is interpreted in the Itō sense.This has a slightly different structure to the CLE commonly found in the literature, which is as presented in Eq. ( 4) and Eq. ( 5).In Eq. ( 9) we have one noise source, ζ i , for each reaction -similar to the conventions found in [26].The more common choice is to have one noise source, η α , for each species.The two forms of the CLE are related via We shall use the form of Eq. ( 9), as it makes clearer the origin of the noise and is more straightforward to extend to delay reactions, both for analytical calculations and numerical simulation.In what follows we will only use the identification shown in Eq. ( 10) when we make the LNA.

III. GAUSSIAN APPROXIMATION FOR CHEMICAL REACTION SYSTEMS WITH DELAYS A. Definitions and notation
We will now introduce our notation for delay reactions.These are reactions triggered at one time, t, in the same way as conventional reactions without delay.Let us focus on one particular type of reaction, say i.As before we assume that the initiation of the reaction occurs with a rate R i [x(t)], and that this rate only depends on the state of the system at time t.Delay reactions can have an immediate effect on particle numbers, indicated by stoichiometric coefficients v i,α as before.At the time a delay reaction of type i triggers, a delay time, τ , is drawn from the distribution K i (•).We always imply K i (τ < 0) = 0 and At a later time, t + τ , a further change of particle numbers may then occur, indicated by a second set of coefficients, w i,α .The special case of fixed delay is recovered when K i (•) is a δ-distribution.Furthermore, reactions without delay are contained in this notation as well, one simply has w i,α = 0 for all α for such reactions, the kernel K i (•) is then irrelevant.Our formalism is technically limited to reactions with delayed effects at one single subsequent time, i.e. particle numbers may change at the initial time and at most at one later time.However it is possible to extend the method to include multiple delay periods, that is reactions for which events take place at more than two distinct times.

B. Derivation of the approximation
Similar to the Markovian case, we proceed by first discretising time into intervals of duration ∆.All times t and τ are then integer multiples of ∆.The total number of reactions of type i triggered between times t and t + ∆ is then a Poissonian random variable, k i,t , with parameter λ i,t ≡ R i (x t )∆ as before.Each of these k i,t reactions may have an additional effect on particle numbers at a later time.The delay periods are independently drawn from the kernel K i (•) for each occurrence of the reaction.The probability that a particular delay of τ ∈ ∆N + is drawn is ∆ × K i (τ ) in the discrete-time model.The rate with which a reaction of type i triggered at t with a delay of precisely τ is λ τ i,t ≡ (R i (x t )∆) × (∆K i (τ )), as the triggering of the reaction and the delay period are independent.The number of reactions of type i triggered at t and with a delay of precisely τ is hence a Poissonian random variable, k τ i,t , with parameter λ τ i,t .One has τ ≥∆ k τ i,t = k i,t , and the discrete-time kernels are normalised such that ∆ τ ≥∆ K i (τ ) = 1.The change of particle concentrations at time step t can then be written as The first sum on the right-hand side captures the instantaneous effects of reactions at the time they are triggered.The second sum represents changes of particle numbers occurring at time t, but resulting from delay reactions triggered at earlier times, t − τ , where τ = ∆, 2∆, . . . .This is illustrated in Fig. 1.Keeping in mind that k i,t = τ k τ i,t , we can write Eq. ( 11) as Similar to the procedure in [25], the Gaussian approximation now consists in replacing the Poissonian random variables k τ i,t by Gaussian noise variables with appropriate first and second moments.The mean of ... the distribution of k τ i,t is λ τ i,t , and the variance is also λ τ i,t .We therefore replace k τ i,t as follows: where the ζ τ i,t are independent Gaussian random variables of mean zero and with second moments = δ i,j δ t,t δ τ,τ λ τ i,t /(Ω∆).Recalling the definition λ τ i,t ≡ Ωr i (x t )K i (τ )∆ 2 this leads to where and where we have used in the limit of small ∆, this is recognised as the discretisation of the CLE (see Appendix A for further details) with the so-called drift term [9] F and where In this formalism reactions without delay have w i,α = 0 for all α, as mentioned before.In this case the noise variables ζ i (t, τ ) only enter Eq. ( 16) when the reaction fires, and we can define This CLE is the starting point for further analytical approximations, and it can also be used for efficient numerical simulation of the stochastic process in the limit of weak noise, see Appendix B for details.

C. Linear-noise approximation
The above CLE represents, in general, nonlinear dynamics with multiplicative noise.This makes analytical progress difficult, further simplifications and approximations are required.We note that the noise terms in the CLE are of order Ω −1/2 .The dynamic variables, x α , are stochastic variables and they enter nonlinearly in the drift term and in the noise amplitudes.This means nonlinear effects of the intrinsic stochasticity are retained, despite the Gaussian approximation of the noise.The LNA takes the approximation one step further, and linearises the effects of noise entirely [8].It retains nonlinearity only at the level of the deterministic limit; the time-evolution of fluctuations about this limit are described by a set of linear Langevin equations, containing additive noise only.The LNA is derived by separating deterministic and stochastic effects, viz.x → x ∞ + Ω −1/2 ξ, and by subsequently discarding all terms of quadratic or higher order in Ω −1/2 .The superscript ∞ here indicates concentrations in the deterministic limit, Ω → ∞.Applying this procedure to Eq. ( 16) one finds ẋ∞ x ∞ ] on the deterministic level, i.e. to lowest order in Ω −1/2 .At subleading order one finds that fluctuations are governed by where We have here collected all noise terms into one single Gaussian source, η α (t), given by compared to the expression in Eq. ( 18).As explained above, in the LNA the noise becomes additive, and the correlation properties of the {η α (t)} only depend on the deterministic variable x ∞ (t).The explicit form of the correlations of the Gaussian additive noise is given by η α (t)η β (t ) = B α,β (t, t , x ∞ ), where [27] While many analytical studies of these linearised dynamics focus on regimes in which the deterministic system reaches a fixed point (see e.g.[3]), it is important to keep in mind that the LNA can be derived for more general classes of model systems.In the general case the resulting Langevin equation has time-dependent coefficients, derived from the underlying deterministic trajectory.Progress can be made not only when expanding about fixed points, but also when the deterministic system approaches a limit cycle [28,29].Recently, expansion techniques have also been applied to chaotic systems in discrete time [30], although there is then relatively little scope for a full analytical characterisation.The LNA does have some limitations, and can not be used for all problems.For example it breaks down when the system is not monostable, when stochasticity induces transitions between attractors, and when non-linear effects are important (an example of which will be studied in Sec.IV D 3).

IV. STOCHASTIC BRUSSELATOR WITH DELAY A. Model definition
The Brusselator is a paradigmatic model for oscillations in chemical reaction systems [31][32][33][34].Previous work on stochastic versions of the Brusselator model include the study of fluctuations about limit cycles [28], pattern formation in spatial variants [5], and a recent study on noise-induced switching between large and small amplitude oscillations [35].
We here study a non-Markovian variation of the Brusselator model.This variant is not chosen with a particular real-world chemical system in mind, but instead in order to demonstrate how the techniques discussed in the previous section apply to a specific example.Our choice of model is motivated by the work of [36], who have studied delay reactions of the same form in the context of pattern formation in developmental biology.This existing work, however, is based on numerical simulations only, no analytical calculations are presented in [36].
The reactions by which we define the Brusselator with delay are the following: The reactions are indexed with i = 1, 2, 3, 4 from top to bottom.We here focus on a well-mixed reactor.
The first reaction is a delay reaction, and is the sole distinction between this model and the conventional Brusselator.The above notation indicates a two-step reaction: at the time the reaction is triggered two molecules of type X 1 react with one particle of type X 2 , to generate a particle of type Y .The X 2 -particle is removed in this process.Following mass-action principles, reactions of this type are triggered with a rate cn 2 1 n 2 /Ω 2 ,where c is a constant, and where n 1 and n 2 indicate the number of particles of types X 1 and X 2 in the reaction container respectively.As before Ω sets the scale of the system size.The intermediate particle Y remains for a duration, τ , drawn from a distribution K(•).The second segment of the delay reaction occurs τ units of time later, when the Y -particle decays into a particle of type X 1 .The remaining reactions are standard, and represent the creation and removal of particles of type X 1 (second and third reaction above, respectively), and the conversion of particles of type X 1 into particles of type X 2 (fourth reaction).

B. Chemical Langevin equation and deterministic limit
The CLE obtained from reactions Eq. ( 22) and using the above formalism are The noise correlators are Note that the subscripts of ζ i refer to the reaction number and the subscripts of x α refer to the species.We notice that the dynamics for x 1 and x 2 is independent of that for y, in the sense that the first two equations in Eq. ( 23) along with Eq. ( 24) constitute a closed set of equations.This is because the concentration of particles of type Y does not enter into any of the reaction rates.Given this decoupling, we will focus on the dynamics of x 1 and x 2 from this point forward.
In the deterministic limit, Ω → ∞, the noise terms in Eq. ( 23) vanish, and one obtains the following deterministic equations, As before, the superscript ∞ is used to indicate that these equations apply to concentrations evaluated in the limit of an infinite system.This deterministic system has a unique fixed point at identical to that of the conventional Brusselator system [32].For systems without delay it is straight-forward to characterise the stability of fixed points from the eigenvalues of the corresponding Jacobian.
Often the corresponding phase diagram can then be computed analytically.For systems with delay such an analysis is more involved.Delay kernels lead to transcendental equations for the equivalent of eigenvalues, and these typically have an infinite number of solutions in the complex plane.The assessment of the local stability for a fixed point of a delay system is hence mostly limited to numerical methods [37].

C. Linear-noise approximation
We will now proceed to make further analytical progress towards describing the stationary state of the stochastic delay Brusselator system.Our starting point is the general expression for the LNA for delay systems in Sec.III C. For the Brusselator system with delay one finds with noise correlators η α (t)η β (t ) = B αβ (t, t , x ∞ ) given by We focus on a parameter regime in which the deterministic system approaches the fixed point (x * 1 , x * 2 ) given in Eq. (26).At asymptotic times the LNA then consists of a set of linear Langevin equations with constant coefficients, which makes further analysis particularly straightforward.One has with From these expressions, a complete statistical characterisation of the stationary state of the linearised dynamics can be obtained.The power spectra, S αβ (ω) = ξ α (ω) ξ β (−ω) , of fluctuations about the deterministic fixed point can be found after Fourier transforming Eq. ( 29) and Eq. ( 30) [9,38].The elements of the matrix S(ω) are the Fourier transforms of the cross-correlation and auto-correlation functions of the {ξ α (t)}, C αβ (τ ) = ξ α (t)ξ β (t + τ ) .
The equal-time covariance matrix Ξ αβ = ξ α (t)ξ β (t) in the stationary state is found as The steady-state probability distribution P s (x) for the concentration vector x is then given by the multi-variate Gaussian In these calculations the only difference between the delay case and the well-studied non-delay case is the fact that the expression for S(ω) is more complicated, in particular it involves the Fourier transform of the delay kernel K(•).For the Brusselator model with delay we find A(ω) ξ(ω) = η(ω), where the matrix A(ω) reads (32) From this one has where † denotes the conjugate transposition and where It follows, for example, that Similar expressions can be found for S 22 (ω) or indeed for cross-correlation spectra.Eq. ( 35) can be numerically integrated (over ω) to find Ξ 11 and hence P s (x 1 ).

Quasi-cycles and stationary distributions
A sample trajectory for the Brusselator with delay kernel K(t) = δ(t − τ ) is shown in Fig. 2. The data is from numerical simulations using the modified next-reaction method (MNRM) [39].The MNRM is an algorithm which simulates the stochastic process exactly, for details see Appendix B. Model parameters are chosen such that the limiting deterministic model reaches a stable fixed point.The trajectory is shown in the steady state and, as seen in the figure, it fluctuates about the deterministic fixed point.The corresponding marginal distribution P s (x 1 ) in the steady state is shown in Fig. 3 for different choices of the system size Ω.Semi-analytical results for this distribution, calculated using Eq.(31), are in good agreement with the results of the numerical simulations, even for the relatively low value of Ω = 10.We stress that the LNA is only valid when effects of higher order than Ω −1 can be safely neglected [9].As expected, the distribution P s (x 1 ) becomes more sharply peaked around the deterministic fixed point when the system size is increased.
Another feature seen in Fig. 2 is a temporal structure to the fluctuations.The concentration variables appear to undergo noisy oscillations with an average period of approximately T = 2. Recall that simulations are carried out with parameter values for which the deterministic system has a stable fixed point.Hence the oscillations are an effect of the intrinsic stochasticity, they are noise-driven quasi-cycles which have been widely discussed for other systems in the literature, see e.g.[3].These cycles can be analysed further by means of the power spectrum S 11 (ω), similar to existing studies of delay systems [20,22].While Eq. ( 33) provides a general expression, applicable to arbitrary delay kernels, we carry out this analysis for a family of Γ-distributed delays with a  36), used for the simulations in Fig. 5.The average delay is fixed to τ = 2.
fixed average delay τ .Specifically we choose where Γ(L) = ∞ 0 du u L−1 e −u .The parameter L controls the shape of the distribution, as shown in Fig. 4. The kernel becomes increasingly sharply distributed about τ as L is increased.For L = 1 the kernel K(t) is an exponential distribution, and for L → ∞ it is the δ-distribution K(t) = δ(t − τ ) used in Figs. 2 and 3.For all L one has ∞ 0 dt tK(t) = τ .Results for the power spectrum, S 11 (ω), are shown in Fig. 5 for different choices of L. In all cases we find very good agreement between the theoretical predictions and simulation data.For L = 1 the power  36)) with τ = 2 and L = 1, 20, 100, ∞.Noise-induced oscillations are observed for large values of L, but not for the exponential delay kernel, L = 1.The remaining parameters are the same as in Fig. 3. Noisy lines show data from simulations using the MNRM with Ω = 100, smooth lines are from the theory, Eq. ( 35).spectrum of fluctuations in the concentration of particles of type X 1 decreases monotonically from its peak value at ω = 0.This is characteristic of a system which does not display noise-induced oscillations.For larger values of L, i.e. more sharply peaked delay distributions, an additional peak emerges at ω ≈ 3 corresponding to oscillations with period of approximately T = 2.1.Higher harmonics are seen as well.As L → ∞ this peak grows and becomes dominant, indicating coherent stochastic oscillations and confirming the observations of Fig. 2. The steady-state probability distribution P s (x 1 ) as a function of L can also be investigated.It turns out that the variation with L is relatively insignificant, so that we do not show results here, and limit ourselves to stating that the semi-analytical theory predicts the simulation outcome to a good accuracy.

Effects of the delay period
The effects of varying the delay period, τ , on the power spectrum of the delay Brusselator with fixed delay are shown in Figs. 6 and 7. When there is no delay (top-left panel of Fig. 6) the power spectrum of fluctuations in the conventional Brusselator is recovered [28].For the model parameters a, b and c chosen here, the system displays noise-induced cycles, as indicated by the peak at a non-zero frequency.As τ is increased this peak shrinks and disappears.New peaks appear located approximately at frequencies ω n = 2nπ/τ , with n = 1, 2, . . . .As τ is increased the peak located at ω n grows and moves to the left, c.f. the first peak in the top-middle and top-right panels of Fig. 6.We find empirically that the peak height of the n-th peak increases with τ so long as ω n 2.6.Once τ is sufficiently large so that ω n 2.6 the height of peak n decreases upon further increase of τ , see e.g.height of the first peak in the spectra for τ = 2, 3, 4, 5 in Fig. 6.Fig. 7 shows the value of ω for which S 11 (ω) is at its maximum.For small τ this maximum is attained at a small angular frequency, and then jumps to a significantly larger value of ω at τ ≈ 1.2.This occurs when the maximum of S 11 (ω) is attained at ω ≈ ω 1 = 2π/τ .Further increasing τ the maximum of the power spectrum remains at the first peak.The location of this peak, ω ≈ 2π/τ , gradually decreases along with peak height.At τ ≈ 3.5 the second peak becomes dominant, and the maximum of the spectrum is now found at ω ≈ ω 2 = 4π/τ .As τ is increased even further this process repeats.This leads to a discontinuous behaviour of the dominant frequency as a function of τ , as shown in Fig. 7.We stress that the peak locations are only approximately given by ω n = 2nπ/τ , careful inspection of the data shown in Fig. 7 reveals quantitative deviations.The behaviour we have described indicates that delay can have two distinct effects: for τ 1.2 the delay dampens the ability of the system to oscillate, whereas for τ 1.2 the delay induces new oscillations in the system.Their characteristic frequency shows a relatively complex dependence on the delay period.

Spiking behaviour
We now focus on a different parameter regime, one in which the Brusselator system with delay becomes an excitable system, and where it shows spiking behaviour [35].In this situation the deterministic system has a stable fixed point, and separated from it in phase space a deterministic 'excursion' orbit, taking the system along a looped trajectory with significant amplitude and then back to the fixed point.These excursions are triggered from starting points separated from the deterministic fixed point by a finite distance.In the face of noise (either intrinsic or extrinsic) excitable systems typically fluctuate about the deterministic fixed point, and if fluctuations take it across a threshold triggering an excursion, a spike will occur.Once the spike is complete the system spends time near the fixed point again, until a new excursion is triggered by noise.To demonstrate this phenomenon in the delay Brusselator system we show a sample trajectory in Fig. 8.The system mostly fluctuates with small amplitude about the a deterministic fixed point, but intermittent spikes are seen as well.This is a phenomenon seen in many models of nonlinear dynamics with and without delay, including e.g.models of plankton bloom [40] or the celebrated FitzHugh-Nagumo model of firing nervous cells [41].We will now use the CLE approach to study this phenomenon in more detail.We denote the time between spikes by σ, and measure the average time between spikes in simulations.Some care needs to be taken here to identify spikes in simulations.We identify the beginning of a spike as a point in time in which the concentration x 1 crosses a lower threshold from above, see Fig. 8, and the spike is taken to end when the trajectory of x 1 has not crossed either the upper or lower thresholds for a fixed time, t s = 2.In the example shown in Fig. 8 the deterministic fixed point corresponds to x * 1 = 1, the upper and lower thresholds are chosen as x l = 0.4 and x u = 1.75.The location of the two thresholds and the value of t s are in principle arbitrary, however the choices we make result in a satisfactory identification of spikes.Fig. 9 shows the mean time between spikes, σ , in the delay Brusselator.Exact simulations using the MNRM are in very good agreement with simulations of the CLE.As seen in the figure the typical time between spikes increases exponentially with Ω.In the regime of large Ω efficient simulations of the CLE can be advantageous.As discussed further in Appendix B the computing time required to perform simulations of the CLE is largely independent of the system size, whereas computational resources to run the MNRM increase quickly with Ω.The exponential behaviour of the time between spikes is in-line with observations that escape rates from locally stable fixed points in dynamical systems subject to noise scale exponentially with decreasing noise strength [42], and the fact that spikes are triggered by excursions of the dynamics from the fixed point crossing a specific threshold [41].In the deterministic limit, Ω → ∞ the time between spikes is infinite, there are no spikes in the deterministic system.At finite noise strength, the duration and amplitude of any given spike remain roughly constant as Ω is varied, the spikes are different from the stochastic quasi-cycles discussed above.Quasicycles are sustained by noise, if the stochasticity were to be removed or reduced in the stationary state the quasi-cycles would disappear (or their amplitude be reduced).The spikes in Fig. 8 are triggered by noise, but once a spike is triggered it completes irrespective of the stochasticity, fundamentally following a deterministic orbit.The noise does not affect the spikes' amplitude, but only the frequency with which they occur.The time between spiking phases is very sensitive to the delay τ , as seen in Fig. 9, where we show results for two values of τ .An increase in the delay leads to significantly larger average inter-spike time σ .It is important to stress that the spiking behaviour results from the combination of intrinsic noise and the nonlinear dynamics of the delay Brusselator system.Neither the deterministic limit, nor the LNA capture the spiking behaviour completely.The spiking orbits are present in the deterministic flow, but the frequency with which they occur cannot be obtained from a deterministic analysis alone.A linear analysis about the fixed point on the other hand neglects all nonlinearity, and hence the spiking orbits are not captured.A recent paper has investigated how the nonlinearity of the CLE can be studied analytically (for systems without delay), and this may provide a starting point for further analytical study of similar effects in delay systems [43].

V. CONCLUSIONS
In summary we have presented an intuitive heuristic derivation for the Gaussian approximation of discrete-particle dynamics with distributed delay.These approximations build on our earlier work [22], in which we have used a more formal approach based on functional integrals.Our main result is a chemical Langevin equation for systems with general delay distributions, ready to be applied for the efficient simulation of stochastic delay systems of a large, but finite size.The chemical Langevin equation can be reduced in a linear-noise approximation, which allows one to make further analytical progress, in particular with a view towards calculating stationary probability distributions and correlation properties of the stationary state.
We have applied our results to the example of a Brusselator system with delay.Comparison against simulations shows that the linear-noise approximation works well.We have characterised the spectra of noise-induced quasi-cycles, and we have shown how simulations of the chemical Langevin equation allow for an efficient computational characterisation of noise-triggered spiking behaviour for parameter choices in which the delay Brusselator becomes an excitable system.
We expect that the procedure for the systematic derivation of Gaussian approximations in systems with distributed delay will be of interest to a number of applications, for example in other chemical reaction systems or in the biological sciences.The approach we have presented here is intuitive and relatively easy to apply and to generalise.Model features that have not been captured so far are the uncertain completion of delay reactions or cases in which a delay reaction can result in several later outcomes, possibly depending on the state of the system at the scheduled completion time.Work is in progress to extend the formalism to such cases [44].
3. Generate an independent random number ρ i for all i, each drawn from a uniform distribution over (0, 1].Set P i = − ln(ρ i ).
4. For each i set ∆t i = (P i − T i )/(Ωr i (x)). is the first (i.e.lowest) entry in s i .We write i 0 for the corresponding reaction, i 0 = arg min i {∆t i , s i (1) − t}.If the minimum is attained at ∆t i0 the selected event corresponds to the initiation of a new reaction.If the minimum is at s i0 (t) − t, then the selected event is the completion of a delay reaction.
7a.If the event selected in 5. corresponds to the completion of a delay reaction then • Update the state of the system: x α ← x α + w i0,α /Ω for all α.
• Delete the first entry in s i0 .
7b.If the event selected in 5. corresponds to the initiation of a new reaction • Update the state of the system x α ← x α + v i0,α /Ω for all α.
• If reaction i 0 is a delay reaction: -Draw a random number τ from K i0 (•).
-Add the entry t + τ to the list s i0 .
-Sort s i0 so that the list remains in ascending order.
10. Recalculate the reaction rates r i (x).
11. Go to step 4, or exit if final time of the simulation is reached.This is a minor change to algorithm 7 presented in [39].The only difference relates to the maintenance of ordered lists for each reaction channel.For systems with fixed delay, as discussed in [39], updating these lists is straightforward, events are completed in the order they are initiated (they all have the same delay period).Extending the algorithm to cover distributed delay requires two steps: (i) drawing a random delay time from the appropriate delay kernel when a delay reaction is initiated, and (ii) sorting the list whenever a a new entry is added, to ensure the entries are always in ascending order.
For systems without delay the time between events scales as Ω −1 as the reaction rates are all proportional to Ω.The computing time to run the MNRM algorithm up to a fixed total time thus grows linearly in the system size Ω.For models with distributed delay the simulation time increases even faster with Ω as additional time is needed for the sorting of the list of delay events, the typical length of which itself increases with Ω.

Simulation of the chemical Langevin equation
The alternative approach is to simulate the CLE, Eq. ( 16).These equations are not an exact representation of the original microscopic model, they are an approximation for large, but finite systems.As such, simulations of the CLE will not be exact.On the other hand, simulation times of the CLE are independent of the system size, Ω.This scale only enters through the noise amplitude, and it does not affect the run time needed to integrate the CLE numerically.For large Ω the simulation of the CLE is not only accurate (in this limit it represents the original process faithfully), but also faster than simulation of the discrete-particle model.In practice the discretisation of the CLE is implemented using Eqs.( 14) and (15).In the case of Markovian systems any one Gaussian random variable only enters at one single time, see Eq. ( 8), and they can hence be removed from the computer memory once they have been used.For delay systems it is necessary to store the ζ τ i,t , generated at time t, for later use at t + τ .

FIG. 1 :
FIG.1:(Colour on-line) Diagram indicating the reaction events contributing to the change of particle numbers at a given time t.Particle numbers are changed due to instantaneous effects of reactions triggering at time t, as indicated by ki,t.These may then have delayed effects at future times, indicated by the filled arrow.Other changes of particle numbers at t are due to delayed effects of reactions triggered at earlier times, as indicated by the arrows on the left.For simplicity we have set ∆ = 1.

2 FIG. 2 :
FIG. 2: (Colour on-line) Sample trajectory of the delay Brusselator with fixed delay.Parameters are Ω = 100, a = 2, b = 2.9, c = 1, τ = 2. Oscillations with period of approximately 2 time units can be identified visually.These observations are validated by the results shown in Fig. 5.The system has been initialised at the deterministic fixed point x1 = 2, x2 = 1.45.

4 FIG. 3 :FIG. 4 :
FIG.3:(Colour on-line) Marginal steady-state probability distribution P s (x1) for the Brusselator with fixed delay and for different system sizes.Parameters are a = 2, b = 2.9, c = 1, and τ = 2. Results obtained from the LNA (black lines) provide a good approximation to the exact simulations of the process using the MNRM, shown using coloured (grey) lines.

FIG. 5 :
FIG.5:(Colour on-line) Power spectra S11(ω) of the delay Brusselator for Γ-distributed delay kernels (see Eq. (36)) with τ = 2 and L = 1, 20, 100, ∞.Noise-induced oscillations are observed for large values of L, but not for the exponential delay kernel, L = 1.The remaining parameters are the same as in Fig.3.Noisy lines show data from simulations using the MNRM with Ω = 100, smooth lines are from the theory, Eq. (35).

5 FIG. 6 :FIG. 7 :
FIG.6:(Colour on-line) Power spectra S11(ω) for the delay Brusselator, calculated from Eq. (35) with fixed delay τ .The remaining parameters are the same as in Fig.3.Exact simulations with the MNRM have also been performed to confirm the spectra (data not shown).

FIG. 9 :
FIG.9:(Colour on-line) Mean time between spikes, σ , as a function of Ω for the delay Brusselator.Blue squares correspond to MNRM simulations, ran up to a final time t f = 81920.Red circles correspond to simulations of the CLE, run up to time 10t f with a time step ∆t = 10 −5 .The black line shows a fit of an exponential function to the data from the MNRM simulations, with which the CLE simulations can be seen to agree.Parameters used are the same as Fig.8, apart from τ and Ω which are as indicated.

5 .
Set ∆ = min i {∆t i , s i (1) − t}, where s i (1) To derive the Gaussian approximation one replaces the Poissonian random variables by Gaussian noise with these first and second moments, i.e. k i,t → R i (x t )∆+ √ ∆Ωζ i,t with ζ i,t ζ j,t = δ i,j δ t,t r i (x t ),