Smoldyn on graphics processing units: massively parallel Brownian dynamics simulations

Space is a very important aspect in the simulation of biochemical systems; recently, the need for simulation algorithms able to cope with space is becoming more and more compelling. Complex and detailed models of biochemical systems need to deal with the movement of single molecules and particles, taking into consideration localized fluctuations, transportation phenomena, and diffusion. A common drawback of spatial models lies in their complexity: models can become very large, and their simulation could be time consuming, especially if we want to capture the systems behavior in a reliable way using stochastic methods in conjunction with a high spatial resolution. In order to deliver the promise done by systems biology to be able to understand a system as whole, we need to scale up the size of models we are able to simulate, moving from sequential to parallel simulation algorithms. In this paper, we analyze Smoldyn, a widely diffused algorithm for stochastic simulation of chemical reactions with spatial resolution and single molecule detail, and we propose an alternative, innovative implementation that exploits the parallelism of Graphics Processing Units (GPUs). The implementation executes the most computational demanding steps (computation of diffusion, unimolecular, and bimolecular reaction, as well as the most common cases of molecule-surface interaction) on the GPU, computing them in parallel on each molecule of the system. The implementation offers good speed-ups and real time, high quality graphics output     

Bridging from molecular simulation to biochemical networks

How can we make the connection between the three-dimensional structures of individual proteins and understanding how complex biological systems involving many proteins work? The modelling and simulation of protein structures can help to answer this question for systems ranging from multimacromolecular complexes to organelles and cells. On one hand, multiscale modelling and simulation techniques are advancing to permit the spatial and temporal properties of large systems to be simulated using atomic-detail structures. On the other hand, the estimation of kinetic parameters for the mathematical modelling of biochemical pathways using protein structure information provides a basis for iterative manipulation of biochemical pathways guided by protein structure. Recent advances include the structural modelling of protein complexes on the genomic level, novel coarse-graining strategies to increase the size of the system and the time span that can be simulated, and comparative molecular field analyses to estimate enzyme kinetic parameters.     

Accelerated maximum likelihood parameter estimation for stochastic biochemical systems

ABSTRACT: BACKGROUND: A prerequisite for the mechanistic simulation of a biochemical system is detailed knowledge of its kinetic parameters. Despite recent experimental advances, the estimation of unknown parameter values from observed data is still a bottleneck for obtaining accurate simulation results. Many methods exist for parameter estimation in deterministic biochemical systems; methods for discrete stochastic systems are less well developed. Given the probabilistic nature of stochastic biochemical models, a natural approach is to choose parameter values that maximize the probability of the observed data with respect to the unknown parameters, a.k.a. the maximum likelihood parameter estimates (MLEs). MLE computation for all but the simplest models requires the simulation of many system trajectories that are consistent with experimental data. For models with unknown parameters, this presents a computational challenge, as the generation of consistent trajectories can be an extremely rare occurrence. RESULTS: We have developed Monte Carlo Expectation-Maximization with Modified Cross-Entropy Method (MCEM2): an accelerated method for calculating MLEs that combines advances in rare event simulation with a computationally efficient version of the Monte Carlo expectation-maximization (MCEM) algorithm. Our method requires no prior knowledge regarding parameter values, and it automatically provides a multivariate parameter uncertainty estimate. We applied the method to five stochastic systems of increasing complexity, progressing from an analytically tractable pure-birth model to a computationally demanding model of yeast-polarization. Our results demonstrate that MCEM2 substantially accelerates MLE computation on all tested models when compared to a stand-alone version of MCEM. Additionally, we show how our method identifies parameter values for certain classes of models more accurately than two recently proposed computationally efficient methods. CONCLUSIONS: This work provides a novel, accelerated version of a likelihood-based parameter estimation method that can be readily applied to stochastic biochemical systems. In addition, our results suggest opportunities for added efficiency improvements that will further enhance our ability to mechanistically simulate biological processes.     

Simulation tools for biochemical networks: evaluation of performance and usability

