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 discretestochastic 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 wellmixed. 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. © 2013 AIP Publishing LLC.[http://dx

Stochastic chemical systems with diffusion are modeled with a reactiondiffusion master equation. On a macroscopic level, the governing equation is a reaction-diffusion equation for the averages of the chemical species. On a mesoscopic level, the master equation for a well stirred chemical system is combined with Brownian motion in space to obtain the reactiondiffusion master equation. The space is covered by an unstructured mesh and the diffusion coefficients on the mesoscale are obtained from a finite element discretization of the Laplace operator on the macroscale. The resulting method is a flexible hybrid algorithm in that the diffusion can be handled either on the meso-or on the macroscale level. The accuracy and the efficiency of the method are illustrated in three numerical examples inspired by molecular biology.

Cartilaginous structures are at the core of embryo growth and shaping before the bone forms. Here we report a novel principle of vertebrate cartilage growth that is based on introducing transversally-oriented clones into pre-existing cartilage. This mechanism of growth uncouples the lateral expansion of curved cartilaginous sheets from the control of cartilage thickness, a process which might be the evolutionary mechanism underlying adaptations of facial shape. In rod-shaped cartilage structures (Meckel, ribs and skeletal elements in developing limbs), the transverse integration of clonal columns determines the well-defined diameter and resulting rod-like morphology. We were able to alter cartilage shape by experimentally manipulating clonal geometries. Using in silico modeling, we discovered that anisotropic proliferation might explain cartilage bending and groove formation at the macro-scale.DOI: http://dx.doi.org/10.7554/eLife.25902.001

Stochastic modeling of reaction-diffusion kinetics has emerged as a powerful theoretical tool in the study of biochemical reaction networks. Two frequently employed models are the particle-tracking Smoluchowski framework and the on-lattice Reaction-Diffusion Master Equation (RDME) framework. As the mesh size goes from coarse to fine, the RDME initially becomes more accurate. However, recent developments have shown that it will become increasingly inaccurate compared to the Smoluchowski model as the lattice spacing becomes very fine. In this paper we give a new, general and simple argument for why the RDME breaks down. Our analysis reveals a hard limit on the voxel size for which no local RDME can agree with the Smoluchowski model.

BackgroundExperiments in silico using stochastic reaction-diffusion models have emerged as an important tool in molecular systems biology. Designing computational software for such applications poses several challenges. Firstly, realistic lattice-based modeling for biological applications requires a consistent way of handling complex geometries, including curved inner- and outer boundaries. Secondly, spatiotemporal stochastic simulations are computationally expensive due to the fast time scales of individual reaction- and diffusion events when compared to the biological phenomena of actual interest. We therefore argue that simulation software needs to be both computationally efficient, employing sophisticated algorithms, yet in the same time flexible in order to meet present and future needs of increasingly complex biological modeling.ResultsWe have developed URDME, a flexible software framework for general stochastic reaction-transport modeling and simulation. URDME uses Unstructured triangular and tetrahedral meshes to resolve general geometries, and relies on the Reaction-Diffusion Master Equation formalism to model the processes under study. An interface to a mature geometry and mesh handling external software (Comsol Multiphysics) provides for a stable and interactive environment for model construction. The core simulation routines are logically separated from the model building interface and written in a low-level language for computational efficiency. The connection to the geometry handling software is realized via a Matlab interface which facilitates script computing, data management, and post-processing. For practitioners, the software therefore behaves much as an interactive Matlab toolbox. At the same time, it is possible to modify and extend URDME with newly developed simulation routines. Since the overall design effectively hides the complexity of managing the geometry and meshes, this means that newly developed methods may be tested in a realistic setting already at an early stage of development.ConclusionsIn this paper we demonstrate, in a series of examples with high relevance to the molecular systems biology community, that the proposed software framework is a useful tool for both practitioners and developers of spatial stochastic simulation algorithms. Through the combined efforts of algorithm development and improved modeling accuracy, increasingly complex biological models become feasible to study through computational methods. URDME is freely available at http://www.urdme.org.

Facial shaping results from oriented divisions and crowd movements of ectomesenchymal cells during morphogenetic events.

The chemical master equation is solved by a hybrid method coupling a macroscopic, deterministic description with a mesoscopic, stochastic model. The molecular species are divided into one subset where the expected values of the number of molecules are computed and one subset with species with a stochastic variation in the number of molecules. The macroscopic equations resemble the reaction rate equations and the probability distribution for the stochastic variables satisfy a master equation. The probability distribution is obtained by the Stochastic Simulation Algorithm due to Gillespie. The equations are coupled via a summation over the mesoscale variables. This summation is approximated by Monte Carlo and Quasi Monte Carlo methods. The error in the approximations is analyzed. The hybrid method is applied to three chemical systems from molecular cell biology.

The mesoscopic reaction-diffusion master equation (RDME) is a popular modeling framework frequently applied to stochastic reaction-diffusion kinetics in systems biology. The RDME is derived from assumptions about the underlying physical properties of the system, and it may produce unphysical results for models where those assumptions fail. In that case, other more comprehensive models are better suited, such as hard-sphere Brownian dynamics (BD). Although the RDME is a model in its own right, and not inferred from any specific microscale model, it proves useful to attempt to approximate a microscale model by a specific choice of mesoscopic reaction rates. In this paper we derive mesoscopic scale-dependent reaction rates by matching certain statistics of the RDME solution to statistics of the solution of a widely used microscopic BD model: the Smoluchowski model with a Robin boundary condition at the reaction radius of two molecules. We also establish fundamental limits on the range of mesh resolutions for which this approach yields accurate results and show both theoretically and in numerical examples that as we approach the lower fundamental limit, the mesoscopic dynamics approach the microscopic dynamics. We show that for mesh sizes below the fundamental lower limit, results are less accurate. Thus, the lower limit determines the mesh size for which we obtain the most accurate results.