## Abstract

We study general stochastic birth and death processes including delay. We develop several approaches for the analytical treatment of these non-Markovian systems, valid, not only for constant delays, but also for stochastic delays with arbitrary probability distributions. The interplay between stochasticity and delay and, in particular, the effects of delay in the fluctuations and time correlations are discussed.

## 1. Introduction

Stochastic modelling plays an important role in many areas of science, such as physics, ecology or chemistry [1]. Stochasticity may appear due to the lack of complete knowledge about all the relevant variables, the precise dynamics of the system or the interactions with the environment. In some cases, one can obtain a compact description of a complicated system considering only a few relevant variables but at the expense of losing deterministic predictability. Often, probabilities for some fundamental processes can be assigned on the basis of symmetries and other considerations, or on empirical analysis, and the dynamics of the processes can be derived bottom-up.

Stochasticity appears together with delay terms in many situations of interest, such as gene regulation [2–4], physiological processes [5] or postural control [6,7]. The combined effects of stochasticity and delay are, however, not completely understood. From the mathematical point of view, stochastic processes including delay are difficult to analyse due to the non-Markovian character. Most of the previous approaches have focused on stochastic differential equations that consider continuous variables [8–12] or random walks in discrete time [13,14], where delay can be taken into account by increasing the number of variables. Models with discrete variables but continuous time are the natural description of many systems, such as chemical reactions, population dynamics or epidemic spreading. In some cases, discreteness can be a major source of fluctuations, not well captured by continuous models [15]. The approach with discrete variables and continuous time was used in [4,16–18]. Most often, the delay time is taken to be a constant with zero fluctuations. This is not very realistic in applications, since it is unusual to have a deterministic delay when the rest of the dynamics is stochastic. We will take this consideration into account by allowing the delay times to be random variables with arbitrary probability density functions.

In this work, we study some simple, yet general, stochastic birth and death processes including delay. We will develop three different approaches to the analytical study of these kinds of non-Markovian processes, in the general case of a stochastically distributed delay: a direct approach in §2*a*, an effective Markovian reduction in §2*b*,*c*, and a master equation approach, together with a time-reversal invariance assumption, in §3. The first direct approach method is interesting for its simplicity, but its application is limited to systems with first-order reactions and without feedback. The second one, effective Markovian reduction, is rather flexible and general and its development is one of the main advances of this paper. The last master equation approach, developed in [18], complements the previous, giving information about the full probability distribution; here we apply it to a particular feedback form and derive the phase diagram of the system. The main limitation of all the approaches is the need to assume that completion times for delayed reactions are independent random variables (independent of each other and of other variables of the system), although the initiation rates may depend on the state of the system, allowing, for example, for feedback and crowding effects, so we do not consider this limitation to be very relevant for practical applications. Although our methodology is rather general, we present it here using specific examples that have been grouped in two categories: delay in the degradation (§2) and delay in the creation (§3). We end the paper with a brief discussion and comments in §4. Some more technical details are left for appendices A and B.

## 2. Delayed degradation

We will start by studying simple stochastic birth and death processes that include delay in the degradation step. A process of this type was proposed in [4] as a model for protein-level dynamics with a complex degradation pathway.

### (a) Simple case

We consider first the simplest possible process including delayed degradation
2.1that is, a particle *X* is created at a rate *C* and disappears (‘dies’ or ‘degrades’) a time *τ* after being created. We allow the delay time *τ* to be randomly distributed, i.e. the lifetimes *τ* of the created particles are random variables, that for simplicity, we consider independent and identically distributed with probability density *f*(*τ*). Although not considered in this paper, the case of non-identically distributed delay times, in particular a probability density that depends on the time from birth, can also be treated. However, as commented earlier, the case of non-independent delay times does not seem to be tractable with the methods that we present later.

