Iain Hepburn

STEPS: efficient simulation of stochastic reaction–diffusion models in realistic morphologies

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/

Stochastic Calcium Mechanisms Cause Dendritic Calcium Spike Variability

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.

The role of dendritic spine morphology in the compartmentalization and delivery of surface receptors

Since AMPA receptors are major molecular players in both short- and long-term plasticity, it is important to identify the time-scales of and factors affecting the lateral diffusion of AMPARs on the dendrite surface. Using a mathematical model, we study how the dendritic spine morphology affects two processes: (1) compartmentalization of the surface receptors in a single spine to retain local chemistry and (2) the delivery of receptors to the post-synaptic density (PSD) of spines via lateral diffusion following insertion onto the dendrite shaft. Computing the mean first passage time (MFPT) of surface receptors on a sample of real spine morphologies revealed that a constricted neck and bulbous head serve to compartmentalize receptors, consistent with previous works. The residence time of a Brownian diffusing receptor on the membrane of a single spine was computed to be ∼ 5 s. We found that the location of the PSD corresponds to the location at which the maximum MFPT occurs, the position that maximizes the residence time of a diffusing receptor. Meanwhile, the same geometric features of the spine that compartmentalize receptors inhibit the recruitment of AMPARs via lateral diffusion from dendrite insertion sites. Spines with narrow necks will trap a smaller fraction of diffusing receptors in the their PSD when considering competition for receptors between the spines, suggesting that ideal geometrical features involve a tradeoff depending on the intent of compartmentalizing the current receptor pool or recruiting new AMPARs in the PSD. The ultimate distribution of receptors among the spine PSDs by lateral diffusion from the dendrite shaft is an interplay between the insertion location and the shape and locations of both the spines and their PSDs. The time-scale for delivery of receptors to the PSD of spines via lateral diffusion was computed to be ∼ 60 s.

Efficient calculation of the quasi-static electrical potential on a tetrahedral mesh and its implementation in STEPS

We describe a novel method for calculating the quasi-static electrical potential on tetrahedral meshes, which we call E-Field. The E-Field method is implemented in STEPS, which performs stochastic spatial reaction-diffusion computations in tetrahedral-based cellular geometry reconstructions. This provides a level of integration between electrical excitability and spatial molecular dynamics in realistic cellular morphology not previously achievable. Deterministic solutions are also possible. By performing the Rallpack tests we demonstrate the accuracy of the E-Field method. Efficient node ordering is an important practical consideration, and we find that a breadth-first search provides the best solutions, although principal axis ordering suffices for some geometries. We discuss potential applications and possible future directions, and predict that the E-Field implementation in STEPS will play an important role in the future of multiscale neural simulations.

Accurate approximation and MPI parallelization of spatial stochastic reaction-diffusion in STEPS.

Spatial stochastic reaction-diffusion simulations have become an important component of molecular modeling in Computational Neuroscience, as shown in a growing number of recent studies including our previous work in which we show that stochastic effects, in particular stochastic calcium dynamics, contribute to Purkinje cell calcium burst variability [1]. However, the computational cost of exact, serial algorithms such as Gillespie’s direct method as implemented in STEPS [2] in which every reaction event and cross-subvolume diffusion event is simulated, impose severe restrictions on model complexity. Serial solutions limit the number of molecular species and their reaction channels that can be represented, as well as imposing restrictive upper-limits on the number of subvolumes in the discretized space. It is desirable to provide a solution that allows investigation of stochastic effects in larger systems, with complex meshes of many hundreds of thousands or millions of subvolumes representing complex neuronal morphology, in projects such as the Human Brain project [3], which is not possible with serial solutions. Therefore, parallel solutions for stochastic reaction-diffusion is vital for the future of this field. There are common challenges of parallelization for reaction-diffusion simulators as well as some unique challenges for STEPS, in which space is discretized into irregular tetrahedrons and not cubes that are commonly used by other reaction-diffusion simulators. Exact solutions are not beneficial, where stochastic diffusion across geometry partitions result in frequent conflict between nodes, requiring regular costly rollbacks. Therefore it is essential that an approximation algorithm is applied. Such an approximation must minimise loss of accuracy whilst maximising performance gain. Existing solutions such as the Gillespie-Multi-Particle method (GMP) [4,5] have been designed for cubic subvolumes, which have the advantage of regularity (yet restrict morphological accuracy [2]). We review such approaches from a theoretical standpoint and find that a systematic diffusion slowing error exists with naive application to irregular subvolumes, which have local scaled diffusion coefficients and therefore localized expected diffusion times. We present a method that is tailored for irregular subvolumes, with a defined theoretical upper limit to the communication time-step. The method could, however, be generalised for different subvolume configurations where adaptable step-sizes could be beneficial. We further expand our method to allow multiple particle leaps, which for some systems this can greatly increase the minimum acceptable step-size and thus benefit performance whilst still ensuring acceptable accuracy. We demonstrate accuracy of our approximation by comparing to analytical solutions on a number of test cases. Our method is implemented on an MPI framework, and we present performance gains under two scenariosglobal application of the method, and application only at partition boundaries. We discuss the potential of this method for future large-scale parallel simulations in STEPS.

