Journal of Biomedicine and Biotechnology

Volume 2011, Article ID 572492, 5 pages

http://dx.doi.org/10.1155/2011/572492

## Exact and Approximate Stochastic Simulation of Intracellular Calcium Dynamics

Medical Biophysics Group, Institute of Physiology and Pathophysiology, University of Heidelberg, 69120 Heidelberg, Germany

Received 1 July 2011; Accepted 25 August 2011

Academic Editor: Aikaterini Kontrogianni-Konstantopoulos

Copyright © 2011 Nicolas Wieder et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. Introduction

Ca^{2+} is a vital second messenger in many cell types. The ubiquitous appearance as well as its often crucial role in functional important signaling pathways are the reasons for the great scientific interest in Ca^{2+} dynamics. Oscillations in the typically low-intracellular Ca^{2+} concentration carry information which is processed (filtered) via signal cascades (e.g., calmodulin-dependent pathways) to induce a variety of cellular responses on different spatiotemporal scales (e.g., muscle contraction [1], fertilization [2], and regulation of gene expression [3]).

Oscillations are induced by sequences of random calcium “puffs” or “sparks”, that are highly localized Ca^{2+} release events from intracellular Ca^{2+} stores (sarcoplasmic and endoplasmatic reticulum) [4]. However, the link between the dynamics of individual molecules involved in microdomain Ca^{2+} dynamics (e.g., Ca^{2+} channels) and the resulting cellular responses are not completely understood yet. This is where computational simulation algorithms come into play. The great number of biological functions as well as the simple structure of Ca^{2+} ions are the reason for its well-described chemical and physical properties. Those circumstances make the in silico examination of Ca^{2+} dynamics very promising. Microdomains controlled by Ca^{2+} channels play an important role in the context of Ca^{2+} oscillations and waves [5]. For an accurate model of these domains, we have to consider stochastic processes at the mesoscopic level. Many regulatory and signaling processes are only inadequately described by deterministic simulation approaches [6].

An often used exemplary cell type in this context is given by striated muscle cells (e.g., cardiac and skeletal muscle cells) where a large amount of electrophysiological, biochemical, and ultrastructural data exists. Furthermore, these cells express many of the ubiquitous calcium-dependent signaling pathways which make it easy to transfer the methods and results to other cell types.

#### 2. Background

Every chemical system is defined by a set of properties. First of all, there are chemical species, , in the system volume . The state vector denotes the number of molecules of each species at time in . Events define transitions between different system states, where an event is any biophysical process which is characterized by a rate constant. It is possible to derive an event propensity by its stochastic rate constant and the stoichiometric coefficients of the underlying process. A chemical system contains events which are represented through a state change vector , where component denotes the change in the number of molecules due to the event . Even more complex processes, which do not affect , but other system properties, can easily be considered, for instance, voltage-gated ion channels have been introduced into stochastic simulation algorithms. To introduce the continuously varying membrane potential regulating the channel, its effect on the subunit transition rates was modeled as a time-varying stochastic event rate [10]. In the following, we exemplify the calculation of some relevant event propensities.

##### 2.1. Reaction Events

, where denotes the number of possible molecular combinations available for reaction at time t.

For first- and second-order reactions, the stochastic rate constants are calculated from the macroscopic reaction rate constants as for monomolecular reactions and as for bimolecular reactions given the system volume .

##### 2.2. Diffusion Events

, where denotes the number of molecules of the diffusing species at time .

, where denotes the specific diffusion coefficient and the diffusion area.

##### 2.3. Channel Gating

Channel gating, that is, subunit transitions, are naturally described by a discrete stochastic process. Their biological properties force us to differentiate between ligand-gated and voltage-gated ion channels.

The interaction of ligands and channel subunits can be modeled as ordinary reaction events. Keeping track of the subunit states and depending on the gating scheme, we can decide whether the channel is in an open state or not.

The introduction of voltage-gated ion channels is more complex and requires a characteristic function, representing the membrane potential. The transition rates between the open/close states of the channels are based on the actual function value.