MOTIVATION: Simulation of dynamic biochemical systems is receiving considerable attention due to increasing availability of experimental data of complex cellular functions. Numerous simulation tools have been developed for numerical simulation of the behavior of a system described in mathematical form. However, there exist only a few evaluation studies of these tools. Knowledge of the properties and capabilities of the simulation tools would help bioscientists in building models based on experimental data. RESULTS: We examine selected simulation tools that are intended for the simulation of biochemical systems. We choose four of them for more detailed study and perform time series simulations using a specific pathway describing the concentration of the active form of protein kinase C. We conclude that the simulation results are convergent between the chosen simulation tools. However, the tools differ in their usability, support for data transfer to other programs and support for automatic parameter estimation. From the experimentalists' point of view, all these are properties that need to be emphasized in the future.     

Quantitative modeling of biochemical networks

Today different database systems for molecular structures (genes and proteins) and metabolic pathways are available. All these systems are characterized by the static data representation. For progress in biotechnology the dynamic representation of this data is important. The metabolism can be characterized as a complex biochemical network. Different models for the quantitative simulation of biochemical networks are discussed, but no useful formalization is available. This paper shows that the theory of Petrinets is useful for the quantitative modeling of biochemical networks.     

Generalized Chemokinetic Method for Gene Network Simulation

Development of methods for mathematical simulation of biological systems and building specific simulations is an important trend of bioinformatics. Here we describe the method of generalized chemokinetic simulation generating flexible and adequate simulations of various biological systems. Adequate simulations of complex nonlinear gene networks—control system of cholesterol by synthesis in the cell and erythrocyte differentiation and maturation—are given as examples. The simulations were expressed in terms of unit processes—biochemical reactions. Optimal sets of parameters were determined and the systems were numerically simulated under various conditions. The simulations allow us to study the possible functional conditions of these gene networks, calculate the consequences of mutations, and define optimal strategies for their correction including therapeutic ones. A graphical user interface for these simulations is available at http://wwwmgs.bionet.nsc.ru/systems/MGL/GeneNet/     

CADLIVE toolbox for MATLAB: automatic dynamic modeling of biochemical networks with comprehensive system analysis

Mathematical modeling has become a standard technique to understand the dynamics of complex biochemical systems. To promote the modeling, we had developed the CADLIVE dynamic simulator that automatically converted a biochemical map into its associated mathematical model, simulated its dynamic behaviors and analyzed its robustness. To enhance the feasibility by CADLIVE and extend its functions, we propose the CADLIVE toolbox available for MATLAB, which implements not only the existing functions of the CADLIVE dynamic simulator, but also the latest tools including global parameter search methods with robustness analysis. The seamless, bottom-up processes consisting of biochemical network construction, automatic construction of its dynamic model, simulation, optimization, and S-system analysis greatly facilitate dynamic modeling, contributing to the research of systems biology and synthetic biology. This application can be freely downloaded from http://www.cadlive.jp/CADLIVE_MATLAB/ together with an instruction.     

Weighted next reaction method and parameter selection for efficient simulation of rare events in biochemical reaction systems

The weighted stochastic simulation algorithm (wSSA) recently developed by Kuwahara and Mura and the refined wSSA proposed by Gillespie et al. based on the importance sampling technique open the door for efficient estimation of the probability of rare events in biochemical reaction systems. In this paper, we first apply the importance sampling technique to the next reaction method (NRM) of the stochastic simulation algorithm and develop a weighted NRM (wNRM). We then develop a systematic method for selecting the values of importance sampling parameters, which can be applied to both the wSSA and the wNRM. Numerical results demonstrate that our parameter selection method can substantially improve the performance of the wSSA and the wNRM in terms of simulation efficiency and accuracy.     

Transition from stochastic to deterministic behavior in calcium oscillations