We first note that a distributed delay is completely equivalent to degradation at a rate that depends on the ‘age’ *a* (time from creation) of the particle, i.e. processes
2.2are equivalent if the rate *γ*(*a*) and the probability density of the delay *f*(*τ*) are related by
2.3with being the cumulative distribution of the delay time. This is so because *γ*(*a*) d*a* is the probability of dying at the time interval (*a*,*a*+d*a*), if the particle is still present at *a*, and so it is nothing but the probability *f*(*a*) d*a* that the delay time *τ* belongs to that same interval conditioned to the particle still being alive at time *a*, an event with probability . In the notation of [19], *γ*(*a*) is nothing but the conditional failure rate. We take *t*=0 as the time origin, so the number of alive particles at time *t* is *n*(*t*)=0 for *t*≤0. Let *P*(*n*,*t*) be the probability of *n* particles being alive at time *t*. In the remainder of this section, we assume that there is no feedback, in the sense that the creation rate *C* is independent of the number of particles *n*, but, for the sake of generality, we do allow it to be a function of time *C*(*t*). The non-feedback assumption allows us to obtain a full analytical solution. As shown in appendix A, independently of the form of the delay distribution, *P*(*n*,*t*) follows a Poisson distribution:
2.4with average . If the creation rate, *C*(*t*), is independent of time, a steady state is reached, in which the average number of particles is 〈*n*〉_{st}=*C*〈*τ*〉, again independently of the form of the delay distribution.

We will now compute the time correlation function. We shall see that its analytical expression does depend on the form of the delay distribution. We start from the relation
2.5with *n*_{new} (*n*_{old}) particles created after (before) *t*. *n*_{new} can be computed exactly as before (now taking *t* as the time origin), so we have
2.6The evolution of the number of particles already present at *t* depends on the age *a* of these particles. Their survival probability until time *t*+*T* can be written as follows:
2.7where we used *P*(*a*|*b*)=*P*(*a*;*b*)/*P*(*b*), so we find
2.8From this, one easily obtains the correlation function: *K*[*n*](*t*,*T*)=〈*n*(*t*)*n*(*t*+*T*)〉−〈*n*(*t*)〉〈*n*(*t*+*T*)〉=〈*n*(*t*)〈*n*(*t*+*T*)|*n*(*t*)〉〉−〈*n*(*t*)〉〈*n*(*t*+*T*)〉
2.9If *C*(*t*)=*C*, independent of time, a steady state can be reached with correlation function . For a constant rate *γ*, which would be equivalent to an exponential delay distribution *f*(*τ*)=*γ* e^{−γτ}, it has the usual exponential decay *K*_{st}[*n*](*T*)=(*C*/*γ*) e^{−γT}. For a fixed delay time *τ*_{0}, corresponding to *f*(*τ*)=*δ*(*τ*−*τ*_{0}), the correlation function is a straight line *K*_{st}[*n*](*T*)=*C*(*τ*_{0}−*T*) for *T*<*τ*_{0} and *K*_{st}[*n*](*T*)=0 for *T*≥*τ*_{0}. For other distributions of delay time, the correlation function adopts different forms, but it is always monotonically decreasing. In figure 1, we plot the correlation function for two different types of distribution of delay, for different values of the variance of the delay. We see that the distribution with fatter tail displays a slower asymptotic decay, and that the decay is slower as the variance of the delay increases. Numerical simulations, performed with a conveniently modified version of the Gillespie algorithm [20,21], are in perfect agreement with this exact result, providing a check of its correctness. We remark that the functional form of the decay of the correlation function depends on the delay distribution and can differ from the exponential decay found in systems without delay.

### (b) More elaborated case

We now consider a process including both instantaneous and delayed degradation steps
2.10that is, particles are created at a rate *C* and each particle can be eliminated by two processes: (i) instantaneous degradation at a rate *γ* and (ii) delayed degradation, initiated at a rate *D* but completed only a time *τ* after initiation. Again, we will allow the delay-degradation times to be random variables that, for simplicity, will be independent and identically distributed with probability density function *f*(*τ*).

For the process to be completely defined, one has to specify if a particle that initiates delayed degradation at time *t* and thus will disappear at *t*+*τ* (this kind of particle will be called ‘infected’), can also disappear before the completion of this reaction, through instantaneous degradation. In the most general case, this can happen at a rate *γ*′, not necessarily equal to *γ*. Note that, in the case of first-order degradation (*γ*′ not dependent on the number of particles *n*), this instantaneous degradation is completely equivalent to a system with *γ*′=0, after modifying the distribution of the delayed-degradation times in the following way:
2.11That is, when instantaneous degradation is added to infected particles, the probability that the lifetime is equal to *τ* has two contributions: (i) a particle initially has a lifetime *τ* (probability density *f*(*τ*)) and survives up to this time (an event with probability e^{−γ′τ}) and (ii) a particle has a lifetime larger than *τ* (probability ), but survives up to *τ* (probability e^{−γ′τ}) and then undergoes instantaneous degradation (at rate *γ*′). The consideration of these two contributions leads straightforwardly to equation (2.11). We see that omitting first-order instantaneous degradation of infected particles involves no loss of generality, given that the treatment is valid for general distributions of delay.

