Numerical integration of the master equation in some models of stochastic epidemiology

The processes by which disease spreads in a population of individuals are inherently stochastic. The master equation has proven to be a useful tool for modeling such processes. Unfortunately, solving the master equation analytically is possible only in limited cases (e.g., when the model is linear), and thus numerical procedures or approximation methods must be employed. Available approximation methods, such as the system size expansion method of van Kampen, may fail to provide reliable solutions, whereas current numerical approaches can induce appreciable computational cost. In this paper, we propose a new numerical technique for solving the master equation. Our method is based on a more informative stochastic process than the population process commonly used in the literature. By exploiting the structure of the master equation governing this process, we develop a novel technique for calculating the exact solution of the master equation--up to a desired precision--in certain models of stochastic epidemiology. We demonstrate the potential of our method by solving the master equation associated with the stochastic SIR epidemic model. MATLAB software that implements the methods discussed in this paper is freely available as Supporting Information S1.     

Weak convergence of a sequence of stochastic difference equations to a stochastic ordinary differential equation

We consider a sequence of discrete parameter stochastic processes defined by solutions to stochastic difference equations. A condition is given that this sequence converges weakly to a continuous parameter process defined by solutions to a stochastic ordinary differential equation. Applying this result, two limit theorems related to population biology are proved. Random parameters in stochastic difference equations are autocorrelated stationary Gaussian processes in the first case. They are jump-type Markov processes in the second case. We discuss a problem of continuous time approximations for discrete time models in random environments.     

Bayesian inference for stochastic kinetic models using a diffusion approximation

This article is concerned with the Bayesian estimation of stochastic rate constants in the context of dynamic models of intracellular processes. The underlying discrete stochastic kinetic model is replaced by a diffusion approximation (or stochastic differential equation approach) where a white noise term models stochastic behavior and the model is identified using equispaced time course data. The estimation framework involves the introduction of m- 1 latent data points between every pair of observations. MCMC methods are then used to sample the posterior distribution of the latent process and the model parameters. The methodology is applied to the estimation of parameters in a prokaryotic autoregulatory gene network.     

A study of some diffusion models of population growth

The stochastic differential equations of many diffusion processes which arise in studies of population growth in random environments can be transformed, if the Stratonovich stochastic calculus is employed, to the equation of the Wiener process. If the transformation function has certain properties then the transition probability density function and quantities relating to the time to first attain a given population size can be obtained from the known results for the Wiener process. Some other random growth processes can be derived from the Ornstein-Uhlenbeck process. These transformation methods are applied to the random processes of Malthusian growth, Pearl-Verhulst logistic growth and a recent model of density independent growth due to Levins.     

A note on the diffusion approximation for single neuron firing problems

Summary This paper presents a well-known stochastic model used to describe the firing or discharge pattern of a single neuron in terms of various input processes, and shows how the potential level of the neuron can be given by means of a diffusion equation approximation. There is a discussion of the adequacy of this approximation, and the paper concludes with a brief discussion of first passage time problems.Supported in part by a grant from the Alfred P. Sloan Foundation to the Committee on Mathematical Biology and by a grant of the Statistical Branch, Office of Naval Research to the Department of Statistics, University of Chicago.     

Dynamic maximum entropy provides accurate approximation of structured population dynamics

Realistic models of biological processes typically involve interacting components on multiple scales, driven by changing environment and inherent stochasticity. Such models are often analytically and numerically intractable. We revisit a dynamic maximum entropy method that combines a static maximum entropy with a quasi-stationary approximation. This allows us to reduce stochastic non-equilibrium dynamics expressed by the Fokker-Planck equation to a simpler low-dimensional deterministic dynamics, without the need to track microscopic details. Although the method has been previously applied to a few (rather complicated) applications in population genetics, our main goal here is to explain and to better understand how the method works. We demonstrate the usefulness of the method for two widely studied stochastic problems, highlighting its accuracy in capturing important macroscopic quantities even in rapidly changing non-stationary conditions. For the Ornstein-Uhlenbeck process, the method recovers the exact dynamics whilst for a stochastic island model with migration from other habitats, the approximation retains high macroscopic accuracy under a wide range of scenarios in a dynamic environment.     

The Evolution of Restricted Recombination and the Accumulation of Repeated DNA Sequences

We suggest hypotheses to account for two major features of chromosomal organization in higher eukaryotes. The first of these is the general restriction of crossing over in the neighborhood of centromeres and telomeres. We propose that this is a consequence of selection for reduced rates of unequal exchange between repeated DNA sequences for which the copy number is subject to stabilizing selection: microtubule binding sites, in the case of centromeres, and the short repeated sequences needed for terminal replication of a linear DNA molecule, in the case of telomeres. An association between proximal crossing over and nondisjunction would also favor the restriction of crossing over near the centromere. The second feature is the association between highly repeated DNA sequences of no obvious functional significance and regions of restricted crossing over. We show that highly repeated sequences are likely to persist longest (over evolutionary time) when crossing over is infrequent. This is because unequal exchange among repeated sequences generates single copy sequences, and a population that becomes fixed for a single copy sequence by drift remains in this state indefinitely (in the absence of gene amplification processes). Increased rates of exchange thus speed up the process of stochastic loss of repeated sequences.     

