Generating functionals and Gaussian approximations for interruptible delay reactions

We develop a generating functional description of the dynamics of non-Markovian individual-based systems, in which delay reactions can be terminated before completion. This generalises previous work in which a path-integral approach was applied to dynamics in which delay reactions complete with certainty. We construct a more widely applicable theory, and from it we derive Gaussian approximations of the dynamics, valid in the limit of large, but finite population sizes. As an application of our theory we study predator-prey models with delay dynamics due to gestation or lag periods to reach the reproductive age. In particular we focus on the effects of delay on noise-induced cycles.


I. INTRODUCTION
In recent years there has been a growth in the understanding of the importance of intrinsic noise in complex systems composed of a finite number of interacting constituents. It has been recognised that there are visible, macroscopic effects due to the stochasticity inherent in the interactions in a variety of systems, including for example metabolic pathways [1], gene regulatory systems [2][3][4][5], predator-prey dynamics [6,7] or models of disease spread [8,9]. This inherent stochasticity is referred to as 'intrinsic noise' or 'demographic stochasticity' [10]. Intrinsic noise can give rise to phenomena such as cyclic dynamics [7,11], patterns and waves [12][13][14][15][16], and extinction events [17][18][19][20][21]. These phenomena are not captured by more traditional deterministic modelling approaches, instead they are purely noise-induced. Deterministic models are built on ordinary or partial differential equations. These equations are formally only valid in the limit of an infinite system, that is the number of interacting constituents is so large that stochastic effects play no role [22]. To take into acount the intrinsic stochasticity of the interactions a full probabilistic description is required. Widely used modelling approaches are drawn from the theory of stochastic processes, most notably the master equation [22,23], describing the time evolution of the probability distribution over the space of states. Formulating the master equation approach relies on the Markov property of the underlying dynamics: transition rates from one state to another must only depend on the target state and the present * tsbrett@umich.edu † tobias.galla@manchester.ac.uk state of the system, but not on the path the dynamics has taken to arrive at the current state. This implies that the system has no memory of previous interactions, all effects of an interaction must be realised instantaneously. While this is a reasonable assumption for many processes in physics, this is typically not the case in biological models. For example birth in a predator-prey system occurs after a period of gestation, transcriptional and translational delays are relevant in gene regulatory systems [24], and recovery occurs a certain period of time after infection in models of disease spread [25,26]. These processes can lead to situations in which effects of an initial interaction (e.g. impregnation, initiation of transcription, infection) materialise with a significant delay. Such dynamics are then no longer Markovian as the change of the state of the system at any one time t may depend on what processes were set in motion in the past and which complete at t. Of course the modelling as a delay system can often be avoided through the construction of higher-dimensional models with intermediate states (e.g. staged models in epidemics [26]), but mathematically it is often convenient to use delay models, as these are relatively easy to set up for arbitrary delay distributions.
Traditionally delay dynamics have been investigated based on deterministic approaches [24,[27][28][29][30][31]. These are subject to the limitations outlined above in that they do not capture the effects of stochasticity. It is only recently that analytical (and also numerical) techniques for individual-based models subject to both intrinsic noise and delay have been developed [32][33][34][35][36][37][38][39][40][41][42]. In previous work we have derived Gaussian approximations of such dynamics, and we have shown that these can capture the effects of delay reactions [40,41]. This analysis is conveniently carried out in terms of generating functionals [43,44], arXiv:1505.05100v1 [cond-mat.stat-mech] 19 May 2015 and the equivalent of the linear-noise approximation (LNA) for delay systems can be derived. An alternative, more informal method to derive the same Gaussian approximations (by-passing the exact description) was highlighted in [41].
In these existing approaches delay reactions fire at an initial time with rates determined by the state of the system at that time. The reaction may then have an immediate effect on the composition of the system, and subsequently a second effect materialises at a later time (after the delay). The delay is either fixed, or it can be drawn from an underlying distribution each time a delay reaction is initiated. Existing work is often restricted to what we will refer to as 'definitive' completion of the delay reaction. Once a delay reaction is initiated (it 'fires') the delayed effect will occur, no matter what the trajectory of the system is between the time the reaction fires and the designated completion time. This is an obvious limitation of the modelling approach. Descriptions in which delay reactions may fail to complete depending on events that occur between initiation and designated completion provide a much more realistic description of many real-world processes. The gestation period in a predator-prey model presents perhaps the most intuitive example of delay reactions which can be interrupted. A pregnancy is initiated at a given time and the birth event will occur at a designated later time. However, birth is not certain, the mother might die during pregnancy. Similarly, in a model of disease spread, an individual may get infected with the disease and would then be scheduled to recover at a later time. However, the individual may die in the interim, or be removed from the system in some other way, so that the completion event (recovery) may not occur. The rate with which such removal occurs may well depend on the state of the system at the time of removal.
The purpose of the present paper is to extend existing approaches to the modelling of stochastic dynamics with delay to the case in which delay reactions may not complete. We will refer to such reactions as 'interruptible delay reactions'. In this setting delay reactions can be of several types: (i) Delay reactions which cannot be interrupted. This is the case considered in [40]; (ii) Delay reactions which can be interrupted, but the probability of interruption does not depend on the state of the system. An example is death of the mother for internal reasons. Another example was also considered in [40] in the context of the susceptible-infective-recovered (SIR) model with birth and death, with constant death rate; (iii) Delay reactions which can be interrupted and the probability of interruption depends on the state of the system at the time of interruption. Continuing the example of a pregnant prey individual, the rate of predation depends on the number of predators present in the system throughout the pregnancy period. This case is not covered by [40]. In this paper we develop a systematic Gaussian approximation to models with interruptible delay reactions. Specifically we extend the generating functional method used in [40] to include interruption effects so that all three cases above are covered. As an application we consider the SIR model of disease spread and two variants of a predator-prey model, one with a delay period due to pregnancy and one with a delay period due to the maturation of juveniles.

II.1. Model definition
We consider a finite population of interacting constituents, each of which is of one of M different types, labelled α = 1, . . . , M . We will refer to the constituents as 'individuals' in the following. The state of the system at any one time is then described by an M -dimensional vector, n(t) = [n 1 (t), . . . , n M (t)], where the non-negative integer n α (t) indicates the number of individuals of type α at time t. We assume a continuous-time evolution. The system is well-mixed and two individuals of the same type are indistinguishable. The individuals interact through a set of R reactions, we will label these i = 1, . . . , R. Each possible reaction i is associated with an initiation rate T i (n), indicating the rate with which reactions of type i fire if the system is in state n. When such a reaction fires an instantaneous change of the state of the system n(t) occurs at the time of firing, described by the vector v i . That is to say the state of the system changes from n to n + v i . Subsequently these reactions may also have a delayed effect. This is implemented as follows: if a delay reaction fires at time t an instantaneous change of the state of the system occurs as described above. In addition to this, a delay time τ is drawn from an underlying delay distribution, K i (τ ). This distribution is specific to the reaction triggered, as indicated by the subscript i. After the delay time has elapsed a second change in the state of the population may occur. This happens at time t + τ , and we denote the delayed effect of the reaction by w i . The key element of the dynamics that we add in the present work is the possibility that a delay reaction may not complete. We assume that a delay reaction of type i, triggered at time t and due to complete at t + τ can be interrupted at any time between t and t + τ . The rate with which this occurs is f i [n(t )], where t < t < t+τ . The termination rate may thus depend on the state of the system n(t ). If this happens the delayed effect at time t + τ does not occur, instead we assume that the state of the system changes by u i at time t . For the time being we will only consider cases in which there is only one way in which each delay reaction can be interrupted. A generalisation to models in which delay reactions can be interrupted in multiple different ways is relatively straightforward, though slightly more cumbersome, see Appendix B for details. As is commonly done in the modelling of interacting particle systems, we will assume that all reaction rates T i (n) scale with a parameter Ω -that is to say T i (n) = O(Ω). This parameter can be seen as setting the scale of the size of the system (the scale of the total number of individuals in the system, or equivalently the volume of the system), and the scaling of the rates reflects the fact that the total (average) number of reactions in the system per unit time is proportional to its size. For later purposes it is convenient to introduce the quantities x α (t) = n α (t)/Ω. We will refer to these as the 'concentrations' of particles of type α. We also introduce the intensive rates

II.2. Discrete-time dynamics
As in [40] we will proceed by discretising time into steps of duration ∆. The continuum limit will be restored at the end. We will write x α,t for the concentration of individuals of type α at time step t. In the discretised model, we will assume that all reaction rates remain constant between t and t + ∆. If the concentration vector is x t at time t then the number of newly triggered reactions of type i during this time step is a Poissonian random variable k i,t with parameter T i (n)∆. In the absence of delay reactions the dynamics of the discrete-time model would hence read In models with delay reactions we have to extend this expression to take into account (i) delay reactions completing at the designated time and (ii) delay reactions which terminate before the designated completion time. We then have In this expression k τ i,t is the number of reactions of type i that fire at time t and which have a delay period τ . The quantity m τ i,t−τ is the number of reactions of type i that fired at time t−τ and successfully complete at time t. Finally, the quantity τ,s i,t−s is the number of reaction of type i that fired at time t − s with a designated delay period τ (i.e. scheduled for completion at t − s + τ ), but which are interrupted at time t (0 < s < τ ). We note that non-delay reactions are included in this descriptions, they would simply have w i = u i = 0. Eq. (2) defines the conditional probability where the notation {x, k, , m} t ≤t indicates variables with indices t ≤ t.
We have already described the Poissonian nature of the variables k i,t = τ k τ i,t , but it remains to define in more detail how the k τ i,t , m τ i,t−τ and τ,s i,t−s are chosen. This is explained in the following section.

II.3. Statistics of reaction numbers
We first discuss the statistics of the variables k τ i,t , indicating the number of reactions of type i triggered in the discrete time model at time step t and with completion due at time t + τ . As discussed in [40] each k τ i,t is a Poissonian random variable with parameter ∆ 2 K i (τ )Ωr i (x t ). This is not affected by possible interruptions of the delay reactions. At this point it is useful to recall the definition r i (x t ) = T i (x t )/Ω. Broadly speaking Ωr i (x t ) reflects the rate with which reactions of type i are initiated (with any delay). Once a reaction is initiated a delay period is drawn independently from a distribution K i (τ ). In the discrete-time model this is reflected in the factor ∆K i (τ ). Ultimately we will be interested in taking the continuous-time limit for the dynamics. Anticipating this, we focus on the case of small ∆, which simplifies the problem significantly. The probability distribution of k τ i,t is of the form Eq. (4) indicates that the probability to observe two or more initiation events of reactions of type i and with delay τ in any one time interval ∆ is at least of order ∆ 4 . In the limit of small ∆ we can therefore restrict ourselves to k τ i,t ∈ {0, 1}. If k τ i,t = 0 then it is clear that τ,s i,t = 0 for all s and also m τ i,t = 0, i.e. we have If k τ i,t = 1 then only one out of the { τ,s i,t } 0<s<τ and m τ i,t can be non-zero. If k τ i,t = 1 and τ,σ i,t = 1 then it follows that τ,s i,t = 0 for all s = σ, and m τ i,t = 0. I.e., we have As explained above, interruption happens with probability ∆ × f i (x t+s ) in the next time interval ∆ if the system is in state x t+s , and so Finally, if the reaction has not been interrupted by time t + τ then the reaction always completes,

III. GENERATING FUNCTIONAL
In discrete time the generating function is where the average is performed over all possible paths. We have introduced a source term ψ, derivatives with respect to it generate correlation functions [43,44] The main task in the calculation that follows is to compute the path-average in the above expression. This entails integrating/summing over all random variables, where P(x, k, , m) is the joint probability distribution of all x t,α , k τ i,t , τ,s i,t , and m τ i,t -i.e. the probability of a path. Our objective is to find an expression for the generating functional after performing the sums over all k, , and m and after subsequently taking the continuous-time limit ∆ → 0. The algebra involved is somewhat lengthy, and so we do not give full details here, they are relegated to Appendix A. The calculation consists of carrying out the following main steps: (i) The path probability P(x, k, , m) is expressed as a product of conditional probabilities using Eqs. their statistical weight. Knowing the weight of each combination enables us to average over the k, , and m; (iv) Finally, the continuous-time limit is restored. The resulting generating functional is with the action There are two contributions to the action from each type of reaction: (i) a contribution from the reactions which are not interrupted, (13) and (ii) a contribution from the reactions which are interrupted, These contributions are both functionals of the path of x between the time at which reaction first fires and the time it either completes successfully or is interrupted. We notice the factor e − τ 0 dσfi[x(t−σ)] in Eq. (13). This exponential is the probability that a delay reaction of type i triggered at t − τ and with completion time t, reaches completion, given a path x between t − τ and t. Similarly, the quantity (14)] indicates the probability that a delay reaction of type i triggered at time t − s is interrupted in the time interval [t, +dt), given a path x between t − s and t. If all delay reactions of type i complete with certainty (i.e. in absence of interruptions) one has f i (x) = 0. This implies R (2) i = 0 (see Eq. (14)). The quantity R (1) i (t) in Eq. (13) reduces to the corresponding object in [40].

