Reaction rates for mesoscopic reaction-diffusion kinetics

Close agreement between deterministic versus stochastic modeling of first-passage time to vesicle fusion

Ca2+-dependent cell processes, such as neurotransmitter or endocrine vesicle fusion, are inherently stochastic due to large fluctuations in Ca2+ channel gating, Ca2+ diffusion, and Ca2+ binding to buffers and target sensors. However, previous studies revealed closer-than-expected agreement between deterministic and stochastic simulations of Ca2+ diffusion, buffering, and sensing if Ca2+ channel gating is not Ca2+ dependent. To understand this result more fully, we present a comparative study complementing previous work, focusing on Ca2+ dynamics downstream of Ca2+ channel gating. Specifically, we compare deterministic (mean-field/mass-action) and stochastic simulations of vesicle exocytosis latency, quantified by the probability density of the first-passage time (FPT) to the Ca2+-bound state of a vesicle fusion sensor, following a brief Ca2+ current pulse. We show that under physiological constraints, the discrepancy between FPT densities obtained using the two approaches remains small even if as few as ∼50 Ca2+ ions enter per single channel-vesicle release unit. Using a reduced two-compartment model for ease of analysis, we illustrate how this close agreement arises from the smallness of correlations between fluctuations of the reactant molecule numbers, despite the large magnitude of fluctuation amplitudes. This holds if all relevant reactions are heteroreaction between molecules of different species, as is the case for bimolecular Ca2+ binding to buffers and downstream sensor targets. In this case, diffusion and buffering effectively decorrelate the state of the Ca2+ sensor from local Ca2+ fluctuations. Thus, fluctuations in the Ca2+ sensor's state underlying the FPT distribution are only weakly affected by the fluctuations in the local Ca2+ concentration around its average, deterministically computable value.

A Gaussian jump process formulation of the reaction–diffusion master equation enables faster exact stochastic simulations

We propose a Gaussian jump process model on a regular Cartesian lattice for the diffusion part of the Reaction–Diffusion Master Equation (RDME). We derive the resulting Gaussian RDME (GRDME) formulation from analogy with a kernel-based discretization scheme for continuous diffusion processes and quantify the limits of its validity relative to the classic RDME. We then present an exact stochastic simulation algorithm for the GRDME, showing that the accuracies of GRDME and RDME are comparable, but exact simulations of the GRDME require only a fraction of the computational cost of exact RDME simulations. We analyze the origin of this speedup and its scaling with problem dimension. The benchmarks suggest that the GRDME is a particularly beneficial model for diffusion-dominated systems in three dimensional spaces, often occurring in systems biology and cell biology.

A multiscale compartment-based model of stochastic gene regulatory networks using hitting-time analysis

Spatial stochastic models of single cell kinetics are capable of capturing both fluctuations in molecular numbers and the spatial dependencies of the key steps of intracellular regulatory networks. The spatial stochastic model can be simulated both on a detailed microscopic level using particle tracking and on a mesoscopic level using the reaction-diffusion master equation. However, despite substantial progress on simulation efficiency for spatial models in the last years, the computational cost quickly becomes prohibitively expensive for tasks that require repeated simulation of thousands or millions of realizations of the model. This limits the use of spatial models in applications such as multicellular simulations, likelihood-free parameter inference, and robustness analysis. Further approximation of the spatial dynamics is needed to accelerate such computational engineering tasks. We here propose a multiscale model where a compartment-based model approximates a detailed spatial stochastic model. The compartment model is constructed via a first-exit time analysis on the spatial model, thus capturing critical spatial aspects of the fine-grained simulations, at a cost close to the simple well-mixed model. We apply the multiscale model to a canonical model of negative-feedback gene regulation, assess its accuracy over a range of parameters, and demonstrate that the approximation can yield substantial speedups for likelihood-free parameter inference.

A novel stochastic simulation approach enables exploration of mechanisms for regulating polarity site movement