A Model of Induction of Cerebellar Long-Term Depression Including RKIP Inactivation of Raf and MEK

We report an updated stochastic model of cerebellar Long Term Depression (LTD) with improved realism. Firstly, we verify experimentally that dissociation of Raf kinase inhibitor protein (RKIP) from Mitogen-activated protein kinase kinase (MEK) is required for cerebellar LTD and add this interaction to an earlier published model, along with the known requirement of dissociation of RKIP from Raf kinase. We update Ca2+ dynamics as a constant-rate influx, which captures experimental input profiles accurately. We improve α-amino-3-hydroxy-5-methyl-4 isoxazolepropionic acid (AMPA) receptor interactions by adding phosphorylation and dephosphorylation of AMPA receptors when bound to glutamate receptor interacting protein (GRIP). The updated model is tuned to reproduce experimental Ca2+ peak vs. LTD amplitude curves at four different Ca2+ pulse durations as closely as possible. We find that the updated model is generally more robust with these changes, yet we still observe some sensitivity of LTD induction to copy number of the key signaling molecule Protein kinase C (PKC). We predict natural variability in this number by stochastic diffusion may influence the probabilistic LTD response to Ca2+ input in Purkinje cell spines and propose this as an extra source of stochasticity that may be important also in other signaling systems.

Accurate reaction-diffusion operator splitting on tetrahedral meshes for parallel stochastic molecular simulations

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, which has led to the development of parallel methods that can take advantage of the power of modern supercomputers in recent years. We systematically test suggested components of stochastic reaction-diffusion operator splitting in the literature and discuss their effects on accuracy. We introduce an operator splitting implementation for irregular meshes that enhances accuracy with minimal performance cost. We test a range of models in small-scale MPI simulations from simple diffusion models to realistic biological models and find that multi-dimensional geometry partitioning is an important consideration for optimum performance. We demonstrate performance gains of 1-3 orders of magnitude in the parallel implementation, with peak performance strongly dependent on model specification.

Improving performance of the STochastic Engine for Pathway Simulation (STEPS)

