A new method for choosing the computational cell in stochastic reaction-diffusion systems

Multi-Grid Reaction-Diffusion Master Equation: Applications to Morphogen Gradient Modelling

The multi-grid reaction-diffusion master equation (mgRDME) provides a generalization of stochastic compartment-based reaction-diffusion modelling described by the standard reaction-diffusion master equation (RDME). By enabling different resolutions on lattices for biochemical species with different diffusion constants, the mgRDME approach improves both accuracy and efficiency of compartment-based reaction-diffusion simulations. The mgRDME framework is examined through its application to morphogen gradient formation in stochastic reaction-diffusion scenarios, using both an analytically tractable first-order reaction network and a model with a second-order reaction. The results obtained by the mgRDME modelling are compared with the standard RDME model and with the (more detailed) particle-based Brownian dynamics simulations. The dependence of error and numerical cost on the compartment sizes is defined and investigated through a multi-objective optimization problem.

A Random Walk Approach to Transport in Tissues and Complex Media: From Microscale Descriptions to Macroscale Models

The biological processes necessary for the development and continued survival of any organism are often strongly influenced by the transport properties of various biologically active species. The transport phenomena involved vary over multiple temporal and spatial scales, from organism-level behaviors such as the search for food, to systemic processes such as the transport of oxygen from the lungs to distant organs, down to microscopic phenomena such as the stochastic movement of proteins in a cell. Each of these processes is influenced by many interrelated factors. Identifying which factors are the most important, and how they interact to determine the overall result is a problem of great importance and interest. Experimental observations are often fit to relatively simple models, but in reality the observations are the output of complicated functions of the physicochemical, topological, and geometrical properties of a given system. Herein we use multistate continuous-time random walks and generalized master equations to model transport processes involving spatial jumps, immobilization at defined sites, and stochastic internal state changes. The underlying spatial models, which are framed as graphs, may have different classes of nodes, and walkers may have internal states that are governed by a Markov process. A general form of the solutions, using Fourier-Laplace transforms and asymptotic analysis, is developed for several spatially infinite regular lattices in one and two spatial dimensions, and the theory is developed for the analysis of transport and internal state changes on general graphs. The goal in each case is to shed light on how experimentally observable macroscale transport coefficients can be explained in terms of microscale properties of the underlying processes. This work is motivated by problems arising in transport in biological tissues, but the results are applicable to a broad class of problems that arise in other applications.

Multiscale Stochastic Reaction-Diffusion Algorithms Combining Markov Chain Models with Stochastic Partial Differential Equations

Two multiscale algorithms for stochastic simulations of reaction-diffusion processes are analysed. They are applicable to systems which include regions with significantly different concentrations of molecules. In both methods, a domain of interest is divided into two subsets where continuous-time Markov chain models and stochastic partial differential equations (SPDEs) are used, respectively. In the first algorithm, Markov chain (compartment-based) models are coupled with reaction-diffusion SPDEs by considering a pseudo-compartment (also called an overlap or handshaking region) in the SPDE part of the computational domain right next to the interface. In the second algorithm, no overlap region is used. Further extensions of both schemes are presented, including the case of an adaptively chosen boundary between different modelling approaches.

Spatial Stochastic Intracellular Kinetics: A Review of Modelling Approaches

Models of chemical kinetics that incorporate both stochasticity and diffusion are an increasingly common tool for studying biology. The variety of competing models is vast, but two stand out by virtue of their popularity: the reaction-diffusion master equation and Brownian dynamics. In this review, we critically address a number of open questions surrounding these models: How can they be justified physically? How do they relate to each other? How do they fit into the wider landscape of chemical models, ranging from the rate equations to molecular dynamics? This review assumes no prior knowledge of modelling chemical kinetics and should be accessible to a wide range of readers.

Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach

