# Endogenous Quasicycles and Stochastic Coherence in a Closed Endemic Model

###### Abstract

We study the role of demographic fluctuations in typical endemics as exemplified by the stochastic SIRS model. The birth-death master equation of the model is simulated using exact numerics and analysed within the linear noise approximation. The endemic fixed point is unstable to internal demographic noise, and leads to sustained oscillations. This is ensured when the eigenvalues () of the linearised drift matrix are complex, which in turn, is possible only if detailed balance is violated. In the oscillatory state, the phases decorrelate asymptotically, distinguishing such oscillations from those produced by external periodic forcing. These so-called quasicycles are of sufficient strength to be detected reliably only when the ratio is of order unity. The coherence or regularity of these oscillations show a maximum as a function of population size, an effect known variously as stochastic coherence or coherence resonance. We find that stochastic coherence can be simply understood as resulting from a non-monotonic variation of with population size. Thus, within the linear noise approximation, stochastic coherence can be predicted from a purely deterministic analysis. The non-normality of the linearised drift matrix, associated with the violation of detailed balance, leads to enhanced fluctuations in the population amplitudes.

###### pacs:

87.18.Tt, 05.10.Gg, 02.50.Ga## I Introduction

Two and a half centuries ago, D. Bernoulli Bernoulli (1760) used a nonlinear ordinary differential equation to study the effect of cow-pox inoculation on the spread of smallpox. This was one of the earliest examples of the mathematical study of epidemics. This field of study continues to hold the interest of the scientific community especially in the light of recent outbreaks of viral pandemics like SARS and H1N1. Kermack and McKendrick in their seminal paper Kermack and McKendrick (1927) put forward the classic Susceptible-Infected-Recovered (SIR) model of the spread of epidemics which, like most early epidemic models, assumes a homogeneously mixed population. More recent work focuses on the geotemporal spread of epidemics, especially on model networks Kuperman and Abramson (2001); Pastor-Satorras and Vespignani (2001). However, homogeneous mixing models still prove to be useful Keeling (2005) and have been used to study various outbreaks of diverse size, fatality and chronology. Examples range from the study of the plague in the village of Eyam in 1665-66 Raggett (1982) to the Bombay plague of 1905-06 Kermack and McKendrick (1927) and the influenza epidemic in an English boarding school Murray (2002).

Mathematical models like the SIR model are usually analysed deterministically and are only exactly valid when the size of the population under consideration is exceedingly large. Fluctuations due to finite population sizes or due to external causes can give rise to phenomena which cannot be captured by deterministic mean-field models and necessitates the use of stochastic models. Bartlett Bartlett (1956, 1957) was one of the first to realise that a stochastic description was necessary to explain the periodic recurrence of measles, a phenomenon which could not be explained by deterministic models Soper (1929); Wilson and Worcester (1945). Bartlett formulated Bartlett (1956) a stochastic version of the SIR model to describe the periodic recurrence of measles.

The mechanism for the generation of sustained oscillations in population dynamics has been analysed within the stochastic framework Nisbet and Gurney (1976) which concentrates on external fluctuations as the noise source. However, finite-sized populations give rise to fluctuations whose relative amplitude is of the order of the inverse of the square root of the size of the population Schrödinger (1992). The role played by this internal noise, arising out of demographic stochasticity, in the generation of sustained oscillations has been studied in a prey-predator model using a master equation approach by McKane and Newman McKane and Newman (2005). They have used the expansion method due to van Kampen van Kampen (2007) in their analysis, which provides a systematic way of deriving the phenomonological equations due to Bartlett Bartlett (1956). Alonso et al. Alonso et al. (2007) used similar techniques in an open model of infectious diseases within the homogeneous mixing assumption, while Rozhnova and Nunes Rozhnova and Nunes (2009) applied systematic expansion to a closed epidemic model on networks, using a pair approximation. The oscillations generated and sustained by internal noise are called endogenous resonant quasicycles and are qualitatively different from stochastic oscillations forced by external periodicities which are exogenous Nisbet and Gurney (1982). The quality or coherence of these oscillations are intuitively expected to vary monotonically with the size of the population or equivalently, the relative noise amplitude. However, it has been observed in various theoretical models including the Fitz Hugh-Nagumo Pikovsky and Kürths (1997) and gene circuit models Hilborn and Erwin (2008) that the regularity or coherence of oscillations is small for low and high noise amplitudes and reaches a maximum for an intermediate value. This phenomenon is called stochastic coherence or coherence resonance and has also been observed in optical laser experiments Giacomelli et al. (2000).

