BackgroundModels of cellular molecular systems are built from components such as biochemical reactions (including interactions between ligands and membrane-bound proteins), conformational changes and active and passive transport. A discrete, stochastic description of the kinetics is often essential to capture the behavior of the system accurately. Where spatial effects play a prominent role the complex morphology of cells may have to be represented, along with aspects such as chemical localization and diffusion. This high level of detail makes efficiency a particularly important consideration for software that is designed to simulate such systems.ResultsWe describe STEPS, a stochastic reaction–diffusion simulator developed with an emphasis on simulating biochemical signaling pathways accurately and efficiently. STEPS supports all the above-mentioned features, and well-validated support for SBML allows many existing biochemical models to be imported reliably. Complex boundaries can be represented accurately in externally generated 3D tetrahedral meshes imported by STEPS. The powerful Python interface facilitates model construction and simulation control. STEPS implements the composition and rejection method, a variation of the Gillespie SSA, supporting diffusion between tetrahedral elements within an efficient search and update engine. Additional support for well-mixed conditions and for deterministic model solution is implemented. Solver accuracy is confirmed with an original and extensive validation set consisting of isolated reaction, diffusion and reaction–diffusion systems. Accuracy imposes upper and lower limits on tetrahedron sizes, which are described in detail. By comparing to Smoldyn, we show how the voxel-based approach in STEPS is often faster than particle-based methods, with increasing advantage in larger systems, and by comparing to MesoRD we show the efficiency of the STEPS implementation.ConclusionSTEPS simulates models of cellular reaction–diffusion systems with complex boundaries with high accuracy and high performance in C/C++, controlled by a powerful and user-friendly Python interface. STEPS is free for use and is available at http://steps.sourceforge.net/

Bursts of dendritic calcium spikes play an important role in excitability and synaptic plasticity in many types of neurons. In single Purkinje cells, spontaneous and synaptically evoked dendritic calcium bursts come in a variety of shapes with a variable number of spikes. The mechanisms causing this variability have never been investigated thoroughly. In this study, a detailed computational model using novel simulation routines is applied to identify the roles that stochastic ion channels, spatial arrangements of ion channels, and stochastic intracellular calcium have toward producing calcium burst variability. Consistent with experimental recordings from rats, strong variability in the burst shape is observed in simulations. This variability persists in large model sizes in contrast to models containing only voltage-gated channels, where variability reduces quickly with increase of system size. Phase plane analysis of Hodgkin-Huxley spikes and of calcium bursts identifies fluctuation in phase space around probabilistic phase boundaries as the mechanism determining the dependence of variability on model size. Stochastic calcium dynamics are the main cause of calcium burst fluctuations, specifically the calcium activation of mslo/BK-type and SK2 channels. Local variability of calcium concentration has a significant effect at larger model sizes. Simulations of both spontaneous and synaptically evoked calcium bursts in a reconstructed dendrite show, in addition, strong spatial and temporal variability of voltage and calcium, depending on morphological properties of the dendrite. Our findings suggest that stochastic intracellular calcium mechanisms play a crucial role in dendritic calcium spike generation and are therefore an essential consideration in studies of neuronal excitability and plasticity.

Spatial stochastic molecular simulations in biology are limited by the intense computation required to track molecules in space either in a discrete time or discrete space framework, meaning that the serial limit has already been reached in sub-cellular models. This calls for parallel simulations that can take advantage of the power of modern supercomputers; however exact methods are known to be inherently serial. We introduce an operator splitting implementation for irregular grids with a novel method to improve accuracy, and demonstrate potential for scalable parallel simulations in an initial MPI version. We foresee that this groundwork will enable larger scale, whole-cell stochastic simulations in the near future.