We develop numerical methods for stochastic reaction-diffusion systems based on approaches used for fluctuatinghydrodynamics (FHD). For hydrodynamicsystems, the FHD formulation is formally described by stochastic partial differential equations (SPDEs). In the reaction-diffusion systems we consider, our model becomes similar to the reaction-diffusion master equation (RDME) description when our SPDEs are spatially discretized and reactions are modeled as a source term having Poissonfluctuations. However, unlike the RDME, which becomes prohibitively expensive for an increasing number of molecules, our FHD-based description naturally extends from the regime where fluctuations are strong, i.e., each mesoscopic cell has few (reactive) molecules, to regimes with moderate or weak fluctuations, and ultimately to the deterministic limit. By treating diffusion implicitly, we avoid the severe restriction on time step size that limits all methods based on explicit treatments of diffusion and construct numerical methods that are more efficient than RDME methods, without compromising accuracy. Guided by an analysis of the accuracy of the distribution of steady-state fluctuations for the linearized reaction-diffusion model, we construct several two-stage (predictor-corrector) schemes, where diffusion is treated using a stochastic Crank-Nicolson method, and reactions are handled by the stochastic simulation algorithm of Gillespie or a weakly second-order tau leaping method. We find that an implicit midpoint tau leaping scheme attains second-order weak accuracy in the linearized setting and gives an accurate and stable structure factor for a time step size of an order of magnitude larger than the hopping time scale of diffusing molecules. We study the numerical accuracy of our methods for the Schlögl reaction-diffusion model both in and out of thermodynamic equilibrium. We demonstrate and quantify the importance of thermodynamicfluctuations to the formation of a two-dimensional Turing-like pattern and examine the effect of fluctuations on three-dimensional chemical front propagation. By comparing stochastic simulations to deterministic reaction-diffusion simulations, we show that fluctuations accelerate pattern formation in spatially homogeneous systems and lead to a qualitatively different disordered pattern behind a traveling wave.

A hybrid continuous-discrete method for stochastic reaction-diffusion processes

Stochastic fluctuations in reaction-diffusion processes often have substantial effect on spatial and temporal dynamics of signal transductions in complex biological systems. One popular approach for simulating these processes is to divide the system into small spatial compartments assuming that molecules react only within the same compartment and jump between adjacent compartments driven by the diffusion. While the approach is convenient in terms of its implementation, its computational cost may become prohibitive when diffusive jumps occur significantly more frequently than reactions, as in the case of rapid diffusion. Here, we present a hybrid continuous-discrete method in which diffusion is simulated using continuous approximation while reactions are based on the Gillespie algorithm. Specifically, the diffusive jumps are approximated as continuous Gaussian random vectors with time-dependent means and covariances, allowing use of a large time step, even for rapid diffusion. By considering the correlation among diffusive jumps, the approximation is accurate for the second moment of the diffusion process. In addition, a criterion is obtained for identifying the region in which such diffusion approximation is required to enable adaptive calculations for better accuracy. Applications to a linear diffusion system and two nonlinear systems of morphogens demonstrate the effectiveness and benefits of the new hybrid method.

How spatial heterogeneity shapes multiscale biochemical reaction network dynamics

Spatial heterogeneity in cells can be modelled using distinct compartments connected by molecular movement between them. In addition to movement, changes in the amount of molecules are due to biochemical reactions within compartments, often such that some molecular types fluctuate on a slower timescale than others. It is natural to ask the following questions: how sensitive is the dynamics of molecular types to their own spatial distribution, and how sensitive are they to the distribution of others? What conditions lead to effective homogeneity in biochemical dynamics despite heterogeneity in molecular distribution? What kind of spatial distribution is optimal from the point of view of some downstream product? Within a spatially heterogeneous multiscale model, we consider two notions of dynamical homogeneity (full homogeneity and homogeneity for the fast subsystem), and consider their implications under different timescales for the motility of molecules between compartments. We derive rigorous results for their dynamics and long-term behaviour, and illustrate them with examples of a shared pathway, Michaelis-Menten enzymatic kinetics and autoregulating feedbacks. Using stochastic averaging of fast fluctuations to their quasi-steady-state distribution, we obtain simple analytic results that significantly reduce the complexity and expedite simulation of stochastic compartment models of chemical reactions.