##### 2.4. Channel Conductance

Ion permeation through a channel protein is modeled by translating the electrical conductance to a rate constant which can then be transformed into a permeation probability.

, where is the elementary charge and is the charge of the permeating ion. denotes the mean current of the channel in .

All further considerations assume a well-stirred chemical system, and the velocities and positions of individual molecules in are ignored. This assumption is based on the fact that nonreactive molecular collisions that lead to a mixing of the system are far more likely than the occurrence of reactive events. This leads to uniformly randomized positions and thermally randomized velocities throughout the system volume (Maxwell-Boltzmann distribution) [7].

All simulation strategies can finally be derived from the chemical master equation (CME). The CME is a differential (or difference) equation that describes the time evolution of the probability density that is the multivariate distribution of the state vector . Using the notation introduced above, the CME reads In words, (1) states that the change in is calculated as the net probability flow conveyed by the flows from state to (via event ) and the reverse flows out of state .

#### 3. The Gillespie Algorithm

While the simulation of the Brownian motion of all individual molecules would provide an accurate model of Ca^{2+} dynamics, the computational effort is very high to examine signaling pathways on an adequate spatiotemporal scale. Increasing system volume and increasing numbers of molecules lead to extreme increases in computation time. Gillespie's algorithm is an alternative to accurately simulate small system volumes under the assumption of a well stirred system. Sample paths generated using this approach are samples of a multivariate Markov process. This means that given a system state , the state at any future time only depends on and not on the previous history of the system. Transitions between system states are performed by generating independent, exponentially distributed waiting times for each possible event (according to the event propensities) and choosing the event with the minimum waiting time. The computational procedure is the same for all event types . It uses the current event propensity and a uniformly distributed random variable :

Equation (2) leads to exponentially distributed waiting times which are characteristic for Markov processes. After selection of the event with the smallest -value , the current simulation time is advanced to , and the system state is changed to . To improve the computational performance, the Gibson-Bruck algorithm introduced a dependency graph which leads to an efficient update rule where only those reaction propensities are updated that are influenced by the last reaction realized. A reaction depends on a reaction if one of the educts of is changed by reaction [8]. A complete overview of improvements of the classical Gillespie algorithm can be found in the literature [9].

The algorithm reads as follows.(0)*Initialization*: to initialize the system, set the initial number of each species to , , and calculate the waiting times for each event . (1)*Event selection*: select the event that minimizes the associated -value (2). (2)*Update time and state*: increment the system time and change the state vector according to . (3)*Exit condition*: if then *Exit*, where is the maximum simulation time. (4)*Update*: recalculate -values of related events as stored in the dependency graph. (5)*Go to step 1*.

The Gillespie algorithm quickly meets its limits when systems with increasing system volumes and large numbers of implemented events are simulated. In particular, the presence of different time scales can decrease the efficiency crucially (e.g., diffusion events outnumber Ca^{2+} buffer association/dissociation events by far [10]). The -leaping method, which is introduced in the next paragraph, partially solves these problems.

#### 4. The -Leaping Method

Frequently, it is not necessary to track all individual events explicitly because event propensities do not change significantly in a certain time interval of length . Using the classical Gillespie algorithm, this fact leads to a high computational effort which is not always necessary in terms of accuracy. The -leaping method takes advantage of this observation by evaluating how many times each event will occur in the discrete time interval . To apply the -leaping method to model systems, it is necessary to find an appropriate value for every simulation step. The so-called leaping condition requires the time interval to be small enough to assure that no propensity function undergoes an appreciable change. This is traditionally accomplished by bounding to the product of a defined error control parameter and the sum of all propensity functions [9]: Considering chemical reaction systems exhibiting a wide range of temporal scales, this approach might lead to inaccuracies due to the high variability of the involved event propensities. To satisfy the leap condition in such cases, Cao et al. [11] suggested a new -selection procedure which bounds the relative changes in by a fixed factor . A detailed discussion of different -selection procedures can be found in the literature [11].

