Exact and Approximate Stochastic Simulation of Intracellular Calcium Dynamics

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²⁺ 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²⁺ dynamics. Oscillations in the typically low-intracellular Ca²⁺ 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²⁺ release events from intracellular Ca²⁺ stores (sarcoplasmic and endoplasmatic reticulum) [4]. However, the link between the dynamics of individual molecules involved in microdomain Ca²⁺ dynamics (e.g., Ca²⁺ 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²⁺ ions are the reason for its well-described chemical and physical properties. Those circumstances make the in silico examination of Ca²⁺ dynamics very promising. Microdomains controlled by Ca²⁺ channels play an important role in the context of Ca²⁺ 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 𝑥_𝑘(𝑡),𝑘=1,…,𝑁 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 𝑐_j and the stoichiometric coefficients of the underlying process. A chemical system contains 𝑀 events 𝑅={𝑅₁,…,𝑅_𝑀} which are represented through a state change vector 𝐯_𝑗=(𝑑_𝑗,1,…,𝑑_𝑗,𝑁),𝑗=1,…,𝑀, 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.

𝑐_𝑗=𝐼_ch/(𝑒⋅𝑧_ion), where 𝑒 is the elementary charge and 𝑧_ion is the charge of the permeating ion. 𝐼_ch 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𝜕_𝑡𝑃(𝐱,𝑡)=^𝑀_𝑗=1𝑎_𝑗𝐱−𝐯_𝑗,𝑡𝑃𝐱−𝐯_𝑗,𝑡−𝑎_𝑗(𝐱,𝑡)𝑃(𝐱,𝑡). (1) 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²⁺ 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 𝑟∼𝑈[0,1]:𝜏_𝑗=−log𝑟𝑎_𝑗().𝐱,𝑡 (2)

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 𝜏_𝑗,𝑗=1,…,𝑀 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 (𝑡+𝜏_𝑗>𝑇_max) then Exit, where 𝑇_max 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²⁺ 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]:Δ_𝜏𝑎_𝑗(𝐱,𝑡)≤𝜖⋅𝑎₀(𝐱,𝑡). (3) 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𝐱(𝑡+𝜏)=𝐱(𝑡)+^𝑀_𝑗=1𝐯_𝑗𝜉_𝑗𝑎_𝑗(𝐱)𝜏. (4) 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) 𝑎_𝑗𝜏≫1,forall𝑗∈[1,𝑀], the Poisson-distributed random variable 𝜉_𝑗 can be approximated by a normally (Gaussian-) distributed random variable with mean and variance equal to 𝑎_𝑗(𝐱)𝜏:𝜉_𝑗𝑎_𝑗(𝐱)𝜏≈𝑁𝑎_𝑗(𝐱)𝜏,𝑎_𝑗(𝐱)𝜏. (5) With the linear combination theorem of random variables and (4), we can derivate the “white noise form” of the Langevin equation:𝐱(𝑡+d𝑡)=𝐱(𝑡)+^𝑀_𝑗=1𝐯_𝑗𝑎_𝑗(𝐱(𝑡))d𝑡+^𝑀_𝑗=1𝐯_𝑗𝑎_𝑗(𝐱(𝑡))d𝐵(𝑡), (6) where at d𝐵(𝑡) 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 d𝑡 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²⁺ blips, describing the state transitions between the subunits of an IP₃R calcium channel with the classical Gillespie algorithm and the surrounding Ca²⁺-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²⁺ 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²⁺-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²⁺-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.