Lyapunov Functions, Stationary Distributions, and Non-equilibrium Potential for Reaction Networks

We consider the relationship between stationary distributions for stochastic models of reaction systems and Lyapunov functions for their deterministic counterparts. Specifically, we derive the well-known Lyapunov function of reaction network theory as a scaling limit of the non-equilibrium potential of the stationary distribution of stochastically modeled complex balanced systems. We extend this result to general birth-death models and demonstrate via example that similar scaling limits can yield Lyapunov functions even for models that are not complex or detailed balanced, and may even have multiple equilibria.

Constant-complexity stochastic simulation algorithm with optimal binning

At the molecular level, biochemical processes are governed by random interactions between reactant molecules, and the dynamics of such systems are inherently stochastic. When the copy numbers of reactants are large, a deterministic description is adequate, but when they are small, such systems are often modeled as continuous-time Markov jump processes that can be described by the chemical master equation. Gillespie's Stochastic Simulation Algorithm (SSA) generates exact trajectories of these systems, but the amount of computational work required for each step of the original SSA is proportional to the number of reaction channels, leading to computational complexity that scales linearly with the problem size. The original SSA is therefore inefficient for large problems, which has prompted the development of several alternative formulations with improved scaling properties. We describe an exact SSA that uses a table data structure with event time binning to achieve constant computational complexity with respect to the number of reaction channels for weakly coupled reaction networks. We present a novel adaptive binning strategy and discuss optimal algorithm parameters. We compare the computational efficiency of the algorithm to existing methods and demonstrate excellent scaling for large problems. This method is well suited for generating exact trajectories of large weakly coupled models, including those that can be described by the reaction-diffusion master equation that arises from spatially discretized reaction-diffusion processes.

Stochastic Turing patterns: analysis of compartment-based approaches

Turing patterns can be observed in reaction-diffusion systems where chemical species have different diffusion constants. In recent years, several studies investigated the effects of noise on Turing patterns and showed that the parameter regimes, for which stochastic Turing patterns are observed, can be larger than the parameter regimes predicted by deterministic models, which are written in terms of partial differential equations (PDEs) for species concentrations. A common stochastic reaction-diffusion approach is written in terms of compartment-based (lattice-based) models, where the domain of interest is divided into artificial compartments and the number of molecules in each compartment is simulated. In this paper, the dependence of stochastic Turing patterns on the compartment size is investigated. It has previously been shown (for relatively simpler systems) that a modeler should not choose compartment sizes which are too small or too large, and that the optimal compartment size depends on the diffusion constant. Taking these results into account, we propose and study a compartment-based model of Turing patterns where each chemical species is described using a different set of compartments. It is shown that the parameter regions where spatial patterns form are different from the regions obtained by classical deterministic PDE-based models, but they are also different from the results obtained for the stochastic reaction-diffusion models which use a single set of compartments for all chemical species. In particular, it is argued that some previously reported results on the effect of noise on Turing patterns in biological systems need to be reinterpreted.

Stochastic analysis of reaction-diffusion processes

Reaction and diffusion processes are used to model chemical and biological processes over a wide range of spatial and temporal scales. Several routes to the diffusion process at various levels of description in time and space are discussed and the master equation for spatially discretized systems involving reaction and diffusion is developed. We discuss an estimator for the appropriate compartment size for simulating reaction-diffusion systems and introduce a measure of fluctuations in a discretized system. We then describe a new computational algorithm for implementing a modified Gillespie method for compartmental systems in which reactions are aggregated into equivalence classes and computational cells are searched via an optimized tree structure. Finally, we discuss several examples that illustrate the issues that have to be addressed in general systems.

Perspective: Stochastic algorithms for chemical kinetics

