Hybrid stochastic simulations of intracellular reaction-diffusion systems

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.

Simulation of Biochemical Reactions with ANN-Dependent Kinetic Parameter Extraction Method

The measurement of thermodynamic properties of chemical or biological reactions were often confined to experimental means, which produced overall measurements of properties being investigated, but were usually susceptible to pitfalls of being too general. Among the thermodynamic properties that are of interest, reaction rates hold the greatest significance, as they play a critical role in reaction processes where speed is of essence, especially when fast association may enhance binding affinity of reaction molecules. Association reactions with high affinities often involve the formation of a intermediate state, which can be demonstrated by a hyperbolic reaction curve, but whose low abundance in reaction mixture often preclude the possibility of experimental measurement. Therefore, we resorted to computational methods using predefined reaction models that model the intermediate state as the reaction progresses. Here, we present a novel method called AKPE (ANN-Dependent Kinetic Parameter Extraction), our goal is to investigate the association/dissociation rate constants and the concentration dynamics of lowly-populated states (intermediate states) in the reaction landscape. To reach our goal, we simulated the chemical or biological reactions as system of differential equations, employed artificial neural networks (ANN) to model experimentally measured data, and utilized Particle Swarm Optimization (PSO) algorithm to obtain the globally optimum parameters in both the simulation and data fitting. In the Results section, we have successfully modeled a protein association reaction using AKPE, obtained the kinetic rate constants of the reaction, and constructed a full concentration versus reaction time curve of the intermediate state during the reaction. Furthermore, judging from the various validation methods that the method proposed in this paper has strong robustness and accuracy.

Multiscale Stochastic Reaction-Diffusion Algorithms Combining Markov Chain Models with Stochastic Partial Differential Equations

Two multiscale algorithms for stochastic simulations of reaction-diffusion processes are analysed. They are applicable to systems which include regions with significantly different concentrations of molecules. In both methods, a domain of interest is divided into two subsets where continuous-time Markov chain models and stochastic partial differential equations (SPDEs) are used, respectively. In the first algorithm, Markov chain (compartment-based) models are coupled with reaction-diffusion SPDEs by considering a pseudo-compartment (also called an overlap or handshaking region) in the SPDE part of the computational domain right next to the interface. In the second algorithm, no overlap region is used. Further extensions of both schemes are presented, including the case of an adaptively chosen boundary between different modelling approaches.

Hybrid deterministic/stochastic simulation of complex biochemical systems

In a biological cell, cellular functions and the genetic regulatory apparatus are implemented and controlled by complex networks of chemical reactions involving genes, proteins, and enzymes. Accurate computational models are indispensable means for understanding the mechanisms behind the evolution of a complex system, not always explored with wet lab experiments. To serve their purpose, computational models, however, should be able to describe and simulate the complexity of a biological system in many of its aspects. Moreover, it should be implemented by efficient algorithms requiring the shortest possible execution time, to avoid enlarging excessively the time elapsing between data analysis and any subsequent experiment. Besides the features of their topological structure, the complexity of biological networks also refers to their dynamics, that is often non-linear and stiff. The stiffness is due to the presence of molecular species whose abundance fluctuates by many orders of magnitude. A fully stochastic simulation of a stiff system is computationally time-expensive. On the other hand, continuous models are less costly, but they fail to capture the stochastic behaviour of small populations of molecular species. We introduce a new efficient hybrid stochastic-deterministic computational model and the software tool MoBioS (MOlecular Biology Simulator) implementing it. The mathematical model of MoBioS uses continuous differential equations to describe the deterministic reactions and a Gillespie-like algorithm to describe the stochastic ones. Unlike the majority of current hybrid methods, the MoBioS algorithm divides the reactions' set into fast reactions, moderate reactions, and slow reactions and implements a hysteresis switching between the stochastic model and the deterministic model. Fast reactions are approximated as continuous-deterministic processes and modelled by deterministic rate equations. Moderate reactions are those whose reaction waiting time is greater than the fast reaction waiting time but smaller than the slow reaction waiting time. A moderate reaction is approximated as a stochastic (deterministic) process if it was classified as a stochastic (deterministic) process at the time at which it crosses the threshold of low (high) waiting time. A Gillespie First Reaction Method is implemented to select and execute the slow reactions. The performances of MoBios were tested on a typical example of hybrid dynamics: that is the DNA transcription regulation. The simulated dynamic profile of the reagents' abundance and the estimate of the error introduced by the fully deterministic approach were used to evaluate the consistency of the computational model and that of the software tool.

