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.


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.

Background
Every chemical system is defined by a set of properties. First of all, there are N chemical species, S = {S 1 , . . . , S N }, in the system volume Ω. The state vector x(t) = (x 1 (t), . . . , x N (t)) denotes the number of molecules x k (t), k = 1, . . . , N of each species S k at time t 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 a j by its stochastic rate constant c j and the stoichiometric coefficients of the underlying process. A chemical system contains M events R = {R 1 , . . . , R M } which are represented through a state change vector v j = (d j,1 , . . . , d j,N ), j = 1, . . . , M, where component d j,i denotes the change in the number of molecules S i due to the event R j . Even more complex processes, which do not affect x(t), 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.

Reaction Events
denotes the number of possible molecular combinations available for reaction R j at time t.
For first-and second-order reactions, the stochastic rate constants c j are calculated from the macroscopic reaction rate constants k j as c j = k j for monomolecular reactions and as c j = k j /Ω for bimolecular reactions given the system volume Ω.

Diffusion
Events. a k (x, t) = c k x k (t), where x k (t) denotes the number of molecules of the diffusing species k at time t. c k = D k /A, where D k denotes the specific diffusion coefficient and A the diffusion area.

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.

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 e is the elementary charge and z ion is the charge of the permeating ion. I ch denotes the mean current of the channel in C/s.
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 P(x, t) that is the multivariate distribution of the state vector x(t). Using the notation introduced above, the CME reads (1) In words, (1) states that the change in P(x, t) is calculated as the net probability flow conveyed by the flows from state x − v j to x (via event R j ) and the reverse flows out of state x.

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 x(t), the state x(t + τ) at any future time t + τ only depends on x(t) and not on the previous history of the system. Transitions between system states are performed by generating independent, exponentially distributed waiting times τ j for each possible event R j (according to the event propensities) and choosing the event with the minimum waiting time. The computational procedure is the same for all event types R j . It uses the current event propensity a j (x, t) and a uniformly distributed random variable r ∼ U[0, 1]: Equation (2) leads to exponentially distributed waiting times τ which are characteristic for Markov processes. After selection of the event R j with the smallest τ-value τ j , the current simulation time t is advanced to t ← t + τ j , and the system state is changed to x(t) ← x(t) + v j . 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 R b depends on a reaction R a if one of the educts of R b is changed by reaction R a [8]. A complete overview of improvements of the classical Gillespie algorithm can be found in the literature [9].
Journal of Biomedicine and Biotechnology 3 The algorithm reads as follows.
(0) Initialization: to initialize the system, set the initial number of each species S i to x(t) = x(t 0 ), t = t 0 , and calculate the waiting times τ j , j = 1, . . . , M for each event R j .
(1) Event selection: select the event R j that minimizes the associated τ-value (2).
(2) Update time and state: increment the system time t ← t + τ j and change the state vector according to R j : (3) Exit condition: if (t + τ j > T max ) then Exit, where T max is the maximum simulation time.
(4) Update: recalculate τ-values of related events as stored in the dependency graph.
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.

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 R j will occur in the discrete time interval [t, t + τ). To apply the τleaping method to model systems, it is necessary to find an appropriate τ value for every simulation step. The socalled leaping condition requires the time interval to be small enough to assure that no propensity function a j undergoes an appreciable change. This is traditionally accomplished by bounding Δ τ a j (x, t) to the product of a defined error control parameter and the sum of all propensity functions a 0 [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 a j by a fixed factor . A detailed discussion of different τ-selection procedures can be found in the literature [11]. In this context, a j (x)τ denotes the average occurrence of R j in [t, t + τ), which is adequately approximated by a Poisson-distributed random variable ξ j (a j (x), τ). 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 [t, t + τ). Simulating critical events with the classical Gillespie algorithm and non-critical events using τ-leaping, the probability of negative molecule counts becomes nearly zero.

The Chemical Langevin Equation
The approximation of the event count occurring in a carefully chosen time interval [t, t + τ) can be taken a step further. If a chemical reaction system not just fulfills (a) the leaping condition but also satisfies (b) a j τ 1, for all j ∈ [1, M], the Poisson-distributed random variable ξ j can be approximated by a normally (Gaussian-) distributed random variable with mean and variance equal to a j (x)τ: With the linear combination theorem of random variables and (4), we can derivate the "white noise form" of the Langevin equation: where at dB(t) denotes a temporally uncorrelated Gaussian white noise process (the increments of a Brownian motion B(t)). The CLE approach implies some notable consequences. First, the integration time interval dt 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 x(t) 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 P(x, t) [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.

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].

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.