We outline our perspective on stochastic chemical kinetics, paying particular attention to numerical simulation algorithms. We first focus on dilute, well-mixed systems, whose description using ordinary differential equations has served as the basis for traditional chemical kinetics for the past 150 years. For such systems, we review the physical and mathematical rationale for a discrete-stochastic approach, and for the approximations that need to be made in order to regain the traditional continuous-deterministic description. We next take note of some of the more promising strategies for dealing stochastically with stiff systems, rare events, and sensitivity analysis. Finally, we review some recent efforts to adapt and extend the discrete-stochastic approach to systems that are not well-mixed. In that currently developing area, we focus mainly on the strategy of subdividing the system into well-mixed subvolumes, and then simulating diffusional transfers of reactant molecules between adjacent subvolumes together with chemical reactions inside the subvolumes.

A comparison of bimolecular reaction models for stochastic reaction-diffusion systems

Stochastic reaction-diffusion models have become an important tool in studying how both noise in the chemical reaction process and the spatial movement of molecules influences the behavior of biological systems. There are two primary spatially-continuous models that have been used in recent studies: the diffusion limited reaction model of Smoluchowski, and a second approach popularized by Doi. Both models treat molecules as points undergoing Brownian motion. The former represents chemical reactions between two reactants through the use of reactive boundary conditions, with two molecules reacting instantly upon reaching a fixed separation (called the reaction-radius). The Doi model uses reaction potentials, whereby two molecules react with a fixed probability per unit time, λ, when separated by less than the reaction radius. In this work, we study the rigorous relationship between the two models. For the special case of a protein diffusing to a fixed DNA binding site, we prove that the solution to the Doi model converges to the solution of the Smoluchowski model as λ→∞, with a rigorous [Formula: see text] error bound (for any fixed ϵ>0). We investigate by numerical simulation, for biologically relevant parameter values, the difference between the solutions and associated reaction time statistics of the two models. As the reaction-radius is decreased, for sufficiently large but fixed values of λ, these differences are found to increase like the inverse of the binding radius.

External noise control in inherently stochastic biological systems

Biological systems are often subject to external noise from signal stimuli and environmental perturbations, as well as noises in the intracellular signal transduction pathway. Can different stochastic fluctuations interact to give rise to new emerging behaviors? How can a system reduce noise effects while still being capable of detecting changes in the input signal? Here, we study analytically and computationally the role of nonlinear feedback systems in controlling external noise with the presence of large internal noise. In addition to noise attenuation, we analyze derivatives of Fano factor to study systems' capability of differentiating signal inputs. We find effects of internal noise and external noise may be separated in one slow positive feedback loop system; in particular, the slow loop can decrease external noise and increase robustness of signaling with respect to fluctuations in rate constants, while maintaining the signal output specific to the input. For two feedback loops, we demonstrate that the influence of external noise mainly depends on how the fast loop responds to fluctuations in the input and the slow loop plays a limited role in determining the signal precision. Furthermore, in a dual loop system of one positive feedback and one negative feedback, a slower positive feedback always leads to better noise attenuation; in contrast, a slower negative feedback may not be more beneficial. Our results reveal interesting stochastic effects for systems containing both extrinsic and intrinsic noises, suggesting novel noise filtering strategies in inherently stochastic systems.

Noise drives sharpening of gene expression boundaries in the zebrafish hindbrain

Morphogens provide positional information for spatial patterns of gene expression during development. However, stochastic effects such as local fluctuations in morphogen concentration and noise in signal transduction make it difficult for cells to respond to their positions accurately enough to generate sharp boundaries between gene expression domains. During development of rhombomeres in the zebrafish hindbrain, the morphogen retinoic acid (RA) induces expression of hoxb1a in rhombomere 4 (r4) and krox20 in r3 and r5. Fluorescent in situ hybridization reveals rough edges around these gene expression domains, in which cells co-express hoxb1a and krox20 on either side of the boundary, and these sharpen within a few hours. Computational analysis of spatial stochastic models shows, surprisingly, that noise in hoxb1a/krox20 expression actually promotes sharpening of boundaries between adjacent segments. In particular, fluctuations in RA initially induce a rough boundary that requires noise in hoxb1a/krox20 expression to sharpen. This finding suggests a novel noise attenuation mechanism that relies on intracellular noise to induce switching and coordinate cellular decisions during developmental patterning.