Asynchronous τ-leaping

Stochastic simulation of cell signaling pathways and genetic regulatory networks has contributed to the understanding of cell function; however, investigation of larger, more complicated systems requires computationally efficient algorithms. τ-leaping methods, which improve efficiency when some molecules have high copy numbers, either use a fixed leap size, which does not adapt to changing state, or recalculate leap size at a heavy computational cost. We present a hybrid simulation method for reaction-diffusion systems which combines exact stochastic simulation and τ-leaping in a dynamic way. Putative times of events are stored in a priority queue, which reduces the cost of each step of the simulation. For every reaction and diffusion channel at each step of the simulation the more efficient of an exact stochastic event or a τ-leap is chosen. This new approach removes the inherent trade-off between speed and accuracy in stiff systems which was present in all τ-leaping methods by allowing each reaction channel to proceed at its own pace. Both directions of reversible reactions and diffusion are combined in a single event, allowing bigger leaps to be taken. This improves efficiency for systems near equilibrium where forward and backward events are approximately equally frequent. Comparison with existing algorithms and behaviour for five test cases of varying complexity shows that the new method is almost as accurate as exact stochastic simulation, scales well for large systems, and for various problems can be significantly faster than τ-leaping.

A hybrid continuous-discrete method for stochastic reaction-diffusion processes

Stochastic fluctuations in reaction-diffusion processes often have substantial effect on spatial and temporal dynamics of signal transductions in complex biological systems. One popular approach for simulating these processes is to divide the system into small spatial compartments assuming that molecules react only within the same compartment and jump between adjacent compartments driven by the diffusion. While the approach is convenient in terms of its implementation, its computational cost may become prohibitive when diffusive jumps occur significantly more frequently than reactions, as in the case of rapid diffusion. Here, we present a hybrid continuous-discrete method in which diffusion is simulated using continuous approximation while reactions are based on the Gillespie algorithm. Specifically, the diffusive jumps are approximated as continuous Gaussian random vectors with time-dependent means and covariances, allowing use of a large time step, even for rapid diffusion. By considering the correlation among diffusive jumps, the approximation is accurate for the second moment of the diffusion process. In addition, a criterion is obtained for identifying the region in which such diffusion approximation is required to enable adaptive calculations for better accuracy. Applications to a linear diffusion system and two nonlinear systems of morphogens demonstrate the effectiveness and benefits of the new hybrid method.

Discrete-continuous reaction-diffusion model with mobile point-like sources and sinks

In many applications in soft and biological physics, there are multiple time and length scales involved but often with a distinct separation between them. For instance, in enzyme kinetics, enzymes are relatively large, move slowly and their copy numbers are typically small, while the metabolites (being transformed by these enzymes) are often present in abundance, are small in size and diffuse fast. It seems thus natural to apply different techniques to different time and length levels and couple them. Here we explore this possibility by constructing a stochastic-deterministic discrete-continuous reaction-diffusion model with mobile sources and sinks. Such an approach allows in particular to separate different sources of stochasticity. We demonstrate its application by modelling enzyme-catalysed reactions with freely diffusing enzymes and a heterogeneous source of metabolites. Our calculations suggest that using a higher amount of less active enzymes, as compared to fewer more active enzymes, reduces the metabolite pool size and correspondingly the lag time, giving rise to a faster response to external stimuli. The methodology presented can be extended to more complex systems and offers exciting possibilities for studying problems where spatial heterogeneities, stochasticity or discreteness play a role.

Stochastic chemical kinetics : A review of the modelling and simulation approaches

A review of the physical principles that are the ground of the stochastic formulation of chemical kinetics is presented along with a survey of the algorithms currently used to simulate it. This review covers the main literature of the last decade and focuses on the mathematical models describing the characteristics and the behavior of systems of chemical reactions at the nano- and micro-scale. Advantages and limitations of the models are also discussed in the light of the more and more frequent use of these models and algorithms in modeling and simulating biochemical and even biological processes.