If *D* and *γ* are independent of *n*, the process is equivalent to the one-variable system discussed in §2*a* with a conveniently modified distribution of delay
2.12This comes from the fact that a particle may disappear at time *τ*, because it did not disappear or was infected before and is degraded instantaneously (probability density e^{−(γ+D)τ}*γ*) or because it got infected at some previous time (*t*′) with an appropriate lifetime (*τ*−*t*′, probability density ). This includes as particular cases the ones studied in [17,22]. The results of §2*a* allow us to obtain the full solution also in the general case of distributed delay. If *D* or *γ* depend on *n* the processes are not anymore equivalent, two variables are necessary and a new approach is needed for the analysis. In the following, we develop this method. We will also consider the case in which the creation rate *C* depends on the number of particles.

The full process corresponds to the following two-variable system:
2.13where we have split the proteins into two types: *X*_{I} are infected particles that will die precisely at a time *τ* (itself a stochastic variable) after being infected and *X*_{A} are non-infected (‘active’) particles (so *X*=*X*_{A}∪*X*_{I}). We allow the rates to depend on *n*_{A}, the number of *X*_{A}, active, particles, but not on *n*_{I}, the number of *X*_{I}, infected, particles which are considered to be ‘inert’; this condition will be relaxed in §2*c*. Following the study of Miekisz *et al*. [17], we have introduced the auxiliary particles *Z* whose number is given by the stochastic variable *n*_{Z}(*t*). Note that each time a particle of type *X*_{I} appears, a particle of type *Z* also appears, so the process *n*_{Z}(*t*) is similar to the process *n*_{I}(*t*) but without ‘deaths’ (that constitute the non-Markovian part of the process). The introduction of *Z* will allow us to obtain the properties of *n*_{I} by using the relation
2.14where the discrete process *n*_{Z}(*t*) is a sequence of step (Heaviside) functions and its derivative must be understood as a series of Dirac-delta functions. Here, we have introduced the family of ‘survival’ stochastic processes *s*(*t*′,*t*) defined in the following way: first, for each *t*′ we obtain a value of *τ*(*t*′) independently drawn from the distribution *f*(*τ*). Next, we set *s*(*t*′,*t*)=1, if *t*∈(*t*′,*t*′+*τ*(*t*′)) and *s*(*t*′,*t*)=0, otherwise. This can be considered as the indicator function of a virtual^{1} particle that is infected at *t*′ and survives up to a time *t*′+*τ*(*t*′). It follows from the definition that
2.15and
2.16

Expressions (2.14)–(2.16) are the main advances of this section and provide us with the necessary tools to derive the main properties of the stochastic process (2.10). In the case considered in [17], there is a fixed delay (*f*(*τ*)=*δ*(*τ*−*τ*_{0})) and no instantaneous degradation of infected particles (*γ*′=0), so one has simply *n*_{I}(*t*)=*n*_{Z}(*t*)−*n*_{Z}(*t*−*τ*). The inclusion of the survival process *s*(*t*′,*t*) allows us to consider the general case of distributed delay and rates depending on the state of the system.