In this work we analyse the generation of quasicycles due to internal noise, as well as the non-trivial variation of the quality of oscillation with respect to population size, in a closed epidemic model under the homogeneous mixing assumption. The closed system is relevant in many epidemiological situations, for instance in boarding houses Murray (2002), or island communities, where no inflows or outfluxes occur. Further, the conservation of populations, as implied by a closed system, allows one to deal with a lower-dimensional problem. We exploit this in a systematic manner and show how the master equation can be marginalised using the conservation constraint. The existence of an endemic fixed point allows a two-stage linearisation procedure to be carried out on the model. The linear noise approximation, followed by a further linearisation about the endemic fixed point, reduces the model to the standard multivariate Ornstein-Uhlenbeck (OU) form. Exploiting the linear and Gaussian character of the multivariate OU process then allows for stochastic behaviour to be predicted from the deterministic part of the dynamics, in a spirit similar to the Onsager regression method of equilibrium statistical mechanics.

Below we review (section II) the deterministic analysis of the SIRS model, emphasising the behaviour around the stable endemic fixed point. Perturbations about this fixed point decay either monotonically or in a damped oscillatory fashion. Stochastic analysis (sections III and IV), however, shows that demographic noise destabilises this endemic fixed point, generating and sustaining oscillations. We show the existence of stochastic coherence (section V) by analytical means and then confirm it numerically. We show from purely deterministic analysis that there is stochastic coherence if the absolute value of the ratio of the imaginary and real parts of the eigenvalues of the linearised drift matrix shows a maximum when scanned against noise amplitude. The position of this maximum gives the population size corresponding to stochastic coherence for the relevant parameter values. We also show that (section VI) it is not possible to observe endogenous quasicycles unless detailed balance is violated. Finally we look at the non-normal aspect of the governing dynamics (section VII) and observe that the fluctuation amplitudes of the populations increase due to non-normality.

## Ii SIRS linear deterministic analysis

The classic SIR model for infectious diseases (S stands for susceptibles, I for infected and R for recovered) considers the population to be homogeneously mixed and constant in total number Kermack and McKendrick (1927). The SIRS model is a variant of the SIR model where the recovered section of the population lose their immunity after a delay and become susceptible. The nonlinear ODE system of the form , where , describing the SIRS model is constrained by the fixed population size and is hence a closed system.

(1) | |||||

The rate of infection is , the rate of recovery is while is the rate of loss of immunity. The fixed point () is given by

(2) |

The steady state with zero infected is not of interest in the present study. The fixed point is endemic with non-zero infected in the steady state () when the condition is satisfied.

Since there is a constraint in the system, , the system is effectively a system with .

(3) | |||||

The dynamics of small perturbations, , about the fixed point are described by the linear ODE system . Here is the Jacobian matrix at the fixed point and is given by

(4) |

Its eigenvalues are

(5) |

the real parts of which are always negative for an endemic steady state since . Hence the endemic fixed point is always asymptotically stable. Perturbations about the fixed point decay monotonically if the eigenvalues are purely real and in an oscillatory fashion if they are complex. These correspond, respectively, to overdamped and underdamped decay. In Figure (1), we plot both time traces and phase portraits of S and I showing the underdamped and overdamped cases. Figure (2) is a state diagram of the model, showing the ratio . The region of complex eigenvalues, corresponding to underdamped decay, is bounded by the contours labelled by .

## Iii SIRS linear stochastic analysis

Relative fluctuations about the deterministic expected values vary as the inverse of the square root of the number of interacting entities and thus become important when the entitites are few in number. Often one finds that this is indeed the case in biological systems Schrödinger (1992). Our present study concerns populations where fluctuations due to demographic stochasticity cannot be ignored and mean-field deterministic analysis fails to capture its non-trivial contributions. It then becomes necessary to employ stochastic methods to reliably understand the role of fluctuations.