Simulation and modeling is becoming more and more important when studying complex biochemical systems. Most often, ordinary differential equations are employed for this purpose. However, these are only applicable when the numbers of participating molecules in the biochemical systems are large enough to be treated as concentrations. For smaller systems, stochastic simulations on discrete particle basis are more accurate. Unfortunately, there are no general rules for determining which method should be employed for exactly which problem to get the most realistic result. Therefore, we study the transition from stochastic to deterministic behavior in a widely studied system, namely the signal transduction via calcium, especially calcium oscillations. We observe that the transition occurs within a range of particle numbers, which roughly corresponds to the number of receptors and channels in the cell, and depends heavily on the attractive properties of the phase space of the respective systems dynamics. We conclude that the attractive properties of a system, expressed, e.g., by the divergence of the system, are a good measure for determining which simulation algorithm is appropriate in terms of speed and realism.     

钙离子振荡对肝糖磷酸化酶激活的随机效应研究

在复杂生化系统的研究过程中,仿真与建模变得越来越重要.对于参与分子数量比较大的生化系统,通常可以采用常微分方程来解决这一问题.对于分子数量比较小的系统,离散粒子基础上的随机模拟方法更精确.然而目前还没有明确的理论方法来确定,对于实际问题用哪种方法能得到更合理的结果.因此需要在一个普遍研究的体系中,通过Ca~(2+)振荡传导信号来研究从随机行为到确定行为的过渡过程.本文以肝细胞中Ca~(2+)振荡对肝糖磷酸化酶激活随机效应为例,讨论了利用随机微分方程来解决分子数量比较小的生化系统的仿真与建模问题,利用细胞内Ca~(2+)有关的Li-Rinzel随机模型,研究了在磷酸化酶降解肝糖的磷酸化-去磷酸化作用循环过程中,三磷酸肌醇受体通道(IP_3R)释放Ca~(2+)的调控作用.结果表明,肝糖磷酸化酶的激活率随受体通道IP_3R的总数增大而减弱,而且三磷酸肌醇浓度比较小时出现相干共振.     

Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology

We present a Systems Biology Toolbox for the widely used general purpose mathematical software MATLAB. The toolbox offers systems biologists an open and extensible environment, in which to explore ideas, prototype and share new algorithms, and build applications for the analysis and simulation of biological and biochemical systems. Additionally it is well suited for educational purposes. The toolbox supports the Systems Biology Markup Language (SBML) by providing an interface for import and export of SBML models. In this way the toolbox connects nicely to other SBML-enabled modelling packages. Models are represented in an internal model format and can be described either by entering ordinary differential equations or, more intuitively, by entering biochemical reaction equations. The toolbox contains a large number of analysis methods, such as deterministic and stochastic simulation, parameter estimation, network identification, parameter sensitivity analysis and bifurcation analysis.     

Stochastic P systems and the simulation of biochemical processes with dynamic compartments

We introduce a sequential rewriting strategy for P systems based on Gillespie's stochastic simulation algorithm, and show that the resulting formalism of stochastic P systems makes it possible to simulate biochemical processes in dynamically changing, nested compartments. Stochastic P systems have been implemented using the spatially explicit programming language MGS. Implementation examples include models of the Lotka-Volterra auto-catalytic system, and the life cycle of the Semliki Forest virus.     

Extended Kalman Filter for Estimation of Parameters in Nonlinear State-Space Models of Biochemical Networks

It is system dynamics that determines the function of cells, tissues and organisms. To develop mathematical models and estimate their parameters are an essential issue for studying dynamic behaviors of biological systems which include metabolic networks, genetic regulatory networks and signal transduction pathways, under perturbation of external stimuli. In general, biological dynamic systems are partially observed. Therefore, a natural way to model dynamic biological systems is to employ nonlinear state-space equations. Although statistical methods for parameter estimation of linear models in biological dynamic systems have been developed intensively in the recent years, the estimation of both states and parameters of nonlinear dynamic systems remains a challenging task. In this report, we apply extended Kalman Filter (EKF) to the estimation of both states and parameters of nonlinear state-space models. To evaluate the performance of the EKF for parameter estimation, we apply the EKF to a simulation dataset and two real datasets: JAK-STAT signal transduction pathway and Ras/Raf/MEK/ERK signaling transduction pathways datasets. The preliminary results show that EKF can accurately estimate the parameters and predict states in nonlinear state-space equations for modeling dynamic biochemical networks.     