Cells polarize their movement or growth toward external directional cues in many different contexts. For example, budding yeast cells grow toward potential mating partners in response to pheromone gradients. Directed growth is controlled by polarity factors that assemble into clusters at the cell membrane. The clusters assemble, disassemble, and move between different regions of the membrane before eventually forming a stable polarity site directed toward the pheromone source. Pathways that regulate clustering have been identified but the molecular mechanisms that regulate cluster mobility are not well understood. To gain insight into the contribution of chemical noise to cluster behavior we simulated clustering using the reaction-diffusion master equation (RDME) framework to account for molecular-level fluctuations. RDME simulations are a computationally efficient approximation, but their results can diverge from the underlying microscopic dynamics. We implemented novel concentration-dependent rate constants that improved the accuracy of RDME-based simulations, allowing us to efficiently investigate how cluster dynamics might be regulated. Molecular noise was effective in relocating clusters when the clusters contained low numbers of limiting polarity factors, and when Cdc42, the central polarity regulator, exhibited short dwell times at the polarity site. Cluster stabilization occurred when abundances or binding rates were altered to either lengthen dwell times or increase the number of polarity molecules in the cluster. We validated key results using full 3D particle-based simulations. Understanding the mechanisms cells use to regulate the dynamics of polarity clusters should provide insights into how cells dynamically track external directional cues.

Cells localize polarity molecules in a small region of the plasma membrane forming a polarity cluster that directs functions such as migration, reproduction, and growth. Guided by external signals, these clusters move across the membrane allowing cells to reorient growth or motion. The polarity molecules continuously and randomly shuttle between the cluster and the cell cytosol and, as a result, the number and distribution of molecules at the cluster constantly changes. Here we present an improved stochastic simulation algorithm to investigate how such molecular-scale fluctuations induce cluster movement across the cell membrane. Unexpectedly, cluster mobility does not correlate with variations in total molecule abundance within the cluster, but rather with changes in the spatial distribution of molecules that form the cluster. Cluster motion is faster when polarity molecules are scarce and when they shuttle rapidly between the cluster and the cytosol. Our results suggest that cells control cluster mobility by regulating the abundance of polarity molecules and biochemical reactions that affect the time molecules spend at the cluster. We provide insights into how cells harness random molecular behavior to perform functions important for survival, such as detecting the direction of external signals.

Quantifying the roles of space and stochasticity in computer simulations for cell biology and cellular biochemistry

Most of the fascinating phenomena studied in cell biology emerge from interactions among highly organized multimolecular structures embedded into complex and frequently dynamic cellular morphologies. For the exploration of such systems, computer simulation has proved to be an invaluable tool, and many researchers in this field have developed sophisticated computational models for application to specific cell biological questions. However, it is often difficult to reconcile conflicting computational results that use different approaches to describe the same phenomenon. To address this issue systematically, we have defined a series of computational test cases ranging from very simple to moderately complex, varying key features of dimensionality, reaction type, reaction speed, crowding, and cell size. We then quantified how explicit spatial and/or stochastic implementations alter outcomes, even when all methods use the same reaction network, rates, and concentrations. For simple cases, we generally find minor differences in solutions of the same problem. However, we observe increasing discordance as the effects of localization, dimensionality reduction, and irreversible enzymatic reactions are combined. We discuss the strengths and limitations of commonly used computational approaches for exploring cell biological questions and provide a framework for decision making by researchers developing new models. As computational power and speed continue to increase at a remarkable rate, the dream of a fully comprehensive computational model of a living cell may be drawing closer to reality, but our analysis demonstrates that it will be crucial to evaluate the accuracy of such models critically and systematically.

The blending region hybrid framework for the simulation of stochastic reaction-diffusion processes

The simulation of stochastic reaction–diffusion systems using fine-grained representations can become computationally prohibitive when particle numbers become large. If particle numbers are sufficiently high then it may be possible to ignore stochastic fluctuations and use a more efficient coarse-grained simulation approach. Nevertheless, for multiscale systems which exhibit significant spatial variation in concentration, a coarse-grained approach may not be appropriate throughout the simulation domain. Such scenarios suggest a hybrid paradigm in which a computationally cheap, coarse-grained model is coupled to a more expensive, but more detailed fine-grained model, enabling the accurate simulation of the fine-scale dynamics at a reasonable computational cost. In this paper, in order to couple two representations of reaction–diffusion at distinct spatial scales, we allow them to overlap in a ‘blending region’. Both modelling paradigms provide a valid representation of the particle density in this region. From one end of the blending region to the other, control of the implementation of diffusion is passed from one modelling paradigm to another through the use of complementary ‘blending functions’ which scale up or down the contribution of each model to the overall diffusion. We establish the reliability of our novel hybrid paradigm by demonstrating its simulation on four exemplar reaction–diffusion scenarios.