In this context, denotes the average occurrence of in , which is adequately approximated by a Poisson-distributed random variable . This approach leads to the -leaping approximation which reads Poisson random variables are unbounded and some reaction networks can produce the impossible scenario of negative copy numbers for certain chemical species. Especially, small reactant populations with copy numbers close to zero are affected by this. Different solutions have been proposed to solve this problem. Tian and Burrage [12] developed an alternative by replacing the Poisson-distributed random variable by a binomial random variable which can be bounded by the actual number of reactants of each species.

Cao et al. [11] introduced the modified (nonnegative) Poisson -leaping approach. It is based on the idea of subdividing the set of all events into critical and noncritical events depending on their risk to induce negative reaction species populations in . Simulating critical events with the classical Gillespie algorithm and non-critical events using -leaping, the probability of negative molecule counts becomes nearly zero.

#### 5. The Chemical Langevin Equation

The approximation of the event count occurring in a carefully chosen time interval can be taken a step further. If a chemical reaction system not just fulfills (a) the leaping condition but also satisfies (b) , the Poisson-distributed random variable can be approximated by a normally (Gaussian-) distributed random variable with mean and variance equal to : With the linear combination theorem of random variables and (4), we can derivate the “white noise form” of the Langevin equation: where at denotes a temporally uncorrelated Gaussian white noise process (the increments of a Brownian motion ).

The CLE approach implies some notable consequences. First, the integration time interval is fixed. Second, due to the transformation of the Poisson to a Gaussian random variable, the resulting molecule counts become real values rather than integers. The change of the state vector is dependent on a deterministic part (first sum of (6)), and a stochastic part (second sum of (6)). When the number of reactants tends to infinity, the stochastic term can be neglected compared to the deterministic part, and (6) reduces to the deterministic *reaction rate equation* approach. This extrapolation provides the link between conventional deterministic chemical kinetics and the stochastic approach. A complete derivation of the CLE can be found in the literature [13].

Furthermore, the theory of stochastic processes allows a transformation of the Langevin equation to an associated *Fokker-Planck equation* (FPE) that describes the temporal evolution of the probability density [14]. The FPE method has been applied to calcium dynamics in the dyadic cleft of cardiomyocytes [15]. Given that researchers are often more interested in explicit sample paths of a model than in probability distributions, the CLE approach is more common.

#### 6. Hybrid Simulations

The simulation of complex chemical systems, including a set of events with a high variability of propensity functions, often cannot be conveniently approximated by the CLE. The resulting loss in accuracy ignores important events on slow temporal scales that affect the evolution of the whole system (e.g., ion channels). A promising way to avoid the high computational effort of the classical Gillespie algorithm while achieving a satisfying degree of accuracy is *Hybrid simulation algorithms* that combine deterministic and stochastic approaches.

Rüdiger et al. [16] successfully proposed a hybrid stochastic and deterministic approach for the simulation of Ca^{2+} blips, describing the state transitions between the subunits of an IP_{3}R calcium channel with the classical Gillespie algorithm and the surrounding Ca^{2+}-dynamics (buffering, diffusion) as a deterministic process. Similar approaches have been introduced by Stern et al. [17] and Greenstein and Winslow [18] using cardiac and skeletal muscle cells as model systems. Kalantzis [19] and Choi et al. [20] proposed hybrid algorithms that switch adaptively between the classical Gillespie algorithm, the -leaping method, and the chemical Langevin equation. Skupin et al. introduced a spatially resolved multiscale model which combines detailed stochastic channel gating on the scale of a channel cluster with a linearized spatial bidomain model to integrate them on a microscopic scale [21].

#### 7. Conclusion

High-resolution confocal laser microscopy and mathematical models of calcium signaling showed that stochastic effects can be essential for intracellular information processing. At the same time, we observe a continuous improvement in data quality through advanced measurement techniques and increasing accuracy and availability of physical and chemical properties of signaling molecules. Thus, stochastic simulation approaches are increasingly used to simulate subcellular signaling pathways. New tendencies in developing integrated simulation algorithms that case-dependently switch between different approaches are very promising and may provide efficient solutions for the simulation of large systems with complex event interactions. Most importantly, future work has to prove the validity of those approaches on different spatiotemporal scales.