Semiclassical approximations of stochastic epidemiological processes towards parameter estimation using as prime example the SIS system with import

In this paper we investigate several schemes to approximate the stationary distribution of the stochastic SIS system with import. We begin by presenting the model and analytically computing its stationary distribution. We then approximate this distribution using Kramers–Moyal approximation, van Kampen's system size expansion, and a semiclassical scheme, also called WKB or eikonal approximation depending on its different applications in physics. For the semiclassical scheme, done in the context of the Hamilton–Jacobi formalism, two approaches are taken. In the first approach we assume a semiclassical ansatz for the generating function, while in the second the solution of the master equation is approximated directly. The different schemes are compared and the semiclassical approximation, which performs better, is then used to analyse the time dependent solution of stochastic systems for which no analytical expression is known. Stochastic epidemiological models are studied in order to investigate how far such semiclassical approximations can be used for parameter estimation.     

Sustained oscillations via coherence resonance in SIR

Sustained oscillations in a stochastic SIR model are studied using a new multiple scale analysis. It captures the interaction of the deterministic and stochastic elements together with the separation of time scales inherent in the appearance of these dynamics. The nearly regular fluctuations in the infected and susceptible populations are described via an explicit construction of a stochastic amplitude equation. The agreement between the power spectral densities of the full model and the approximation verifies that coherence resonance is driving the behavior. The validity criteria for this asymptotic approximation give explicit expressions for the parameter ranges in which one expects to observe this phenomenon.     

Models of Darwinian processes and evolutionary principles

Evolutionary processes are described as stochastic motions in a genotype space (set of sequences with a Hamming distance) and a phenotype space (vector space of phenotypic properties). Real value functions are introduced which form a landscape over these spaces; smoothness postulates are formulated. Evolution is considered as a kind of hill climbing on these adaptive landscapes. A rather simple diffusion approximation for the phenotypic processes is proposed which leads to similar mathematical problems as the Schrödinger equation for disordered potential distributions.     

A Jump-Growth Model for Predator–Prey Dynamics: Derivation and Application to Marine Ecosystems

This paper investigates the dynamics of biomass in a marine ecosystem. A stochastic process is defined in which organisms undergo jumps in body size as they catch and eat smaller organisms. Using a systematic expansion of the master equation, we derive a deterministic equation for the macroscopic dynamics, which we call the deterministic jump-growth equation, and a linear Fokker–Planck equation for the stochastic fluctuations. The McKendrick–von Foerster equation, used in previous studies, is shown to be a first-order approximation, appropriate in equilibrium systems where predators are much larger than their prey. The model has a power-law steady state consistent with the approximate constancy of mass density in logarithmic intervals of body mass often observed in marine ecosystems. The behaviours of the stochastic process, the deterministic jump-growth equation, and the McKendrick–von Foerster equation are compared using numerical methods. The numerical analysis shows two classes of attractors: steady states and travelling waves.     

Modeling stochastic noise in gene regulatory systems

The Master equation is considered the gold standard for modeling the stochastic mechanisms of gene regulation in molecular detail, but it is too complex to solve exactly in most cases, so approximation and simulation methods are essential. However, there is still a lack of consensus about the best way to carry these out. To help clarify the situation, we review Master equation models of gene regulation, theoretical approximations based on an expansion method due to N.G. van Kampen and R. Kubo, and simulation algorithms due to D.T. Gillespie and P. Langevin. Expansion of the Master equation shows that for systems with a single stable steady-state, the stochastic model reduces to a deterministic model in a first-order approximation. Additional theory, also due to van Kampen, describes the asymptotic behavior of multistable systems. To support and illustrate the theory and provide further insight into the complex behavior of multistable systems, we perform a detailed simulation study comparing the various approximation and simulation methods applied to synthetic gene regulatory systems with various qualitative characteristics. The simulation studies show that for large stochastic systems with a single steady-state, deterministic models are quite accurate, since the probability distribution of the solution has a single peak tracking the deterministic trajectory whose variance is inversely proportional to the system size. In multistable stochastic systems, large fluctuations can cause individual trajectories to escape from the domain of attraction of one steady-state and be attracted to another, so the system eventually reaches a multimodal probability distribution in which all stable steady-states are represented proportional to their relative stability. However, since the escape time scales exponentially with system size, this process can take a very long time in large systems.     

Weak convergence of discrete time non-markovian processes related to selection models in population genetics