Hierarchical algorithm for the reaction-diffusion master equation

We have developed an algorithm coupling mesoscopic simulations on different levels in a hierarchy of Cartesian meshes. Based on the multiscale nature of the chemical reactions, some molecules in the system will live on a fine-grained mesh, while others live on a coarse-grained mesh. By allowing molecules to transfer from the fine levels to the coarse levels when appropriate, we show that we can save up to three orders of magnitude of computational time compared to microscopic simulations or highly resolved mesoscopic simulations, without losing significant accuracy. We demonstrate this in several numerical examples with systems that cannot be accurately simulated with a coarse-grained mesoscopic model.

Stochastic self-tuning hybrid algorithm for reaction-diffusion systems

Many biochemical phenomena involve reactants with vastly different concentrations, some of which are amenable to continuum-level descriptions, while the others are not. We present a hybrid self-tuning algorithm to model such systems. The method combines microscopic (Brownian) dynamics for diffusion with mesoscopic (Gillespie-type) methods for reactions and remains efficient in a wide range of regimes and scenarios with large variations of concentrations. Its accuracy, robustness, and versatility are balanced by redefining propensities and optimizing the mesh size and time step. We use a bimolecular reaction to demonstrate the potential of our method in a broad spectrum of scenarios: from almost completely reaction-dominated systems to cases where reactions rarely occur or take place very slowly. The simulation results show that the number of particles present in the system does not degrade the performance of our method. This makes it an accurate and computationally efficient tool to model complex multireaction systems.

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.

Learning moment closure in reaction-diffusion systems with spatial dynamic Boltzmann distributions

Many physical systems are described by probability distributions that evolve in both time and space. Modeling these systems is often challenging due to their large state space and analytically intractable or computationally expensive dynamics. To address these problems, we study a machine-learning approach to model reduction based on the Boltzmann machine. Given the form of the reduced model Boltzmann distribution, we introduce an autonomous differential equation system for the interactions appearing in the energy function. The reduced model can treat systems in continuous space (described by continuous random variables), for which we formulate a variational learning problem using the adjoint method to determine the right-hand sides of the differential equations. This approach can be used to enforce a reduced physical model by a suitable parametrization of the differential equations. The parametrization we employ uses the basis functions from finite-element methods, which can be used to model any physical system. One application domain for such physics-informed learning algorithms is to modeling reaction-diffusion systems. We study a lattice version of the Rössler chaotic oscillator, which illustrates the accuracy of the moment closure approximation made by the method and its dimensionality reduction power.

Efficient reactive Brownian dynamics

We develop a Split Reactive Brownian Dynamics (SRBD) algorithm for particle simulations of reaction-diffusion systems based on the Doi or volume reactivity model, in which pairs of particles react with a specified Poisson rate if they are closer than a chosen reactive distance. In our Doi model, we ensure that the microscopic reaction rules for various association and dissociation reactions are consistent with detailed balance (time reversibility) at thermodynamic equilibrium. The SRBD algorithm uses Strang splitting in time to separate reaction and diffusion and solves both the diffusion-only and reaction-only subproblems exactly, even at high packing densities. To efficiently process reactions without uncontrolled approximations, SRBD employs an event-driven algorithm that processes reactions in a time-ordered sequence over the duration of the time step. A grid of cells with size larger than all of the reactive distances is used to schedule and process the reactions, but unlike traditional grid-based methods such as reaction-diffusion master equation algorithms, the results of SRBD are statistically independent of the size of the grid used to accelerate the processing of reactions. We use the SRBD algorithm to compute the effective macroscopic reaction rate for both reaction-limited and diffusion-limited irreversible association in three dimensions and compare to existing theoretical predictions at low and moderate densities. We also study long-time tails in the time correlation functions for reversible association at thermodynamic equilibrium and compare to recent theoretical predictions. Finally, we compare different particle and continuum methods on a model exhibiting a Turing-like instability and pattern formation. Our studies reinforce the common finding that microscopic mechanisms and correlations matter for diffusion-limited systems, making continuum and even mesoscopic modeling of such systems difficult or impossible. We also find that for models in which particles diffuse off lattice, such as the Doi model, reactions lead to a spurious enhancement of the effective diffusion coefficients.

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.