The degree of stochasticity necessary to model Ca^{2+} microdomains is still discussed. Tanskanen et al. [22] as well as Hake and Lines [23] both proposed two different models of the cardiac dyadic cleft, using different algorithms. The dyadic cleft is a functionally important microdomain in cardiac myocytes where the process of excitation-contraction coupling takes place. This microdomain is located between the sarcolemma and the endoplasmic reticulum membrane and is characterized by steep Ca^{2+}-concentration gradients and atoliter volumes. The interaction between L type calcium channels (LCC) and clusters of Ryanodine receptor calcium channels is the key event of this process. Tanskanen et al. [22] introduced a model using a spatially resolved Markov process, taking into account the positions of individual Ca^{2+}-ions, dyadic proteins, the membrane surface charges, and channel gating [22]. They conclude that the full stochasticity of their approach is necessary for an exact description of the dyad. On the other hand, Hake and Lines [23] also used the dyadic cleft as the structural basis of their work. Comparing a hybrid model, combining a deterministic continuous model with stochastic receptor gating and a fully stochastic *Random Walk* (RW) approach, they conclude that under certain limitations of the simulation parameters, the deterministic approach reproduces the RW results sufficiently well. Thurley and Falcke [24] recently used a hybrid simulation approach to study the relation of robustness and sensitivity of calcium-dependent subcellular pathways based on the statistical properties of interspike intervals. Even though there is a high diversity of deterministic, stochastic, and hybrid simulation strategies, the specifying system parameters are equal. The open source community has consistently developed the *System Biology Markup Language* (SBML) [25] which defines a universal *XML* file format to ensure the interchangeability of model definitions between different software packages. Therefore, one of the major goals of future efforts should be a support of the SBML or an equivalent file format to help with the advance of a flexible, barrier-free system biology community.

#### References