We begin by writing down the birth-death master equation (ME) of the SIRS model. Let the state of the system at any time be given by the vector where is the number of individuals in each class ( for S, for I and for R). The general birth-death ME is then Gardiner (2004)

(6) |

Here is the conditional probability for the system to be in the state given some fixed initial state, and are the birth and death rate terms and is the vector denoting the change in the number of entities in the -th reaction. For the SIRS model we have

(7) |

The variable is constrained by the fixed population size . Incorporating this constraint within the ME allows us to work with a system. The partial time derivative of vanishes if for all . We marginalise with respect to one of the variables (here we choose ) taking the population size as parameter: . Modifying the birth terms and the state change vectors appropriately, i.e. replacing by and writing the as vectors, we get the marginalised ME.

(8) |

The transition probability for the infection step is non-linear and as such the ME is not solvable analytically. However, it is possible to simulate the ME using the Doob-Gillespie stochastic simulation algorithm (SSA) Doob (1945); Gillespie (1976, 1977). This generates an exact sampled trajectory of the jump stochastic process described by the ME. We non-dimensionalise time by working in units of . Figure (3) shows the numerical simulation of the susceptibles using the SSA, compared with a deterministic solution of the ODE system. The demographic fluctuations induce and sustain approximate cycles in the populations, a feature absent in the deterministic model.

In the absence of exact solutions, we try to characterise these fluctuations within an approximation method due to van Kampen van Kampen (2007) which replaces the jump process with a stationary multivariate Ornstein-Uhlenbeck process. The Gaussian nature of the process can then be utilised to obtain analytical solutions for the fluctuation properties, while the linear nature of the process can be utilised to make connections between the deterministic and fluctuating dynamics.

We expand the variables in the population size (the large parameter of the approximation method) so that the size of the jumps decreases as the population is increased,

(9) |

where is the mean value of and denotes the fluctuations around the mean. Assuming the fluctuations obey a diffusion process about the mean yields a Fokker-Planck equation (FPE) for the fluctuations,

(10) |

where repeated indices indicate summation, and . This is the linear noise approximation. The elements of the drift vector and the diffusion matrix are given, following the prescription in Gardiner Gardiner (2004), as

(11) |

( being the component indices). Linearising a second time about the endemic fixed point we get the FPE of a stationary multivariate Ornstein-Uhlenbeck process

(12) |

where and are the elements of the linearised drift and diffusion matrices. For the SIRS model, their values are (from Equation (11) after putting )

(13) |

(14) |

We note that this linearised drift matrix is identical to the linearised Jacobian matrix (Equation 4) obtained from the deterministic analysis and hence the two matrices share the same spectrum. This allows us to predict, under the two-stage linearisation procedure, the existence of non-trivial stochastic phenomena like noise-induced quasicycles and stochastic coherence purely from a deterministic analysis of the spectral structure of the linearised Jacobian. We shall discuss this important point in greater detail in sections IV and V. This also allows us to use the terms “linearised drift matrix” and “linearised Jacobian matrix” interchangeably.

The multivariate Ornstein-Uhlenbeck process has exact solutions for both stationary and transition probability densities. Both are multivariate Gaussians, fixed by the equal time covariance matrix and the matrix of time correlations , where the double angular brackets denote the cumulant van Kampen (2007). can be obtained by solving the steady state Einstein relation Gardiner (2004); van Kampen (2007).

(15) |

This has the form of a matrix Lyapunov equation, and can be solved using a method first proposed by Barnett and Storey Barnett and Storey (1966) in the context of linear control systems. We note that Equation (15) can be written as the sum of a matrix and its transpose where is the anti-symmetric matrix . We can solve for in terms of and using the relation

(16) |

which is obtained by eliminating from the Einstein relation and using the definition of . Since is anti-symmetric, it is specified by a single parameter when it is of size . This parameter can be obtained directly from Equation (16), since both and are two-dimensional matrices and are known. For higher dimensions, matrix decompositions are convenient when solving for .