Note that the process followed by {*n*_{A},*n*_{Z}} is a simple Markovian birth and death process:
2.17This is so because the variable *n*_{I} does not affect the creation or degradation of other variables. The case in which variables that undergo delay degradation affect other variables will be considered in §3. In the present case, the properties of *n*_{Z} can be obtained using Markovian methods, and the properties of the variable *n*_{I} can be derived afterwards using (2.14)–(2.16). In particular, the first moments follow:
2.18and
2.19Using standard Markovian methods [1], one can prove that the process {*n*_{A},*n*_{Z}} is described by the master equation
2.20with *E*_{i} the step operator, . In this section, we allow the creation rate *C* to depend on the number of *X*_{A} particles, constituting a feedback term on the number of ‘active’ particles. From the master equation, one easily derives the equations for the moments, the first of them read
2.21
2.22
2.23
2.24
and
2.25
In the case that *C*(*n*_{A}) is a linear function of *n*_{A} and *γ* and *D* do not depend on *n*_{A} (and none of them depend on *n*_{I} or *n*_{Z}), the system of equations is closed and can be solved. For nonlinear systems, we will make use of van Kampen's expansion [1]. This is a standard systematic expansion of the master equation, that consists of assuming a deterministic, and a stochastic part for the variables that scale differently with a large parameter, *Ω* (typically the volume or system size), i.e. *n*_{A}=*Ωϕ*_{A}(*t*)+*Ω*^{1/2}*ξ*_{A}, *n*_{Z}=*Ωϕ*_{z}(*t*)+*Ω*^{1/2}*ξ*_{Z}. One can then write the master equation for the new variables *ξ*_{A}, *ξ*_{Z} and expand in powers of *Ω*^{−1/2}. The method is generically valid, provided that the rates depend on the variables only through *n*_{α}/*Ω* (plus higher orders in *Ω*^{−1}; a common factor depending on *Ω* multiplying all rates is also acceptable), which is fulfilled by most systems of interest, and that the macroscopic equations have a steady state as a single attractor. The equations for the macroscopic components are
2.26
and
2.27
The stochastic contributions, to first order in *Ω*^{−1/2}, read
2.28
2.29
2.30
2.31
and
2.32
with , . Usually, for the ansatz about the scaling of the variables to work (and so the expansion), the equations for the macroscopic components must have a single stable fixed point. In this case, however, the equation for *ϕ*_{z} does not have a fixed point, and *ϕ*_{z}(*t*) and grow without bound. This growth, nevertheless, is consistent with 〈*n*_{Z}(*t*)〉/*σ*^{2}[*n*_{Z}](*t*)=*O*(*Ω*^{0}), and the expansion can still be applied.

Equations (2.28)–(2.32) are a system of closed linear equations and so can always be solved. To compute the time correlations of *n*_{I} from equation (2.19), we need the time correlations of *n*_{Z}. We note that
2.33
2.34
and that 〈*n*_{Z}(*t*_{2})|*n*_{Z}(*t*_{1}),*n*_{A}(*t*_{1})〉 (for *t*_{2}>*t*_{1}) can be obtained integrating (2.21)–(2.22) or (2.28)–(2.29). In the general, nonlinear, case, using first-order van Kampen's expansion, one obtains, over the steady state
2.35
with and *ϕ*_{A,st} the solution of . The derivative that appears in (2.19) is
2.36
with .

Putting all the pieces together, one finally obtains
2.37
Proceeding in a similar way, one can derive
2.38
2.39
and
2.40
with *K*_{st}[*n*_{u},*n*_{v}](*t*)≡〈*n*_{u}(*t*_{0}+*t*)*n*_{v}(*t*_{0})〉_{st}−〈*n*_{u}〉_{st}〈*n*_{v}〉_{st}. This finally allows one to express the correlation function for the total number of particles, *n*=*n*_{A}+*n*_{I}, as
2.41
In this case, the average of *n* again depends only on the average delay, 〈*n*〉_{st}=*Ωϕ*_{A}(1+*D*〈*τ*〉), but the second moment depends on the delay distribution in a more complicated way, through factors involving the integral of .

In figure 2, this result is compared with numerical simulations, showing a very good agreement. Note that the treatment of the delayed reactions is exact, the only approximation coming from the use of van Kampen's expansion, which is needed when nonlinearities are present, but whose error scales as *Ω*^{−1/2}. Similar to the previous case, the process in which the distribution of delay has a fatter tail shows slower decay for the correlation function.

### (c) Full feedback

We now consider the case in which the creation rate depends on all present particles
2.42
with *n* the total (inert+active) number of *X* particles. As noted before, this single-variable model can account for instantaneous plus delayed degradation, in the case that the degradation and ‘contagion’ rates, *γ* and *D* in the previous notation, do not depend on the state of the system. For simplicity, we restrict our attention to this case. This process can be treated with the approach of §2*b* by introducing the additional variable *Z*,
2.43
with *n*_{Z}(*t*) the corresponding random variable giving the number of *Z* particles. We see that
2.44
with *s*(*t*′,*t*) the same as in §2*b*. The probability distribution for {*n*,*n*_{Z}} follows a master equation of the form
2.45
Details of the derivation of the master equation in systems with delay are given in appendix B. Here, , with *f*(*t*) the probability density of the delay distribution, although, since we are only interested in the properties of the variable *n*, we will not be using this expression. The key step in this case is to note that equation (2.45) allows us to derive the statistical properties (moments and correlations) of *n*_{Z}(*t*) as a function of those of *n*(*t*). Then, using (2.44) we will be able to self-consistently derive the properties of *n*. More specifically, the approach proceeds as follows.

Summing equation (2.45) over *n*, we can obtain an equation for the evolution of *P*(*n*_{Z},*t*), but that still depends on *n* (in this step the contribution of the second term in equation (2.45) vanishes):
2.46
The two times probability distribution *P*(*n*_{Z1},*t*_{1};*n*_{Z2},*t*_{2}) follows a similar equation. Conditioning carefully, summing over the variable *n* and considering separately the case *t*_{1}=*t*_{2} (which turns out to be singular), we find
2.47
From (2.44), we easily obtain
2.48
and
2.49
while (2.46) and (2.47) imply
2.50
and
2.51
And we finally obtain the following set of integral equations for the moments:
2.52
and
2.53
In the case of linear feedback, *C*(*n*)=*a*+*bn*, this system of equations is closed. For nonlinear systems, one can use van Kampen's expansion as explained above. In the steady state, one finds
2.54
and
2.55
Equation (2.54) shows that the steady-state number of particles depends only on the average delay. Equation (2.55) shows that the correlations depend on the delay distribution in a non-trivial way. The analysis of this equation is left for future work.

## 3. Delayed creation

We now turn our attention to the case in which the creation reaction, that is initiated stochastically, takes a finite time to be completed. For simplicity, we assume that the degradation reaction is instantaneous. Schematically, we have
3.1
In this case, if the creation rate does not depend on the number of particles, *n*, then the delay in the creation is completely irrelevant, since the probability that a new particle appears at time *t* is equal to the probability that its creation started at a time *t*−*τ*, but this is equal to the probability that a particle starts its creation at time *t* (with a shift in the time if *C* is time-dependent), so the process is completely equal to one with instantaneous creation.

Following [18], we will adopt here an approach different from that of the previous sections, that, besides the moments, will allow us to obtain an expression for the full probability distribution. For completeness, we will explain here the method in some detail. For additional considerations, the reader is referred to [18]. In appendix B, it is shown that the master equation of the process (3.1) is
3.2
The master equation (3.2) can be written as
3.3
where the effective creation rate, , is given by
3.4
The conditional probability *P*(*n*,*t*|*n*_{0},*t*_{0}) follows a master equation identical to (3.2) with all the probabilities conditioned to *n*_{0} at time *t*_{0}. From it, and using that , we obtain the following evolution equation for the conditional average:
3.5
for *t*≥*t*_{0}, with initial condition 〈*n*(*t*_{0})|*n*(*t*_{0})〉=*n*(*t*_{0}).

The knowledge of the steady value allows the calculation of the steady-state probabilities *P*_{st}(*n*), obtained by imposing ∂*P*(*n*,*t*)/∂*t*=0 in equation (3.3), as [1]
3.6
where *P*_{st}(0) is fixed by the normalization condition. All that is left to do now is to compute the effective creation rate .

The effective creation rate will be computed using expression (3.5). In the general case of nonlinear creation rate, we will use van Kampen's expansion to linearize *C*(*n*) around the macroscopic component of *n*. We have *C*(*n*)=*ΩC*(*ϕ*)+*Ω*^{1/2}*C*′(*ϕ*)*ξ*, so
3.7
and using (3.5) we obtain
3.8
and
3.9
*t*≥*t*_{0}. Equation (3.8) is in general a nonlinear integro-differential equation that can be difficult to solve. Here, however, we will focus on the cases in which (3.8) has a stable steady state as a single attractor, which is the solution of *γϕ*=*C*(*ϕ*). This is the regimen in which the validity of van Kampen's expansion is guaranteed.

We reach now a delicate point. Equation (3.9) is a (linear) integro-differential equation. To solve it, we would need an initial condition in the whole interval but we only know a one-time condition 〈*ξ*(*t*=*t*_{0})|*ξ*(*t*_{0})〉=*ξ*(*t*_{0}). We will circumvent this difficulty by assuming that, over the steady state, the system is statistically invariant under time inversion, which implies 〈*ξ*(*t*_{0}+*t*_{1})|*ξ*(*t*_{0})〉=〈*ξ*(*t*_{0}−*t*_{1})|*ξ*(*t*_{0})〉. This condition, together with the value of *ξ* at time *t*_{0}, allows us to find the solution of (3.9). This solution, together with (3.7) and the definition (3.4), yield the effective creation rate . The time-reversal invariance assumption in the steady state is fulfilled by any Markovian system that follows detailed balance. Our system follows detailed balance (as any one-step process [1]), but, owing to the presence of delay, it is not Markovian. So the time-reversal invariance is an assumption, whose validity needs to be checked. It was shown in the study of Lafuerza & Toral [18] that in this system the assumption is approximately valid.

In the case of constant delay, *f*(*t*)=*δ*(*t*−*τ*), the time-reversal symmetric solution of (3.9) is [4,18]
3.10
and
3.11
and using (3.7) we finally obtain
3.12
where *ϕ*_{st} is the steady-state solution of (3.8). From equation (3.6), one can obtain the steady-state probabilities *P*_{st}(*n*). The mean value and variance are given by
3.13
and
3.14
From (3.14), one can see that, interestingly, in the case of negative feedback (*C*′(*ϕ*_{st})<0), as the delay is increased the fluctuations change from sub-Poissonian (*σ*^{2}<〈*n*〉) to super-Poissonian (*σ*^{2}>〈*n*〉). This is illustrated in figure 3, where we also show the line of the Hopf bifurcation that the deterministic system undergoes due to the delayed negative feedback [23]. It is usually obtained that a negative feedback reduces the magnitude of the fluctuations [24], but when delay is present we see that this negative feedback can change totally its effect, giving rise to an increase of the fluctuations.

The time correlation function can also be obtained from (3.10), as 3.15 In the case of negative feedback, it becomes non-monotonic, developing peaks of alternating sign at approximately multiples of the delay, signalling the presence of stochastic oscillations. For positive feedback, the time correlation is always positive, but not necessarily monotonic.

We will finish by noting that the ‘effective Markovian reduction’ method used in §2 can also be used for the case of delay in the creation with feedback. To be completely general, we allow two delays, one in the creation (with probability density *f*_{c}(*t*)) and one in the degradation (with probability density *f*_{d}(*t*)). The process is schematized as follows:
3.16
with *τ*_{c/d} random variables distributed according to *f*_{c/d}(*t*). With the addition of two new variables, the process can be rewritten as
3.17
which allows us to note that
3.18
In this case, the survival function is defined as: , if *t*∈(*t*′+*τ*_{c}(*t*′),*t*′+*τ*_{c}(*t*′)+*τ*_{d}(*t*′)) and , otherwise, *τ*_{c}(*t*′) and *τ*_{d}(*t*′) being random times obtained from the corresponding probability distribution functions *f*_{c}(*τ*_{c}) and *f*_{d}(*τ*_{d}). is equal to one if a virtual particle that initiated its creation at time *t*′ finished it at some intermediate time *t*′′<*t* and since then had a lifetime greater that *t*−*t*′′, so that it is still alive at *t*, being zero otherwise. It follows that
3.19
and
3.20
In the case that the creation rate *C*(*n*) does not depend on the number of *X*-particles, the number of *Z*-particles follows a Markovian process (Poisson process), and the properties of *n* can be derived from (3.18). If the creation rate depends on the number of *X*-particles, i.e. if feedback is present, the properties of *n*_{Z} can be derived formally as a function of *n* and then the properties of *n* can be derived self-consistently through (3.18), as done in §2*c*.

## 4. Comments and conclusions

In this paper, we have analysed general stochastic birth and death models that include delay. We have presented three different methods that together constitute a general toolbox to study stochastic models including delay.

In §2*a*, we have shown that when the creation rate is independent of the state of the system (no feedback) and the initiation of the delayed degradation and the instantaneous degradation are first-order reactions (rate not depending on the state of the system), the process can be solved fully in an exact fashion for general distributions of delay, showing always Poissonian character and a monotonically decreasing time correlation function given by (2.9).

In §2*b*,*c*, we have considered a more general process with delay in the degradation step, allowing the initiation of the delay degradation and the instantaneous degradation to be higher order reactions, as well as the presence of feedback in the creation rate. The method allows one to reduce the system to a Markovian one, where usual techniques can be used. Explicit expressions for the time correlation for general delay distributions were obtained. In this case, the correlation might be non-monotonic, if feedback is present, but typically decreases monotonically.

Section 3 shows that when the delay appears in the creation reaction and feedback is present, the delay typically has more dramatic consequences. In the case of fixed delay, it is shown that for negative feedback, the fluctuations are amplified as the delay increases, going beyond the level found when no feedback is present, and the time correlation function becomes oscillatory, alternating between positive and negative values at approximately multiples of the delay. In the positive feedback case, again for fixed delay, the fluctuations are reduced with increased delay and the time correlation function remains always positive.

## Funding statement

We acknowledge financial support by MINECO (Spain), Comunitat Autonoma de les Illes Balears, FEDER, and the European Commission under projects FIS2007-60327 and FIS2012-30634. L.F.L. is supported by the JAE Predoc Program of CSIC.

## Appendix A. Calculation of *P*(*n*,*t*) in the simple case of delayed degradation

We start by considering the case *n*=0. For the sake of simplicity, we focus on the case with creation rate, *C*, independent of time, but the generalization to time-dependent *C* is straightforward. Since the time origin is taken at *t*=0, the probability of observing zero particles at time *t*>0 is equal to the following limit:
A1
with *Δt*≡*t*/*M* playing the role of a small time increment and *t*_{i}≡*iΔt*. This expression follows from the fact that, in order to find the system with zero particles at time *t*, in every previous infinitesimal time interval (*t*′∈[*t*_{i},*t*_{i+1}),*i*=0,…,*M*−1) one of the following two (incompatible) events must take place: either a particle is not created (probability 1−*CΔt*) or a particle is created with a lifetime smaller that *t*−*t*_{i} (probability *CΔtF*(*t*−*t*_{i})). We now have
A2
with , so we find
A3