Reaction rates for reaction-diffusion kinetics on unstructured meshes

The reaction-diffusion master equation is a stochastic model often utilized in the study of biochemical reaction networks in living cells. It is applied when the spatial distribution of molecules is important to the dynamics of the system. A viable approach to resolve the complex geometry of cells accurately is to discretize space with an unstructured mesh. Diffusion is modeled as discrete jumps between nodes on the mesh, and the diffusion jump rates can be obtained through a discretization of the diffusion equation on the mesh. Reactions can occur when molecules occupy the same voxel. In this paper, we develop a method for computing accurate reaction rates between molecules occupying the same voxel in an unstructured mesh. For large voxels, these rates are known to be well approximated by the reaction rates derived by Collins and Kimball, but as the mesh is refined, no analytical expression for the rates exists. We reduce the problem of computing accurate reaction rates to a pure preprocessing step, depending only on the mesh and not on the model parameters, and we devise an efficient numerical scheme to estimate them to high accuracy. We show in several numerical examples that as we refine the mesh, the results obtained with the reaction-diffusion master equation approach those of a more fine-grained Smoluchowski particle-tracking model.

Bridging intracellular scales by mechanistic computational models

The impact of intracellular spatial organization beyond classical compartments on processes such as cell signaling is increasingly recognized. A quantitative, mechanistic understanding of cellular systems therefore needs to account for different scales in at least three coordinates: time, molecular abundances, and space. Mechanistic mathematical models may span all these scales, but corresponding multi-scale models need to resolve mechanistic details on small scales while maintaining computational tractability for larger ones. This typically results in models that combine different levels of description: from a microscopic representation of chemical reactions up to continuum dynamics in space and time. We highlight recent progress in bridging these model classes and outline current challenges in multi-scale models such as active transport and dynamic geometries.

Mesoscopic-microscopic spatial stochastic simulation with automatic system partitioning

The reaction-diffusion master equation (RDME) is a model that allows for efficient on-lattice simulation of spatially resolved stochastic chemical kinetics. Compared to off-lattice hard-sphere simulations with Brownian dynamics or Green's function reaction dynamics, the RDME can be orders of magnitude faster if the lattice spacing can be chosen coarse enough. However, strongly diffusion-controlled reactions mandate a very fine mesh resolution for acceptable accuracy. It is common that reactions in the same model differ in their degree of diffusion control and therefore require different degrees of mesh resolution. This renders mesoscopic simulation inefficient for systems with multiscale properties. Mesoscopic-microscopic hybrid methods address this problem by resolving the most challenging reactions with a microscale, off-lattice simulation. However, all methods to date require manual partitioning of a system, effectively limiting their usefulness as "black-box" simulation codes. In this paper, we propose a hybrid simulation algorithm with automatic system partitioning based on indirect a priori error estimates. We demonstrate the accuracy and efficiency of the method on models of diffusion-controlled networks in 3D.

Multiscale Modeling of Diffusion in a Crowded Environment

We present a multiscale approach to model diffusion in a crowded environment and its effect on the reaction rates. Diffusion in biological systems is often modeled by a discrete space jump process in order to capture the inherent noise of biological systems, which becomes important in the low copy number regime. To model diffusion in the crowded cell environment efficiently, we compute the jump rates in this mesoscopic model from local first exit times, which account for the microscopic positions of the crowding molecules, while the diffusing molecules jump on a coarser Cartesian grid. We then extract a macroscopic description from the resulting jump rates, where the excluded volume effect is modeled by a diffusion equation with space-dependent diffusion coefficient. The crowding molecules can be of arbitrary shape and size, and numerical experiments demonstrate that those factors together with the size of the diffusing molecule play a crucial role on the magnitude of the decrease in diffusive motion. When correcting the reaction rates for the altered diffusion we can show that molecular crowding either enhances or inhibits chemical reactions depending on local fluctuations of the obstacle density.

Stochastic Simulation Service: Bridging the Gap between the Computational Expert and the Biologist