For the SIRS model (using Equations 13 and 14), we have

(17) |

Knowing , and we can now write down the covariance matrix

(18) |

which for the SIRS model is

(19) |

Having obtained the matrix , the matrix of time correlations follows as

(20) |

The stochastic SIRS model, in the linear noise approximation, is completely specified by and . In the next section we use quantities derived from these to examine the model for signatures of oscillatory behaviour.

## Iv Noise-induced oscillations: Endogenous Quasicycles

The trace of the variation of the populations with time shown in Figure (3) is strongly suggestive of sustained oscillations. This can be verified quantitatively by measuring the power spectral density (PSD) of the population time series. A peak in the PSD indicates the presence of oscillations. The PSD matrix, in terms of the linearised drift and diffusion matrices for a multivariate Ornstein-Uhlenbeck process is

(21) |

where is the identity matrix. The diagonal elements of this matrix give an estimate of the periodicity in the relevant variables (here S and I). The for the SIRS PSD are

(22) |

where , , , and .

In Figure (4) we plot the PSD for both S and I, comparing numerical simulation with Equation (22). A peak is clearly visible for parameters corresponding to underdamped dynamics. The peak disappears for overdamped dynamics as shown in the inset. The peak frequency (around ) corresponds to the period () of the numerical time-trace (Figure 3). The excellent agreement between numerics and analytics provides a post-facto justification of the linear noise approximation for this problem.

The PSD has peaks at real frequencies if and only if the extremum condition has real roots. The regions of parameter space for which this occurs are bounded by contours labelled “” and “” in Figure (2). These are sufficient conditions for the existence of quasicycles. This approach has been used previously in the literature to detect quasicycles McKane and Newman (2005); Alonso et al. (2007); Rozhnova and Nunes (2009).

While Fourier analysis of a signal is a natural tool for studying oscillatory behaviour, a corresponding time-domain analysis must yield equivalent results. The time-correlation function forms the basis of a time-domain analysis, which for the multivariate Ornstein-Uhlenbeck process is given by Equation (20). The temporal variation of the time correlation is fixed entirely by the drift which is the deterministic part of the dynamics, while its scale is set by which involves the stochastic part of the dynamics through . Defining a normalised time correlation , we find that . This is of the form . This observation motivates the use of the ratio to reliably detect quasicycles within the linear noise approximation, where . If the decay time scale, fixed by Re(), is too short compared to the oscillatory time scale fixed by Im(), the decay will dominate and oscillatory effects will not be discernible. This will be so even when the extremum condition has real roots. We thus propose a condition for clearly discernible quasicycles, namely . In Figure (2) we plot the contour . The region is bounded on the right by this contour. As this is more stringent than the extremum condition , it is entirely contained by the regions where the PSD has a peak. In Figure (5) we emphasise this point by comparing the ACF when the PSD has peaks at finite frequencies. When is small the oscillations are barely discernible as seen from the rapid decay of the ACF. For of order unity clear signatures of oscillation are visible in the ACF.

## V Quality of noise-induced oscillations: Stochastic coherence

Noise-induced oscillations, unlike genuine oscillations, are not phase coherent and as such are called quasicycles. The coherence or regularity of quasicycles can be quantified by several measures. Here we use the quality factor, which measures the sharpness of the peak of the PSD, and the coefficient of variation which measures the regularity of the zero crossing of the signals.

The quality factor, , is a dimensionless parameter that characterizes an oscillator’s bandwidth relative to its peak frequency,

(23) |

where is the peak frequency and is the bandwidth. A high corresponds to oscillations of greater regularity. We calculate the for the diagonal entries of the PSD matrix. Let be half the maximal power for each . We calculate the bandwidth or the full-width at half-maximum (FWHM) using the and Equation (22) to get

(24) |

Using Equation (22), the peak frequencies are given by the positive square roots of the positive roots of the two quadratic equations (for ) where and , and are as defined in the previous section. The peak frequency and the FWHM together give the . Figure 6(a) shows a scan of the quality factor against population size and against the inverse of population size (inset). As one would expect, is low for high noise amplitudes and starts increasing as the noise decreases, keeping in mind that the relative noise amplitude varies as the inverse of the square root of the size of the population. However, the graph then has a maximum and then decreases for high amplitudes of noise. This is stochastic coherence.