Simulation strategies for calcium microdomains and calcium-regulated calcium channels

In this article, we present an overview of simulation strategies in the context of subcellular domains where calcium-dependent signaling plays an important role. The presentation follows the spatial and temporal scales involved and represented by each algorithm. As an exemplary cell type, we will mainly cite work done on striated muscle cells, i.e. skeletal and cardiac muscle. For these cells, a wealth of ultrastructural, biophysical and electrophysiological data is at hand. Moreover, these cells also express ubiquitous signaling pathways as they are found in many other cell types and thus, the generalization of the methods and results presented here is straightforward.The models considered comprise the basic calcium signaling machinery as found in most excitable cell types including Ca(2+) ions, diffusible and stationary buffer systems, and calcium regulated calcium release channels. Simulation strategies can be differentiated in stochastic and deterministic algorithms. Historically, deterministic approaches based on the macroscopic reaction rate equations were the first models considered. As experimental methods elucidated highly localized Ca(2+) signaling events occurring in femtoliter volumes, stochastic methods were increasingly considered. However, detailed simulations of single molecule trajectories are rarely performed as the computational cost implied is too large. On the mesoscopic level, Gillespie's algorithm is extensively used in the systems biology community and with increasing frequency also in models of microdomain calcium signaling. To increase computational speed, fast approximations were derived from Gillespie's exact algorithm, most notably the chemical Langevin equation and the τ-leap algorithm. Finally, in order to integrate deterministic and stochastic effects in multiscale simulations, hybrid algorithms are increasingly used. These include stochastic models of ion channels combined with deterministic descriptions of the calcium buffering and diffusion system on the one hand, and algorithms that switch between deterministic and stochastic simulation steps in a context-dependent manner on the other. The basic assumptions of the listed methods as well as implementation schemes are given in the text. We conclude with a perspective on possible future developments of the field.

Multi-scale stochastic simulation of diffusion-coupled agents and its application to cell culture simulation

Many biological systems consist of multiple cells that interact by secretion and binding of diffusing molecules, thus coordinating responses across cells. Techniques for simulating systems coupling extracellular and intracellular processes are very limited. Here we present an efficient method to stochastically simulate diffusion processes, which at the same time allows synchronization between internal and external cellular conditions through a modification of Gillespie's chemical reaction algorithm. Individual cells are simulated as independent agents, and each cell accurately reacts to changes in its local environment affected by diffusing molecules. Such a simulation provides time-scale separation between the intra-cellular and extra-cellular processes. We use our methodology to study how human monocyte-derived dendritic cells alert neighboring cells about viral infection using diffusing interferon molecules. A subpopulation of the infected cells reacts early to the infection and secretes interferon into the extra-cellular medium, which helps activate other cells. Findings predicted by our simulation and confirmed by experimental results suggest that the early activation is largely independent of the fraction of infected cells and is thus both sensitive and robust. The concordance with the experimental results supports the value of our method for overcoming the challenges of accurately simulating multiscale biological signaling systems.

Exact and approximate stochastic simulation of intracellular calcium dynamics

In simulations of chemical systems, the main task is to find an exact or approximate solution of the chemical master equation (CME) that satisfies certain constraints with respect to computation time and accuracy. While Brownian motion simulations of single molecules are often too time consuming to represent the mesoscopic level, the classical Gillespie algorithm is a stochastically exact algorithm that provides satisfying results in the representation of calcium microdomains. Gillespie's algorithm can be approximated via the tau-leap method and the chemical Langevin equation (CLE). Both methods lead to a substantial acceleration in computation time and a relatively small decrease in accuracy. Elimination of the noise terms leads to the classical, deterministic reaction rate equations (RRE). For complex multiscale systems, hybrid simulations are increasingly proposed to combine the advantages of stochastic and deterministic algorithms. An often used exemplary cell type in this context are striated muscle cells (e.g., cardiac and skeletal muscle cells). The properties of these cells are well described and they express many common calcium-dependent signaling pathways. The purpose of the present paper is to provide an overview of the aforementioned simulation approaches and their mutual relationships in the spectrum ranging from stochastic to deterministic algorithms.