STEPS (http://steps.sourceforge.net) is a GNU-licensed simulation platform that uses an extension of Gillespies SSA [1] to deal with reactions and diffusion of molecules in 3D reconstructions of neuronal morphology and tissue [2]. In STEPS, the diffusion of molecules is simulated as diffusive fluxes between tetrahedral elements in the mesh, represented by a series of first-order reactions. STEPS has been used in various research projects where its overall performance was considered adequate. However, as more complex models are being developed and investigated, total model simulation times became an issue leading to requests for a faster reaction-diffusion simulator. Here we discuss a number of strategies we have followed to improve STEPS performance at different levels. The previous implementation of STEPS adapted Gibson and Bruck’s enhancement of the Direct method [3] with a k-ary tree structure. The overall complexity of searching and updating of this method is O(logkN), giving N as the total number of kinetic processes in the system and k as the branch width of the tree. Although this is a great improvement comparing to the O(N) complexity of Gillespie’s original method, the logarithmic dependency of N becomes a critical limitation of performance, particularly in simulations with high amount of tetrahedral elements and consequently large N. Recently a constant-time version of the Direct method with a slightly more complex data structure has been introduced [4], providing an attractive searching and updating complexity of O(1), independent of N. This method is implemented in the new version of STEPS. In the poster, I will present the validation of the new simulator by comparing the stochastic simulation results with analytical solutions. I will also present the improvement made by adapting the new method. Several approaches to parallelizing STEPS have been investigated. In a realistic simulation, the state of the system, including molecule distribution, reaction/diffusion rates, may be frequently read and/or modified during the simulation. Some of these operations could be greatly parallelized in shared memory architecture, efficiently reducing the cost of data accessing. Previous studies [5] also suggested that the Direct SSA could be parallelized as a Parallel Discrete-Event Simulation (PDES). However, the performance of such a parallelization significantly depends on the model itself as well as the approach of simulation decomposition. Some examples of above methods will be present and discussed in the poster.

Tetrahedral mesh generation and visualization for stochastic reaction-diffusion simulation

The generation of high quality three-dimensional meshes is important in voxel-based computational modeling and simulation of stochastic reaction-diffusion systems. Compared with cubical meshes used in some reaction-diffusion simulators such as MesoRD [1], tetrahedral meshes are more flexible in representing complex geometry boundaries, for instance, the spherical head of a dendritic spine. We describe the interaction between tetrahedral mesh generators (CUBIT [2], TetGen [3]) and our STochastic Engine for Pathway Simulation (STEPS) [4], including the generation and importing of meshes, as well as the visualization of the simulations in a 3D graphical environment. STEPS is a GNU-licensed simulation platform that uses an extension of Gillespies SSA [5] to deal with reactions and diffusion of molecules in 3D reconstructions of neuronal morphology and tissue [4]. In order to support a range of research goals, it provides a flexible Python frontend together with an efficient C++ backend. In STEPS, the diffusion of molecules is simulated as diffusive fluxes between tetrahedral elements in the mesh, represented by a series of first-order reactions. Since many mesh generating methods and software packages are available in public, instead of implementing our own mesh generator, STEPS provides a generic, extendable Python-based import interface for different tetrahedral mesh formats. Thus the first step is to create a suitable mesh, which can be challenging. In our research [6] it involves the reconstruction of anatomical/morphological structures from combinations of simple geometry primitives like cylinders and spheres. Once the mesh file is generated, it can be imported to the simulations using either pre-implemented import functions (currently supporting Abaqus and TetGen mesh formats) or user-developed functions by adapting the STEPS element proxy interface. The element proxy interface performs automatic mapping between mesh elements in external files and the ones used by the STEPS simulation. Monitoring mesh quality is possible via mesh generators or the STEPS internal quality metric. This metric calculates the Radius-Edge Ratio (RER) for each tetrahedron, given by dividing the radius of the tetrahedrons circumsphere with the length of the shortest edge. The smaller this value, the more regular the tetrahedron. The minimum value of this metric is given by computing the RER for a fully regular tetrahedron, that is, sqrt(6)/4. Although the visualization of stochastic reaction-diffusion simulations with 3D meshes is optional due to the high number of reactions/diffusions that take place during a simulation, it is sometimes required for demonstration and education. Animations of such a simulation can be made by a script-controlled 3D graphical environment, for example the python frontend of CUBIT, with user-specified data extracted from the STEPS simulation. Alternatively, STEPS provides a preliminary Python-based visual toolkit for general purposes of displaying meshes and simple reaction-diffusion simulations. Overall, the collaborative workflow between STEPS and external mesh generators discussed here has shown its flexibility and will facilitate a broad range of applications in computational pathway modeling and simulation.

Implementation of parallel spatial stochastic reaction-diffusion simulation in STEPS

Spatial stochastic reaction-diffusion simulation has been recognized as an essential modeling tool in computational neuroscience studies of signaling pathways, as shown in an increasing number of recent studies such as our previous work on the stochastic effects of calcium dynamics in Purkinje cells [1]. A critical performance issue arises when attempting to model large pathway models with complex morphologies, for example the one proposed in Human Brain Project [2], due to the serial nature of Gillespies direct method [3], the fundamental algorithm of many spatial stochastic reaction-diffusion simulators including STEPS [4]. Various solutions have been proposed to improve the computational efficiency of the Gillespie method, including the tau-leaping approximation [5], most of which, however, remain serial implementations. The need of parallel implementation of stochastic reaction-diffusion simulators has become urgent, as the scale and complexity of model being studied surpasses the speedup gained from hardware upgrade and algorithm improvement. However, such a task is not trivial as the original Gillespie SSA is known to be extremely serial. In CNS2014 we proposed a parallel solution to approximate diffusion events in STEPS [6], which significantly improves the performance while maintaining high accuracy. We now further improve this solution by introducing a multinomial algorithm for fast diffusion direction selection of multiple molecules in a single subvolume. We also combine this diffusion approximation with a new operator splitting solution for reaction events in the SSA system. The combined solution is implemented in STEPS as its first parallel solver named TetOpSplit. Current implementation of TetOpSplit uses MPI as its parallel protocol and aims to provide solutions for large scale simulations such as whole cell reaction-diffusion models, in modern supercomputers like Blue Gene. In this poster we discuss the difficulties we encountered during the transformation from the serial TetOpSplit algorithm to its parallel counterpart, as well as our solutions. We also provide performance results of the new solver via different examples, and compare them with the results gained from our original serial SSA implementation.