The coefficient of variation, , is the variance over the mean of the times between succesive zeros of a temporal signal. A sharp peak in the histogram of the intervals between zero crossings, then, indicates a strongly coherent signal. is a dimensionless measure of this,

(25) |

A low indicates a high degree of coherence in the signal. Similar measures are used in the literature (for example Pikovsky and Kürths (1997) and Hilborn and Erwin (2008)). Figure 6(b) shows the for the mean crossing time of the numerical signal of the susceptibles scanned against population size and (inset) against its inverse . The plot shows a minimum which indicates stochastic coherence and hence numerically supports the analytical result given by the .

Although this non-intuitive variation of the coherence of the quasicyles with the size of the population has a stochastic origin, it is controlled purely by the deterministic part of the dynamics. The analysis using the and the require a knowledge of the diffusion matrix . However, after the two-step linearisation procedure, the entire non-trivial dependance on the population size is contained only in the spectrum of the linearised drift matrix, while the diffusion matrix scales linearly with , as given by Equations (13) and (14). Thus, any non-monotonicities in the fluctuations arise purely from the deterministic part of the dynamics, while the noise merely excites these modes. For a system which can be reduced to a standard multivariate Ornstein-Uhlenbeck process, the linearised drift matrix is identical to the linearised Jacobian matrix. This motivates the use of the ratio in determining the size of the population at which stochastic coherence is observed. This allows us to study stochastic coherence from the deterministic part of the dynamics.

We observe that in Figure (7) the ratio , when scanned against the size of the population (), shows a peak which occurs at for the parameters , . We see that this value matches well with the peaks in Figures 6(a) and 6(b), within numerical errors. We have calculated the peak value of the ratio in terms of the model parameters. If is the population size at stochastic coherence, then

(26) |

Since the ratio is always positive, there is stochastic coherence for all values of parameters for which quasicycles exist.

## Vi Detailed balance violation necessary for quasicycles

Quasicyles and stochastic coherence are not possible unless detailed balance is violated in the master equation. Typically, variables characterising biological systems are even under time-reversal, , and this is true for the variables of the SIRS model. If the system is in state at time (say) and is in state at some later time , then the joint probability of the forward transition () is while that of the reverse () is provided all time-reversal parities are even. Microscopic reversibility implies that at equilibrium the steady state forward and reverse joint probabilities must be equal. This is the condition for detailed balance and for a Markov process can be written as

(27) |

where the subscript ‘s’ denotes steady state. Expressing this condition macroscopically in terms of the correlation function, expanding in Taylor series, and keeping the first order terms one obtains Lax (1960) the Onsager relations

(28) |

which is the macroscopic condition for equilibrium. This condition requires the drift matrix to be related by a similarity transformation to a symmetric matrix Gardiner (2004) and hence restricts its spectrum to the real axis. Since it is not possible to have quasicycles without having complex eigenvalues, the violation of detailed balance is a necessary condition for the existence of noise-induced oscillations.

Recalling that and using the symmetry properties of and we can write down the following expression for Tomita and Tomita (1974)

(29) |

which then is a measure of the deviation from detailed balance. The SIRS matrix (Equation 17) can never be zero for any choice of parameters under the endemic condition . Thus the SIRS model always violates detailed balance and therefore allows for quasicycles for any choice of parameters.

## Vii Non-normality increases Variance

We have already noted that the violation of detailed balance is necessary for quasicycles. Here we further note that detailed balance violation has another consequence, that of enhancement of fluctuation amplitudes. With detailed balance the drift matrix is similar to a symmetric matrix, and is therefore normal (). In the absence of detailed balance, the drift matrix is no longer symmetric, and in this case is also non-normal.