- A. Tsugorka, E. Rios, and L. A. Blatter, “Imaging elementary events of calcium release in skeletal muscle cells,”
*Science*, vol. 269, no. 5231, pp. 1723–1726, 1995. View at Google Scholar · View at Scopus - S. A. Stricker, “Repetitive calcium waves induced by fertilization in the nemertean worm Cerebratulus lacteus,”
*Developmental Biology*, vol. 176, no. 2, pp. 243–263, 1996. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - R. E. Dolmetsch, K. Xu, and R. S. Lewis, “Calcium oscillations increase the efficiency and specificity of gene expression,”
*Nature*, vol. 392, no. 6679, pp. 933–936, 1998. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - G. Dupont, A. Abou-Lovergne, and L. Combettes, “Stochastic aspects of oscillatory Ca2+ dynamics in hepatocytes,”
*Biophysical Journal*, vol. 95, no. 5, pp. 2193–2202, 2008. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - G. Dupont and L. Combettes, “Modelling the effect of specific inositol 1,4,5-trisphosphate receptor isoforms on cellular Ca2+ signals,”
*Biology of the Cell*, vol. 98, no. 3, pp. 171–182, 2006. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - D. T. Gillespie, “Stochastic simulation of chemical kinetics,”
*Annual Review of Physical Chemistry*, vol. 58, pp. 35–55, 2007. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,”
*Journal of Physical Chemistry*, vol. 81, no. 25, pp. 2340–2361, 1977. View at Google Scholar - M. A. Gibson and J. Bruck, “Efficient exact stochastic simulation of chemical systems with many species and many channels,”
*Journal of Physical Chemistry A*, vol. 104, no. 9, pp. 1876–1889, 2000. View at Google Scholar - D. T. Gillespie, “Approximate accelerated stochastic simulation of chemically reacting systems,”
*Journal of Chemical Physics*, vol. 115, no. 4, pp. 1716–1733, 2001. View at Publisher · View at Google Scholar - F. Von Wegner and R. H. A. Fink, “Stochastic simulation of calcium microdomains in the vicinity of an L-type calcium channel,”
*European Biophysics Journal*, vol. 39, no. 7, pp. 1079–1088, 2010. View at Publisher · View at Google Scholar · View at PubMed - Y. Cao, D. T. Gillespie, and L. R. Petzold, “Efficient step size selection for the tau-leaping simulation method,”
*The Journal of chemical physics*, vol. 124, no. 4, Article ID 044109, 2006. View at Publisher · View at Google Scholar · View at PubMed - T. Tian and K. Burrage, “Binomial leap methods for simulating stochastic chemical kinetics,”
*Journal of Chemical Physics*, vol. 121, no. 21, pp. 10356–10364, 2004. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - D. T. Gillespie, “Chemical Langevin equation,”
*Journal of Chemical Physics*, vol. 113, no. 1, pp. 297–306, 2000. View at Publisher · View at Google Scholar · View at Scopus - C. Gardiner,
*Handbook of Stochastic Methods: For Physics, Chemistry and the Natural Sciences*, Springer, 3rd edition, 2004. - R. L. Winslow, A. Tanskanen, M. Chen, and J. L. Greenstein, “Multiscale modeling of calcium signaling in the cardiac dyad,”
*Annals of the New York Academy of Sciences*, vol. 1080, pp. 362–375, 2006. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - S. Rüdiger, J. W. Shuai, W. Huisinga et al., “Hybrid stochastic and deterministic simulations of calcium blips,”
*Biophysical Journal*, vol. 93, no. 6, pp. 1847–1857, 2007. View at Publisher · View at Google Scholar · View at PubMed - M. D. Stern, G. Pizarro, and E. Ríos, “Local control model of excitation-contraction coupling in skeletal muscle,”
*Journal of General Physiology*, vol. 110, no. 4, pp. 415–440, 1997. View at Publisher · View at Google Scholar · View at Scopus - J. L. Greenstein and R. L. Winslow, “An integrative model of the cardiac ventricular myocyte incorporating local control of Ca2+ release,”
*Biophysical Journal*, vol. 83, no. 6, pp. 2918–2945, 2002. View at Google Scholar · View at Scopus - G. Kalantzis, “Hybrid stochastic simulations of intracellular reaction-diffusion systems,”
*Computational Biology and Chemistry*, vol. 33, no. 3, pp. 205–215, 2009. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - T. Choi, M. R. Maurya, D. M. Tartakovsky, and S. Subramaniam, “Stochastic hybrid modeling of intracellular calcium dynamics,”
*Journal of Chemical Physics*, vol. 133, no. 16, Article ID 165101, 2010. View at Publisher · View at Google Scholar · View at PubMed - A. Skupin, H. Kettenmann, and M. Falcke, “Calcium signals driven by single channel noise,”
*PLoS Computational Biology*, vol. 6, no. 8, Article ID e1000870, 2010. View at Publisher · View at Google Scholar · View at PubMed - A. J. Tanskanen, J. L. Greenstein, A. Chen, S. X. Sun, and R. L. Winslow, “Protein geometry and placement in the cardiac dyad influence macroscopic properties of calcium-induced calcium release,”
*Biophysical Journal*, vol. 92, no. 10, pp. 3379–3396, 2007. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - J. Hake and G. T. Lines, “Stochastic binding of Ca2+ ions in the dyadic cleft; Continuous versus random walk description of diffusion,”
*Biophysical Journal*, vol. 94, no. 11, pp. 4184–4201, 2008. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus - K. Thurley and M. Falcke, “Derivation of Ca2+ signals from puff properties reveals that pathway function is robust against cell variability but sensitive for control,”
*Proceedings of the National Academy of Sciences of the United States of America*, vol. 108, no. 1, pp. 427–432, 2011. View at Publisher · View at Google Scholar · View at PubMed - M. Hucka, A. Finney, B. J. Bornstein et al., “Evolving a lingua franca and associated software infrastructure for computational systems biology: the Systems Biology Markup Language (SBML) project,”
*Systems biology*, vol. 1, no. 1, pp. 41–53, 2004. View at Publisher · View at Google Scholar · View at Scopus