We present StochSS: Stochastic Simulation as a Service, an integrated development environment for modeling and simulation of both deterministic and discrete stochastic biochemical systems in up to three dimensions. An easy to use graphical user interface enables researchers to quickly develop and simulate a biological model on a desktop or laptop, which can then be expanded to incorporate increasing levels of complexity. StochSS features state-of-the-art simulation engines. As the demand for computational power increases, StochSS can seamlessly scale computing resources in the cloud. In addition, StochSS can be deployed as a multi-user software environment where collaborators share computational resources and exchange models via a public model repository. We demonstrate the capabilities and ease of use of StochSS with an example of model development and simulation at increasing levels of complexity.

Stochastic Simulation Service: Bridging the Gap between the Computational Expert and the Biologist

We present StochSS: Stochastic Simulation as a Service, an integrated development environment for modeling and simulation of both deterministic and discrete stochastic biochemical systems in up to three dimensions. An easy to use graphical user interface enables researchers to quickly develop and simulate a biological model on a desktop or laptop, which can then be expanded to incorporate increasing levels of complexity. StochSS features state-of-the-art simulation engines. As the demand for computational power increases, StochSS can seamlessly scale computing resources in the cloud. In addition, StochSS can be deployed as a multi-user software environment where collaborators share computational resources and exchange models via a public model repository. We demonstrate the capabilities and ease of use of StochSS with an example of model development and simulation at increasing levels of complexity.

A framework for discrete stochastic simulation on 3D moving boundary domains

We have developed a method for modeling spatial stochastic biochemical reactions in complex, three-dimensional, and time-dependent domains using the reaction-diffusion master equation formalism. In particular, we look to address the fully coupled problems that arise in systems biology where the shape and mechanical properties of a cell are determined by the state of the biochemistry and vice versa. To validate our method and characterize the error involved, we compare our results for a carefully constructed test problem to those of a microscale implementation. We demonstrate the effectiveness of our method by simulating a model of polarization and shmoo formation during the mating of yeast. The method is generally applicable to problems in systems biology where biochemistry and mechanics are coupled, and spatial stochastic effects are critical.

Reaction rates for a generalized reaction-diffusion master equation

It has been established that there is an inherent limit to the accuracy of the reaction-diffusion master equation. Specifically, there exists a fundamental lower bound on the mesh size, below which the accuracy deteriorates as the mesh is refined further. In this paper we extend the standard reaction-diffusion master equation to allow molecules occupying neighboring voxels to react, in contrast to the traditional approach, in which molecules react only when occupying the same voxel. We derive reaction rates, in two dimensions as well as three dimensions, to obtain an optimal match to the more fine-grained Smoluchowski model and show in two numerical examples that the extended algorithm is accurate for a wide range of mesh sizes, allowing us to simulate systems that are intractable with the standard reaction-diffusion master equation. In addition, we show that for mesh sizes above the fundamental lower limit of the standard algorithm, the generalized algorithm reduces to the standard algorithm. We derive a lower limit for the generalized algorithm which, in both two dimensions and three dimensions, is of the order of the reaction radius of a reacting pair of molecules.

Hybrid approaches for multiple-species stochastic reaction-diffusion models

Reaction-diffusion models are used to describe systems in fields as diverse as physics, chemistry, ecology and biology. The fundamental quantities in such models are individual entities such as atoms and molecules, bacteria, cells or animals, which move and/or react in a stochastic manner. If the number of entities is large, accounting for each individual is inefficient, and often partial differential equation (PDE) models are used in which the stochastic behaviour of individuals is replaced by a description of the averaged, or mean behaviour of the system. In some situations the number of individuals is large in certain regions and small in others. In such cases, a stochastic model may be inefficient in one region, and a PDE model inaccurate in another. To overcome this problem, we develop a scheme which couples a stochastic reaction-diffusion system in one part of the domain with its mean field analogue, i.e. a discretised PDE model, in the other part of the domain. The interface in between the two domains occupies exactly one lattice site and is chosen such that the mean field description is still accurate there. In this way errors due to the flux between the domains are small. Our scheme can account for multiple dynamic interfaces separating multiple stochastic and deterministic domains, and the coupling between the domains conserves the total number of particles. The method preserves stochastic features such as extinction not observable in the mean field description, and is significantly faster to simulate on a computer than the pure stochastic model.