As has been noted by Ioannou Ioannou (1995), the variance of a non-normal system driven by diagonal white noise is larger than its normal counterpart. Consider two stationary multivariate Ornstein-Uhlenbeck processes with drift and diffusion matrices () and () where is non-normal but shares the same eigenvalues as the normal . Then, Schur decompositions of the two matrices gives and where is unitary, is diagonal matrix of eigenvalues and is strictly upper triangular. Restricting the forcing to be diagonally correlated white noise (), Ioannou shows that , where and are the respective covariance matrices and denotes the trace of a matrix. For a general which is not necessarily diagonal, and . We have calculated the ratio of the traces of and .

(30) |

where is the determinant of . This expression is valid only when the spectrum of is purely real. For the SIRS model, this ratio is greater than unity.

Individual time-traces also show an increase in variance. Figure (8) shows time-traces of S and I where the fluctuations are seen to be higher than the expected standard deviation values (, where is the mean) marked by the black lines.

## Viii Conclusion

In this paper, we have analysed a closed endemic model for sustained, though asymptotically incoherent, oscillations in the population classes. These oscillations are generated through fluctuations brought about by internal demographic stochasticity which destabilise the endemic fixed point. The closed nature of the problem allows one to deal with a simplified lower-dimensional problem, an aspect we have exploited systematically by showing how the master equation can be marginalised using the constraint. This model also lends itself to a two-stage linearisation procedure, at the end of which it is reduced to a multivariate Ornstein-Uhlenbeck form. This results in the identification of the linearised drift matrix with the deterministic Jacobian matrix linearised about the endemic fixed point and permits the analysis of stochastic behaviour from the deterministic behaviour.

Noise-induced oscillations or quasicycles are possible only if the eigenvalues of the linearised Jacobian matrix are complex. These oscillations are distinct from those produced by external periodic agencies because their phases decorrelate asymptotically. Quasicycles can be reliably detected only if the oscillation time period is at least of the same order as the decorrelation time scale, as otherwise the decay dominates over the oscillation. Strong quasicycles are seen when the imaginary parts of the eigenvalues are larger than the real parts.

Stochastic coherence, or the non-trivial maximisation of regularity of the oscillations at intermediate relative noise amplitudes (or equivalently at intermediate population sizes), is a striking aspect of the SIRS quasicycles. We have seen this both analytically from the relative strength of the peak of the power spectral density and numerically by directly computing the signal-to-noise ratio of the time-traces of each population class. This analysis requires a knowledge of the intrinsic noise in the system, namely the diffusion matrix . However, we find that, for systems which can be reduced to a standard multivariate Ornstein-Uhlenbeck form by the two-stage linearisation procedure mentioned earlier, it is possible to predict stochastic coherence purely from the deterministic analysis. Any non-trivial dependance on population size is contained only in the eigenvalues of the linearised drift matrix or equivalently the linearised Jacobian matrix, while the diffusion matrix scales linearly with the population size. Thus, any non-monotonicities in fluctuations arise entirely from the deterministic part of the dynamics, i.e. the spectrum of the drift matrix. The noise merely excites these modes. This motivates the maximisation of the ratio in the investigation of the population size value at which stochastic coherence is observed. Numerical results support this procedure. Therefore, we conclude that it is possible to make predictions about non-trivial behaviour of such systems in the stochastic regime by simply analysing the linearised deterministic dynamics.

The violation of detailed balance is a necessary condition for the existence of quasicycles. Analysis of the drift, diffusion and matrices indicates that the population system described by the SIRS model is always out of equilibrium and allows for quasicycles about the endemic fixed point for any choice of model parameters. Violation of detailed balance due to the non-normal nature of the system dynamics is manifest in the enhancement of fluctuation amplitudes of the populations. We have given an expression for the ratio of the trace of the non-normal covariance matrix over its normal counterpart, restricted to parameter values where the Jacobian spectrum is purely real. Numerics indicate that this ratio is greater than unity for the SIRS model.

The analysis of this paper shows that the phenomenon of noise-induced oscillations and stochastic coherence can generically be expected in non-equilibrium birth-death jump Markov processes which can be reduced to the standard multivariate Ornstein-Uhlenbeck form by a successive application of two linearisation procedures: the linear noise approximation followed by a linearisation about the fixed point of the system. This may therefore explain the appearance of asymptotically incoherent oscilations in other systems described by such equations, as for instance, in the repressilator Ghose et al. (2010).