Following a similar line of reasoning, *P*(*n*,*t*) can be computed as
A4
This expression results from the consideration of choosing the times (*t*_{i1},…,*t*_{in}) at which the *n* particles are created and survive up to *t*. The *l*th particle is created with probability *CΔt* and survives up to *t* with probability . The other factor comes from the fact that at the other time intervals either a particle is not created or it is created but dies before *t*.

Using
A5
and replacing the sums by integrals in the limit
A6
we finally obtain
A7
that is, a Poisson distribution with average . In the steady state (found as the limit ), the average becomes 〈*n*(*t*)〉=*C*〈*τ*〉. Remarkably, this Poissonian character is completely independent of the form of the delay distribution. As commented above, this result can be easily generalized to the case in which the creation rate depends on time, , obtaining again a Poisson distribution with average .

## Appendix B. Derivation of the master equation in a system with delay

Here, we derive the master equation of the process (3.1). We consider first the case of fixed delay *τ*. We start with the following identity:
B1
It is immediately seen that *P*(*n*,*t*+*Δ*;*n*+1,*t*)=*γ*(*n*+1)*ΔP*(*n*+1,*t*). In the case of fixed delay, the second sum can be evaluated introducing a tree-times probability as
B2
Now, *P*(*n*,*t*+*Δ*|*n*−1,*t*;*n*′,*t*−*τ*)=*C*(*n*′)*Δ*+*o*(*Δ*). Expanding in a similar way the term *P*(*n*,*t*+*Δ*;*n*,*t*), and taking the limit , we can obtain the master equation of the process
B3
In the case of distributed delay, we start considering a discrete distribution of delays, i.e. *τ*=*τ*_{1},…,*τ*_{M} with corresponding probabilities *f*(*τ*_{1}),…,*f*(*τ*_{M}). The continuum limit can then be obtained making . The creation term in (B1) can be written as
B4
Now, , that is, the probability that a particle started its creation at time *t*−*τ*_{i} with a creation time equal to *τ*_{i}. Replacing in the previous equation and performing the appropriate sums, we obtain
B5
that in the continuum limit reduces to . Considering in a similar way the other terms in (B1) and taking the limit one can obtain the master equation for distributed delay (3.2).

## Footnotes

One contribution of 15 to a Theme Issue ‘Dynamics, control and information in delay-coupled systems’.

↵1

*s*(*t*′,*t*) is defined for all*t*′, regardless of whether a particle is actually infected a time*t*′. However, it contributes only to (2.14) if a particle is actually infected at time*t*′, since only then d*n*_{Z}(*t*′)/d*t*′≠0.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.