IV. APPLICATION TO THE SIR MODEL WITH BIRTH AND DEATH
The SIR model with birth and death can be studied using the approach developed in Sec. II. This model describes the dynamics of an infectious disease in a population of individuals. Each individual can be in one of three states: susceptible, infectious or recovered. We consider an infection dynamics with delayed recovery. Upon contact with an infectious individual a susceptible member of the population may become infectious (with rate β), and then they recover at a time τ after infection. The delay time τ is drawn from a recovery-time distribution K(τ ). We write this process as follows The first arrow indicates the effect the reaction has at the time it is initiated. The double arrow indicates the delayed effect τ units of time after the reaction is triggered. In addition to the recovery process, any individual may die, this occurs at constant rate µ, independent of the infection status of the individual. In order to keep the population size, N , constant any dying individual is immediately replaced by an individual of type S. This leads to the following additional reactions Crucially, an individual may die (and be replaced by an S) after becoming infected, but before the scheduled recovery time. This represents a delay reaction with possible interruption in the formalism described in the previous sections. In this introductory example, however, the rate of interruption does not depend on the state of the system and so this model can be studied using the simpler method presented [40]. There we mapped the above reaction scheme onto the model with The first line in Eq. (17) describes the standard death and replacement reaction R µ −→ S. The second and third lines represent two reaction sequences which may be triggered when an infection event occurs. Each sequence starts with an infection event (occurring with rate βn I n S /N ). With probability χ the newly infected individual is scheduled for recovery at a later time; the time-to-recovery, τ , is drawn from the distributionK(τ ). With complementary probability 1 − χ the newly infected individual is scheduled for death (and replacement by an individual of type S). The time-to-replacement, s, is drawn from the distributionQ(s). This mapping is possible because the interruption rate, µ, is independent of the state of the system. Both distributions,K andQ, are normalised. The reaction scheme in Eq. (17) and the distributions in Eq. (18) were formulated in [40], and from this scheme the generating functional of the process was derived. Now that we have put the more general formalism of the previous sections in place we can recognise these existing results as a special case. Inserting the model specifications into the general formalism of the previous sections the action in Eq. (12) is indeed found as where we have used the definitions of χ,K(τ ), and Q(s) as given in Eqs. (18). The total size of the population, N = n S + n I + n R is constant in time, so there are only two independent degrees of freedom. In formulating Eq. (19) we have therefore eliminated species R via n R = N − n S − n I (equivalently It is important to stress that we have used Ω = N for simplicity, given the constant size of the population. The result in Eq. (19), derived from the general formalism of the previous sections, is the same action as the one found in [40] using the mapping of Eq. (17). The exact generating functional can be approximated by Taylor expanding the action in powers of the inverse system size. To leading order in Ω −1 the general action Eq. (12) reduces to and This action is linear in p (by construction), and it describes the deterministic dynamicṡ with the drift terms and featuring due to the delay. The deterministic approximation is valid for large (formally infinite) system sizes, and neglects all stochastic effects. We have denoted the dynamical variables in the deterministic limit by x ∞ .

V.2. Linear-noise approximation
The linear-noise approximation is used to study fluctuations about the solution to the deterministic equations of motion, Eq. (23). We write x = x ∞ +Ω −1/2 ξ, and formulate the generating functional in terms of the variable ξ. In doing this we re-scale the conjugate variables via p → √ Ωq. The resulting generating functional is Differentiating Eq. (26) with respect to ψ(t) still gives the moments and correlation functions of x(t). To find a generating functional for ξ we divide Eq. (26) by the factor e − dt ψ(t)·x ∞ (t) , which is independent of ξ and q, and rescale the source term, ψ → √ Ωφ. The generating functional for ξ is The moments of ξ can now be found by differentiating The terms in the action are again expanded in powers of Ω −1/2 , and the expansion is curtailed after subleading order (i.e. keeping terms O(Ω 0 ) and above). In order to illustrate the procedure we will consider the case where the deterministic dynamics are at a stable fixed point, x ∞ (t) = x * . This is not necessary for the linear-noise approximation, however it will allow us to keep the amount of algebra under control. The action is In this expressionẋ ∞ α stands for the right-hand side of Eq. (23). The action is now quadratic in q and ξ, and it describes a process with additive Gaussian coloured noise [44,45]. This process is of the forṁ (29) with noise correlator η α (t)η β (t ) = B α,β [x * ](t, t ). The explicit expressions for A α,β [x * ](t, t ) and B α,β [x * ](t, t ) are somewhat lengthy, and before presenting them it is helpful to make a few simplifying definitions.

V.3. Further simplications
To simplify the above expressions further we define the distributions for τ, s > 0. We setK i (τ ) = 0 for τ ≤ 0, andQ i (s) = 0 for s ≤ 0, and we introduce Within the linear-noise approximation and assuming that the deterministic dynamics are at the fixed point x * this is the probability that a delay reaction of type i, once triggered, completes. The functionK i (τ ) is the conditional time to completion. The quantity 1 − χ i (x * ), on the other hand is the probability that the delay reaction is interrupted before completion, andQ i (τ ) describes the conditional distribution of times to interruption. Both distributionsK i (τ ) andQ i (τ ) are normalised. One notices that the expressions in Eq. (18) are special cases of Eq. (30), with f i (x * ) ≡ µ. We stress though that our analysis below will explicitly account for (linear) fluctuations of the interruption rate as the state of the system varies in time. These fluctuations of the interruption rate are obviously not present in the system with constant f i (x) = µ, considered in [40].
Additionally it is helpful to make the definition This quantity can be interpreted as the weighted delay effects (due to either successful completion or interruption) on species α of reactions of type i a time period τ after such a reaction was triggered. Using these definitions we can write A α,β [x * ](t, t ) concisely as Each of the four terms in the curly brackets has a clear interpretation. As the system fluctuates around the fixed point, the number of reactions triggered, interrupted and completed will fluctuate as well. The first term is the contribution toξ α (t) of initial effects of reactions (at the time of triggering) due to fluctuations. The second term is the contribution of interruption effects due to fluctuations at the time of interruption. The third term is the contribution of delayed effects (either successful completion or interruption) due to fluctuations between the reaction firing and the delayed effects occurring. Finally, the fourth term is the contribution of delayed effects due to fluctuations at the time the reaction first fired.
If the rate of interruption, f i , is independent of the state of the system for all i then only the first and last terms in Eq. (32) are non-zero. We now turn to the correlation B α,β [x * ](t, t ) of the Gaussian noise variables in the linear-noise approximation. Given that delay reactions can have effects on the composition of the population at multiple times, this noise is not white. Instead one finds As L i,α (t) = 0 for t ≤ 0 only one of the three terms in Eq. (33) is non-zero for any pair of times t and t and any given reaction i.

V.4. Spectrum of fluctuations
The linear-noise approximation can be used to calculate the Fourier spectrum of fluctuations about the deterministic fixed point [7]. This is particularly useful to study noise-induced quasi cycles, as we will explain in further detail below. Fourier transforming Eq. (29), gives This equation is linear in ξ, and it can be written as η(ω) = M(ω) ξ(ω), with a suitable matrix M(ω). The Fourier transform of Eq.
The spectrum of fluctuations about the deterministic fixed point is then characterised by the matrix S(ω) = ξ(ω) ξ † (ω) and it can be found from (see [46])

VI.1. Model definitions
As an application of the above formalism we will now consider stochastic effects in a predator-prey model with different types of delay. Specifically we write X to denote prey-individuals, and Y for predators. We focus on the dynamics governed by the following reactions We write n X and n Y for the number of individuals of each type, and x = n X /Ω, y = n Y /Ω. As before Ω is a parameter, setting the scale of the population size. The first reaction describes reproduction of prey with a logistic birth rate, dependent on the concentration x. The constant k represents a carrying capacity, more precisely, the system can contain a most kΩ individuals of the prey type. The logistic birth rate for prey distinguishes this model from other stochastic predator-prey models [7]. Within our stylised approach we assume that predation events result in the birth of a predator (second reaction), the third reaction finally describes a death process for predators. This is obviously a minimalist model, but it is in line with previous stylised modelling approaches for birthdeath processes [47], and it serves as a helpful testbed. There are various ways in which delay processes can be introduced into this model. One is to include gestation periods, in which a prey enters a pregnant state before giving birth. After the gestation period the pregnant individual returns to the regular (nonpregnant) state and a new prey individual is created. Another possibility is to assume that newly born prey individuals are initially in a 'juvenile' state and that they cannot immediately reproduce. After a maturation period they become full adult prey individuals and acquire the ability to reproduce. We will refer to the first modification as the 'gestation model' and the second modification as the 'juvenile model'. In both models there is an additional intermediary class of individuals, denoted X -pregnant prey in the gestation model and juvenile prey in the juvenile model. The class X corresponds to non-pregnant prey in the gestation model and adult prey in the juvenile model. We write n X for the number of individuals of type X and and x = n X /Ω. Of course one could introduce a similar state in the context of predators, but in-line with the above stylised approach we keep the complexity of the model to a minimum. For the gestation model the reactions are where we have introduced x tot = x + x , the total concentration of prey. In the first reaction, a nonpregnant individual becomes pregnant with rate b(1− x tot /k). They then give birth after a delay drawn from K(τ ) to a non-pregnant individual and return to the non-pregnant state. For the juvenile model we use In this model an adult prey gives birth to a juvenile individual with rate b(1−x tot /k). After a delay drawn from K(τ ) the juvenile becomes an adult individual. Additionally, in both models, the intermediary individuals, X , can be predated upon, via It is this predation reaction that leads to the possibility of interrupting a delay reaction before its scheduled completion. If a pregnant individual is eliminated in a predation event, it will obviously not give birth at the designated time, and similarly if a juvenile individual is removed during predation it will not reach its reproductive age. Crucially, the rate with which such events happen depends on the state of the system (e.g. on the number of predators in the population at that time). The methods of [40] are hence not applicable, and we will instead use the extended and more general formalism developed earlier in this paper. Both models have the same reaction rates. The only difference between the two models lies in the changes of n X and n X when a delay reaction triggers and completes, i.e. in the numerical values of the stoichiometric coefficients v 1,x and w 1,x (we label the delay reaction by i = 1 in both models). It is therefore convenient to carry out the analysis for both models together, keeping v 1,x and w 1,x general. After performing the calculations we will substitute in for these parameters and compare the two models.
Simulations of the models can be performed using the Modified Next Reaction Method (MNRM) [34,35]. The interruption reaction shown in Eq. (40) fires like a conventional reaction in the MNRM algorithm with rate is the number of 'active' delay reactions at time t, i.e. reactions of type i = 1 which have fired but for which the delayed effects have not yet occurred. When the interruption reaction occurs an element of the list of queued delay reactions is selected at random with uniform weights, and removed from the list. The effects of the interruption reaction are then applied to the state of the system according to u 1,α . Aside from this the MNRM is unchanged.

VI.2. Deterministic dynamics
Both models describe three species, X, X and Y . The primary reactions are listed in Eq. (38) and (39) respectively. These are labelled i = 1, 2, 3 from top to bottom. The reactions rates for both systems are with h(x tot ) = b(1 − x tot /k). The first reaction, i = 1, is the only delay reaction. This reaction can be interrupted due to predation on the intermediaries, with rate f 1 [x(t)] = py(t) per intermediary (i.e. any instance of the delay reaction which is active at time t is subject to interruption with rate py(t)). The delay distribution is K(τ ). The stoichiometric coefficients, v i,α , describing changes to the system when reactions trigger can be summarised as follows, where the rows each stand for one reaction (i = 1, 2, 3) and the columns represent the three types of individuals (X, X and Y ). The only non-zero stochiometric coefficients for interruption events [see Eq. The deterministic equations of motion are found aṡ In principle we can also write an equation forẋ ∞ , however we know from the reactions that the number of X at any one time t is equal to the total number of X created up to that point (through birth) and which have not been removed through maturation or predation. We find

VI.3. Fixed point analysis
At the fixed point the left-hand sides of Eqs. (43) and (44) are zero, the concentration of intermediary individuals can be found from Eq. (45). We obtain where the stars indicate the fixed-point values of the relevant variables. We have also introduced the quantity χ(y * ) = ∞ 0 dτ e −py * τ K(τ ), representing the probability for a delay reaction to reach completion when the system is at the fixed point. For the gestation model this is the fraction of pregnant individuals which successfully give birth, for the juvenile model this is the fraction of juveniles which successfully mature and reach the reproductive age. There are three solutions to the above fixed point equations: all extinct (x = x tot = y = x = 0), only prey (x = x tot = k, y = x = 0), and coexistence. We will focus on the coexistence fixed point, and we will simply refer to it as the fixed point. We see that x * tot = d/p for both models and for any delay distribution.  (48), and find with φ = w x χ(y * ) + v x + 1 − χ(y * ). We will now restrict the further discussion to the case of constant delayτ , i.e. K(τ ) = δ(τ −τ ). In this case we can proceed with the analysis, and find χ(y * ) = e −py * τ . Eqs. (49)-(51) are a transcendental set of equations for x * , x * and y * , which cannot easily be simplified any further. However we can reparametrise the model, and treat χ as a parameter in Eqs. (49)-(51). The delay period,τ , is then a function of the model parameters throughτ = − ln(χ)/(py * ). After substituting in for y * we find If χ = 1 then there is no delay (τ = 0) in either model. Furthermore χ decreases with increasing delayτ . In the gestation model χ → 0.5 asτ → ∞, whereas χ → 0 in the juvenile model, as shown in Fig. 1. We stress that we have not formally established stability of fixed point. However for the parameters used in Fig. 1 the fixed point has numerically been seen to be stable up to at leastτ = 5.
In Fig. 1 we compare these theoretical predictions against data from simulations of the microscopic dynamics in finite populations (Ω = 1000). As seen in the figure the above deterministic analysis is in good agreement with stochastic simulations of the microscopic process in a finite system, Ω = 1000. The concentration of predator and prey individuals x * and y * , determined by Eqs. (49)-(51) decrease with the delay period in both models, whereas the concentration of intermediary individuals, x * , naturally increases when the delay period becomes longer. This is summarised in Fig. 2, obtained as a parametric plot of Eqs. (49)-(51) and Eq. (52). Comparison with exact stochastic simulations at finite Ω again shows good agreement. We notice that the concentration of predators is much more sensitive to the choice of model than the concentration of prey for non-zero delay.

VI.4. Linear-noise approximation
The above deterministic analysis is valid only in the limit of infinite populations. We next study the effects of noise induced by the stochastic dynamics when the number of individuals in the system is finite. Representative trajectories generated from stochastic simulations in Fig. 3 show that the effects of noise can The total prey concentration, x * tot is not affectected by the delay (black line), and is the same for both models. Symbols are from simulations using the MNRM averaged over 100 realisations with Ω = 1000. be quite profound. As seen in numerous Markovian systems [7,9,15,16,48] noise can generate sustained oscillations in parameter regimes in which the purely deterministic system approaches a stable fixed point. This effect has also been observed in stochastic dynamics with fixed and distributed delay [32,[39][40][41], the time series shown in Fig. 3 show that this phenomenon extends to delay dynamics with uncertain completion of delay reactions. The oscillations can be characterised by the power spectra of fluctuations about the deterministic fixed point, more precisely by the diagonal elements od S(ω) = ξ(ω) ξ † (ω) , introduced before Eq. (36). Fig. 4 shows how the power spectrum responds to an increase in the delay period. In both models the effect of delay is to increase the period of the quasicycles, as can be seen by a shift of the peak towards lower frequencies with increasing delay. Delay also increases the amplitude of the oscillations. For the gestation model the effect is more pronounced, the amplitude is much more sensitive to the delay. To calculate the power spectrum of fluctuations analytically we start from the general expression for the Fourier transform of the Langevin equation found in the linear-noise approximation, Eq. (34). For the predator-prey models we have where L α (ω) = ∞ 0 dτ e −iωτ L α (τ ) and G α (ω) =  where c = x * h(x * tot ). The remaining elements follow from the Hermitian property, B(ω) = B † (ω). The spectrum matrix S(ω) can be found by inserting Eqs. (53)-(56) into Eq. (36). Comparison of these theoretical predictions against data from simulation shows good agreement, see Fig. 4, and confirms the viablity of the linear-noise approximation for both models. If we return back to the definition of the models, the delay reaction for the gestation model is whereas for the juvenile model it is The total effect of the delay reaction (i.e. of the initial effects together with the delayed effects) is the same for both models: the number of individuals of type X increases by one. The difference between the models is in the intermediary changes.
To probe the effects of delay, let us compare the outcome of the delay models (τ > 0) with the corresponding Markovian models (obtained in the limit τ → 0). In this limit the reactions in Eq. (57) and Eq. (58) both reduce to The study of the fixed point behaviour and of the power spectrum of the models with delay shows that the location of the fixed point and the peak of the power spectrum at any fixedτ > 0 are closer to the corresponding quantities atτ = 0 in the juvenile model than in the gestation model. This indicates that it is not only the duration of the delay which determines whether or not it can be neglected, but also the details of the delay reaction. We note that the survival probability at the fixed point (shown in Fig. 1) is more sensitive to the delay in the juvenile model than in the gestation model.

VII. CONCLUSIONS
The main focus of this paper has been to extend the generating-functional approach for stochastic interacting particle systems to dynamics in which delay reactions can be interrupted. In particular we consider cases in which the interruption rate depends on the state of the system at the time of interruption. The probability with which the delay effects ultimately occur is then a functional of the path taken by the system during the delay period. We successfully set up a path-integral description of such systems. The action of the generating functional for delay reactions with interruption has a clear relationship with the actions found previously in models in which delay reaction complete with certainty. The example in Sec. IV shows how our extended theory includes results previously derived in simpler cases where the interruption rate is constant and independent of the state of the system at the time of interruption. The generating functional provides a starting point for further analytical calculations, specifically we derive Gaussian and linear-noise approximations from it to characterise noise-induced effects in finite populations. We apply these techniques to an example inspired by predator-prey dynamics. We study two variants of the model, one with gestation delay and one with maturation delay. The deterministic approximation can be used to give a good estimate of the probability that any given delay reaction completes. In the gestation model this describes the probability that a pregnant individual successfully gives birth, in the juvenile model it is the probability that a newborn individual successfully reaches adulthood. Starting from the linear-noise approximation we show how the formalism can be used to analytically study persistent noisy cycles, induced by demographic stochasticity. Finally, the overall effect of delay reactions is identical in the gestation and in the juvenile predatorprey model, that is to say the changes of delay reactions at the time they trigger taken together with the effects at completion are the same in both models. Our analysis shows that different quantities (e.g. survival probability, fixed point location, peak of the power spectrum) have different sensitivities to the delay in the two models. For both the deterministic fixed point and the characteristic frequency of noiseinduced oscillations, the duration of the delay period has a much stronger effect on the outcome in the gestation model than in the juvenile model. Conversely, the probability that the delayed effects successfully occur decreases much faster with increasing delay in the juvenile model than in the gestation model. These observations have implications for the construction of population models. Whether a delay can safely be neglected depends not just on the duration of the delay period itself but also on the details of the changes to the population at the start and end of the delay period. This indicates that care must be taken when approximating models with delay or intermediate processes by effective Markovian dynamics. In summary, our previous work [40] together with the present paper provides a systematic and coherent approach to studying stochastic effects in a wide range of delay systems, including those with distributed delay, delay reactions that can fail to complete with path-dependent rates and systems with multiple possible interruption channels. The generating functional is the natural mathematical entity which with to address these systems, and is the analog of the master equation in the context of Markovian dynamics. We have here only applied this method to a small set of examples, inspired by ecological dynamics, but we anticipate that the formalism that is now in place can be of value in the context of many other applications, including those in epidemiology, metabolic dynamics, gene expression and other fields which delayed dynamics.

ACKNOWLEDGMENTS
TB would like to thank the EPSRC (UK) for support.
Appendix A: The generating functional for interruptible delay

Conditional probabilities
To derive the generating functional we focus on the model in discretised time. The continuous-time limit is taken at the end of the calculation. We adopt the following notation: k t denotes all k τ i,t for which t = t, t denotes all τ,s i,t for which t + s = t, and m t denotes all m τ i,t for which t + τ = t. Using this notation the equation of motion for x t,α can be expressed in the form x t+∆,α = x t,α + g α (k t , t , m t ). The path probability P (x, k, , m) can be written as the product For a particular t, all of the variables contained in {k t , t , m t } are independent of each other (conditioned on the history of the system up to time t), so their joint probability distribution factorises All of these conditional probabilities are known from the definition of the process. We have The explicit expressions for the probabilities in terms of f i (x t ), r i (x t ) and K i (τ ) are given by Eqs. (4)- (8).
The probability distribution P (x t+∆ |{x, k, , m} t ≤t ) is the product of δ-distributions given in Eq. (3), which for clarity we repeat here, We have To make the expressions that follow more compact we can define the product so that the path probability factorises into We will refer to Φ i,t,τ as the statistical weight associated with the combination of random variables It is not itself a probability, as the probabilities which compose it are conditioned on different random variables. All of the conditional probabilities which make up Φ i,t,τ in Eq. (A6) are only conditioned on variables which are either i) {x t } t≤t <t+τ or ii) other members of {k τ i,t , { τ,σ i,t } 0<σ<τ , m τ i,t }. The weight Φ i,t,τ is therefore only a function of {k τ i,t , { τ,σ i,t } 0<σ<τ , m τ i,t } and {x t } t≤t <t+τ . To simplify the notation we suppress these arguments. Later on this definition of Φ i,t,τ will allow us to factorise the generating functional and perform the sums for each factor separately. Inspecting Eqs. (4)-(8) we can identify the only combinations of values of k τ i,t , τ,0<s<τ i,t , and m τ i,t for which Φ i,t,τ is non-zero. They are as follows: (i) A reaction of type i with delay period τ does not fire at time t. In this case we have (ii) A reaction of type i with delay period τ fires at time t and is interrupted at time t + s. The values of the random variables for this case are and (iii) A reaction of type i with delay period τ fires at time t and is not interrupted. The random variables take the values where we have introduced the short-hand Dx ≡ t dx t . Using the exponential representation of the δ-function we can write with Dx ≡ The quantity Ψ i,t,τ is a function of {x t } t≤<t <t+τ , for clarity we have suppressed these arguments. In defining Ψ i,t,τ we have collected all occurrences of the variables indicated in the summation, {k τ i,t , { τ,σ i,t } 0<σ<τ , m τ i,t }. As explained beneath Eq. (A7), Φ i,t,τ is only a function of {k τ i,t , { τ,σ i,t } 0<σ<τ , m τ i,t } and of {x t } t≤<t <t+τ . This construction allows us to evaluate each Ψ i,t,τ separately to all other terms in the generating functional, which can be written as which for small ∆ can also be written as an exponential, cf. 1 + ∆f = e ∆f + O(∆ 2 ). Repeating this procedure for all Ψ i,t,τ leaves a generating functional for the discrete-time dynamics. We can then take the continuous-time limit, and arrive at Although the notation is involves a number of indices, it is straightforward to see how the generating functional can be extended to the case in which any delay reaction can be interrupted in multiple different ways. If reaction i is a delay reaction with which can be interrupted in Λ i possible ways, indexed by µ = 1, . . . , Λ i , then by extension of Eq. (2), We also have to modify the conditional probabilities of the random variables to take into account the different interruption reactions. If k τ i,t = 0 then τ,s i,µ,t = 0 for all s and µ, i.e. P ( τ,s i,µ,t = 0|k τ i,t = 0) = 1 ∀ 0 < s < τ , µ = 1, . . . , Λ i .