###### Acknowledgements.

We thank Prof. Rama Govindarajan for bringing non-normality to our attention and Sandeep K Goyal for discussions. We thank Profs. I. Bose and G. I. Menon for a critical reading of the manuscript. We thank PRISM, DAE for funding.## References

- Bernoulli (1760) D. Bernoulli, Histoires et Mémoires de l’Académie Royale des Sciences de Paris pp. 1–45 (1760).
- Kermack and McKendrick (1927) W. O. Kermack and A. G. McKendrick, Proc. R. Soc. Lond. A 115, 700 (1927).
- Kuperman and Abramson (2001) M. Kuperman and G. Abramson, Phys. Rev. Lett. 86, 2909 (2001).
- Pastor-Satorras and Vespignani (2001) R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001).
- Keeling (2005) M. Keeling, Theor. Popul. Biol. 67, 1 (2005).
- Raggett (1982) G. F. Raggett, Bull. Inst. Math. and its Applic. 18, 221 (1982).
- Murray (2002) J. D. Murray, Mathematical Biology: An introduction, vol. 1 (Springer, 2002), 3rd ed.
- Bartlett (1956) M. S. Bartlett, Proc. Third Berkeley Symp. on Math. Statist. and Prob. 4, 81 (1956).
- Bartlett (1957) M. S. Bartlett, J. R. Stat. Soc. A 120, 48 (1957).
- Soper (1929) H. E. Soper, J. R. Stat. Soc. 92, 34 (1929).
- Wilson and Worcester (1945) E. B. Wilson and J. Worcester, Proc. Natl. Acad. Sci. 31, 294 (1945).
- Nisbet and Gurney (1976) R. M. Nisbet and W. S. C. Gurney, Nature 263, 319 (1976).
- Schrödinger (1992) E. Schrödinger, What is life? (Cambridge University Press, 1992).
- McKane and Newman (2005) A. J. McKane and T. J. Newman, Phys. Rev. Lett. 94, 218102 (2005).
- van Kampen (2007) N. G. van Kampen, Stochastic processes in physics and chemistry (Elsevier, 2007), 3rd ed.
- Alonso et al. (2007) D. Alonso, A. J. McKane, and M. Pascual, J. R. Soc. Interface 4, 575 (2007).
- Rozhnova and Nunes (2009) G. Rozhnova and A. Nunes, Phys. Rev. E 79, 041922 (2009).
- Nisbet and Gurney (1982) R. M. Nisbet and W. S. C. Gurney, Modelling fluctuating populations (Wiley, 1982).
- Pikovsky and Kürths (1997) A. S. Pikovsky and J. Kürths, Phys. Rev. Lett. 78, 775 (1997).
- Hilborn and Erwin (2008) R. C. Hilborn and J. D. Erwin, J. Theor. Biol. 253, 349 (2008).
- Giacomelli et al. (2000) G. Giacomelli, M. Giudici, S. Balle, and J. R. Tredicce, Phys. Rev. Lett. 84, 3298 (2000).
- Gardiner (2004) C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences (Springer, 2004), 3rd ed.
- Doob (1945) J. L. Doob, Transactions of the American Mathematical Society pp. 455–473 (1945).
- Gillespie (1976) D. T. Gillespie, J. Comput. Phys. 22, 403 (1976).
- Gillespie (1977) D. T. Gillespie, J. Phys. Chem 81, 2340 (1977).
- Barnett and Storey (1966) S. Barnett and C. Storey, Electron. Lett. 2, 165 (1966).
- Lax (1960) M. Lax, Rev. Mod. Phys. 32, 25 (1960).
- Tomita and Tomita (1974) K. Tomita and H. Tomita, Prog. Theor. Phys. 51, 1731 (1974).
- Ioannou (1995) P. J. Ioannou, J. Atmos. Sci. 52, 1155 (1995).
- Ghose et al. (2010) S. Ghose, V. Umaiyal, J. Abharna, and R. Adhikari (2010), in preparation.