A gradual update method for simulating the steady-state solution of stiff differential equations in metabolic circuits

Numerical simulation of differential equation systems plays a major role in the understanding of how metabolic network models generate particular cellular functions. On the other hand, the classical and technical problems for stiff differential equations still remain to be solved, while many elegant algorithms have been presented. To relax the stiffness problem, we propose new practical methods: the gradual update of differential-algebraic equations based on gradual application of the steady-state approximation to stiff differential equations, and the gradual update of the initial values in differential-algebraic equations. These empirical methods show a high efficiency for simulating the steady-state solutions for the stiff differential equations that existing solvers alone cannot solve. They are effective in extending the applicability of dynamic simulation to biochemical network models. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Biochemical switching device realizing McCulloch-Pitts type equation

We analyzed cyclic enzyme systems, one of the best candidates for biochemical switching devices, especially focusing on their control mode against external perturbations. Since these systems have the reliability of ON-OFF types of operation (McCulloch-Pitts' neuronic equation), we shall present here the mechanical difference between these systems and electronic switching circuit, especially on the mnemonic mechanism of biochemical switching devices.     

The dynamic systems approach to control and regulation of intracellular networks

Systems theory and cell biology have enjoyed a long relationship that has received renewed interest in recent years in the context of systems biology. The term 'systems' in systems biology comes from systems theory or dynamic systems theory: systems biology is defined through the application of systems- and signal-oriented approaches for an understanding of inter- and intra-cellular dynamic processes. The aim of the present text is to review the systems and control perspective of dynamic systems. The biologist's conceptual framework for representing the variables of a biochemical reaction network, and for describing their relationships, are pathway maps. A principal goal of systems biology is to turn these static maps into dynamic models, which can provide insight into the temporal evolution of biochemical reaction networks. Towards this end, we review the case for differential equation models as a 'natural' representation of causal entailment in pathways. Block-diagrams, commonly used in the engineering sciences, are introduced and compared to pathway maps. The stimulus-response representation of a molecular system is a necessary condition for an understanding of dynamic interactions among the components that make up a pathway. Using simple examples, we show how biochemical reactions are modelled in the dynamic systems framework and visualized using block-diagrams.     

Modeling adaptive biological systems

During the evolution of many systems found in nature, both the system composition and the interactions between components will vary. Equating the dimension with the number of different components, a system which adds or deletes components belongs to a class of dynamical systems with a finite dimensional phase space of variable dimension. We present two models of biochemical systems with a variable phase space, a model of autocatalytic reaction networks in the prebiotic soup and a model of the idiotypic network of the immune system. Each model contains characteristic meta-dynamical rules for constructing equations of motion from component properties. The simulation of each model occurs on two levels. On one level, the equations of motion are integrated to determine the state of each component. On a second level, algorithms which approximate physical processes in the real system are employed to change the equations of motion. Models with meta-dynamical rules possess several advantages for the study of evolving systems. First, there are no explicit fitness functions to determine how the components of the model rank in terms of survivability. The success of any component is a function of its relationship to the rest of the system. A second advantage is that since the phase space representation of the system is always finite but continually changing, we can explore a potentially infinite phase space which would otherwise be inaccessible with finite computer resources. Third, the enlarged capacity of systems with meta-dynamics for variation allows us to conduct true evolution experiments. The modeling methods presented here can be applied to many real biological systems. In the two studies we present, we are investigating two apparent properties of adaptive networks. With the simulation of the prebiotic soup, we are most interested in how a chemical reaction network might emerge from an initial state of relative disorder. With the study of the immune system, we study the self-regulation of the network including its ability to distinguish between species which are part of the network and those which are not.