We consider discrete time stochastic processes defined by solutions to some non-linear difference equations whose coefficients are autocorrelated random sequences. It is proved that these processes converge weakly in D[0, T] to diffusion processes, under the assumption that the random sequences satisfy some mixing condition. Diffusion approximation for stochastic selection models in population genetics is discussed, as the application of this limit theorem.     

Formulas for intrinsic noise evaluation in oscillatory genetic networks

The linear noise approximation is a useful method for stochastic noise evaluations in genetic regulatory networks, where the covariance equation described as a Lyapunov equation plays a central role. We discuss the linear noise approximation method for evaluations of an intrinsic noise in autonomously oscillatory genetic networks; in such oscillatory networks, the covariance equation becomes a periodic differential equation that provides generally an unbounded covariance matrix, so that the standard method of noise evaluation based on the covariance matrix cannot be adopted directly. In this paper, we develop a new method of noise evaluation in oscillatory genetic networks; first, we investigate structural properties, e.g., orbital stability and periodicity, of the solutions to the covariance equation given as a periodic Lyapunov differential equation by using the Floquet-Lyapunov theory, and propose a global measure for evaluating stochastic amplitude fluctuations on the periodic trajectory; we also derive an evaluation formula for the period fluctuation. Finally, we apply our method to a model of circadian oscillations based on negative auto-regulation of gene expression, and show validity of our method by comparing the evaluation results with stochastic simulations.     

First passage time problem for a drifted Ornstein-Uhlenbeck process

We consider a continuous stochastic process defined as a drifted Ornstein-Uhlenbeck, for which the first passage time is of interest. The process being non-homogeneous, the first passage time probability density function cannot be found analytically, but numerical methods enable to find its estimate. Estimating the first passage time implies solving an unsteady convection-diffusion equation, with variable coefficients, and we use an implicit Euler scheme to solve it. This work is applied to simulated data, and the continuous process is inspired from recent work on biological marker modelling for HIV-positive patients. The first passage time probability density function can be useful to compare the marker progression in different groups. Numerical results show that the first passage time is highly dependent from the process perturbation, and is then more relevant than methods not considering the stochastic process directly to compare the progression.     

Logistic population growth under random dispersal

Various diffusion processes employed for modelling logistic growth are briefly summarized. A discrete-time, discrete-state space stochastic process for population growth is proposed and analyzed with either Bose-Einstein or Maxwell-Boltzmann statistics for the distribution of offspring in available sites in a restricted region. A diffusion approximation is constructed, which differs from those previously employed. The logistic law is a natural deterministic analog of the diffusion process.     

Rate of excitation propagation in a reduced Hodgkins-Huxley model. III. Integrodifferential equations]

The Hodgkin - Huxley system of equations is reduced to single integral-differential equation in neglection of slow variables dynamics. Two limiting cases of fast and slow sodium activation processes are considered. The first case leads to a nonlinear differential equation for the potential, the second one - to an ordinary differential equation with a known source as a function of coordinate. Such a simplification is due to approximation of steady-state sodium activation variable with the help of Heviside function. The validity of this approximation is discussed; the corresponding error is estimated by calculation of the second approximation for the source function.     

Investigating the two-moment characterisation of subcellular biochemical networks

While ordinary differential equations (ODEs) form the conceptual framework for modelling many cellular processes, specific situations demand stochastic models to capture the influence of noise. The most common formulation of stochastic models for biochemical networks is the chemical master equation (CME). While stochastic simulations are a practical way to realise the CME, analytical approximations offer more insight into the influence of noise. Towards that end, the two-moment approximation (2MA) is a promising addition to the established analytical approaches including the chemical Langevin equation (CLE) and the related linear noise approximation (LNA). The 2MA approach directly tracks the mean and (co)variance which are coupled in general. This coupling is not obvious in CME and CLE and ignored by LNA and conventional ODE models. We extend previous derivations of 2MA by allowing (a) non-elementary reactions and (b) relative concentrations. Often, several elementary reactions are approximated by a single step. Furthermore, practical situations often require the use of relative concentrations. We investigate the applicability of the 2MA approach to the well-established fission yeast cell cycle model. Our analytical model reproduces the clustering of cycle times observed in experiments. This is explained through multiple resettings of M-phase promoting factor (MPF), caused by the coupling between mean and (co)variance, near the G2/M transition.     

Stochastic epidemics: major outbreaks and the duration of the endemic period

A study is made of a two-dimensional stochastic system that models the spread of an infectious disease in a population. An asymptotic expression is derived for the probability that a major outbreak of the disease will occur in case the number of infectives is small. For the case that a major outbreak has occurred, an asymptotic approximation is derived for the expected time that the disease is in the population. The analytical expressions are obtained by asymptotically solving Dirichlet problems based on the Fokker-Planck equation for the stochastic system. Results of numerical calculations for the analytical expressions are compared with simulation results.