Microdomain calcium fluctuations as a colored noise process

Calcium ions play a key role in subcellular signaling as localized transients of the intracellular calcium concentration modify the activity of ion channels, enzymes and transcription factors, among others. The intracellular calcium concentration is inherently noisy, as diffusion, the transient binding to and dissociation from buffer molecules and stochastically gating calcium channels contribute to the fluctuations of the local copy number of Ca²⁺ ions. We study the properties of the fluctuating calcium concentration in sub-femtoliter volumes using an exact stochastic simulation algorithm and approximations to the exact stochastic solution. It is shown that the time course of the local calcium concentration represents a colored noise process whose autocorrelation time is a function of buffer kinetics and diffusion constants. Using the chemical Langevin description and the excess buffer approximation of the process, fast approximative algorithms and theoretical connections to the Ornstein-Uhlenbeck process are obtained. In a generic example, we show how calcium noise can couple to the dynamics of a single variable moving in a double-well potential, leading to a colored noise induced transition. Our work shows how a multitude of intracellular signaling pathways may be influenced by the inherent stochasticity of calcium signals, a key messenger in virtually any cell type, and how the calcium signal can be implemented efficiently in cellular signaling models.

1. Introduction

Calcium ions regulate intracellular signaling pathways in virtually all known cell types as many regulating proteins offer calcium binding sites by exhibiting negatively charged amino acid residues. Among these proteins, we can find calcium-regulated ion channels, enzymes and transcription factors, some of them clustered in families such as the EF-hand proteins (Bhattacharya et al., 2004). Due to the strong compartmentalization of the intracellular space, the local calcium concentration can differ significantly between two spots within the same cell. The most obvious example may be that of a small intracellular domain next to a calcium channel located within the plasma membrane or the endoplasmic reticulum membrane. Given the strong extra-/intracellular concentration gradient for calcium, channel opening leads to a quick rise of the local calcium concentration. The spatiotemporal extent of this calcium increase is determined by the kinetic constants and the diffusibility of the buffers involved (Smith et al., 1998; Jiang et al., 1999; Uttenweiler et al., 2002). Moreover, the non-linear interplay of different calcium release und re-uptake mechanisms generate more complex global calcium patterns such as traveling waves and oscillations. All of the aforementioned phenomena can be modeled deterministically using ordinary and partial differential equations. However, closer inspection of the conditions encountered in calcium microdomains shows that the subcellular calcium concentration must contain a certain amount of stochasticity as induced by the stochastic nature of diffusion, transient chemical binding and the gating of calcium channel proteins. The magnitude of the fluctuations, and therefore their possible physiological relevance, depends on the reaction volumes considered. Ignoring contributions from calcium channels and given an intracellular calcium concentration of approximately 100 nM in a resting cell, we expect to find approximately 60 Ca²⁺ ions in a volume of 10⁻¹⁵ liter (1 fl). Fluctuations around the mean value will have a SD given by the square root of the mean value, i.e., the number of Ca²⁺ ions will be approximately $60\pm \sqrt{60}$ ions. Femtoliter volumes are of special interest for several reasons. First, localized calcium transients such as calcium sparks, puffs and quarks occupy femtoliter volumes, as well as their functional targets such as ion channel clusters and mitochondria (Cheng and Lederer, 2008). Second, current measurement devices for calcium fluorescence microsopy, such as confocal and multiphoton laser microscopes, collect fluorescence from volumes of approximately 1 fl. Therefore, computational models of reactions occurring within these volumes yield important data that can help to design new experiments and to interpret experimental results.

Until recently, partial differential equations were the main tool to build reaction-diffusion models of calcium microdomains (Smith et al., 1998; Jiang et al., 1999; Uttenweiler et al., 2002). The PDE approach is completely deterministic and ignores deviations from equilibrium concentrations and deterministic solutions. In a previous paper, we have shown computationally that fluctuations of calcium and other reactants lead to highly variable responses of calcium-sensitive signaling pathways (von Wegner and Fink, 2010), suggesting an important role of the stochastic aspects of calcium dynamics. Other authors have studied calcium signaling on a whole cell level using approximations of the exact stochastic model (Zhang et al., 2004; Li et al., 2005; Manninen et al., 2006; Zhu et al., 2007; Choi et al., 2010). The use of different stochastic models in the context of microdomain calcium signaling has been summarized in recent review articles (Wieder et al., 2011; von Wegner et al., 2012). In the present paper, we study the properties of the fluctuating calcium concentration in sub-femtoliter volumes using Gillespie's exact stochastic simulation algorithm as well as stochastic approximations. It is shown that the time course of the local calcium concentration represents a colored noise process whose noise color, i.e., its autocorrelation time is a function of the kinetic constants of the buffer species considered. Using the chemical Langevin equation, theoretical connections to the Ornstein-Uhlenbeck process are obtained. A possible role of calcium noise color in the gating behavior of a calcium-sensitive ion channel is illustrated using a colored noise driven bistable dynamical system as a generic representation of a calcium noise-driven ion channel. Our results show how a multitude of intracellular signaling pathways may be influenced by the inherent stochasticity of calcium signals, a key messenger in virtually any cell type, and how the stochastic calcium signal can be implemented efficiently in computational models of cellular calcium signaling.

2. Deterministic and Stochastic Description of Calcium Dynamics

In the following, we give a mathematical description of deterministic and stochastic calcium dynamics. We will denote the time-dependent concentrations of calcium ions, free buffer molecules and calcium-buffer complexes as [Ca²⁺]_t, [B]_t, and [CaB]_t, respectively. For the sake of readability, we substitute these bracketed expressions by x_t, y_t, and z_t in formulas. When a stochastic framework is used, we will work with molecule counts rather than concentrations and denote these as capital letters, e.g., the number of free calcium ions in a given volume V is denoted as X_t = x_t × V (analogously for Y_t and Z_t). Equilibrium concentrations (or molecule counts) are denoted by a tilde, e.g., the free calcium concentration at chemical equilibrium is written as $\tilde{x}$, and the corresponding number of calcium ions at equilibrium is denoted $\tilde{X}$. When a system contains M different buffer species B₁, …, B_M, we write their time-dependent concentrations with parenthesized superscripts y⁽¹⁾_t, …, y^(M)_t, and their molecule counts as Y⁽¹⁾_t, …, Y^(M)_t. Initial concentrations receive a 0-subscript, e.g., x₀ is the initial concentration of free calcium ions at t = 0. The total calcium concentration is [Ca]_T = x_T = [Ca²⁺] + [CaB], and the total buffer concentration is B_T = y_T = [B] + [CaB].

2.1. Deterministic Dynamics

Consider a simple reaction system containing calcium ions, Ca²⁺, and a single calcium buffer B. In real biological systems and in experimental settings, the buffer often is a protein, calmodulin for instance, or a small organic molecule, e.g., EGTA. The following reaction scheme describes the system :

$\begin{array}{ll}{\text{Ca}}^{2+}+\text{B}\stackrel{k^{+}}{\underset{k^{-}}{\rightleftharpoons }}\text{CaB}\hfill & (1)\hfill \end{array}$

with reaction velocities v⁺_t = k⁺[Ca²⁺]_t[B]_t and v⁻ = k⁻[CaB]_t, defined via the reaction rate constants k⁺, k⁻.

The equilibrium condition $k^{+}[\widetilde{{\text{Ca}}^{2+}}][\tilde{\text{B}}]=k^{-}[\widetilde{\text{CaB}}]$ leads to:

$\begin{array}{ll}[\widetilde{\text{CaB}}]=\frac{{[\text{B}]}_T}{1+\frac{K_D}{[\widetilde{{\text{Ca}}^{2+}}]}}.\hfill & (2)\hfill \end{array}$

Thus, all equilibrium concentrations depend solely on the total calcium and buffer concentrations x_T and y_T, the buffer dissociation constant $K_D=\frac{k^{-}}{k^{+}}$, and the nominal free calcium concentration at equilibrium ($[\widetilde{{\text{Ca}}^{2+}}]=\tilde{x}$).

The following ordinary differential equation describes the deterministic kinetics of the free calcium concentration:

$\begin{array}{ll}{\dot{x}}_t=-k^{+}{x}_t{y}_t+k^{-}{z}_t.\hfill & (3)\hfill \end{array}$

Given fixed total concentrations [Ca²⁺]_T and [B]_T, the single buffer system has only one independent variable, and the dynamics of the other two variables is given by ẏ_t = ẋ_t and ż_t = −ẋ_t. Using the definitions of y_T and x_T, we can rewrite the dynamics as the non-linear ODE:

$\begin{array}{l}{\dot{x}}_t=-k^{+}{x}_t[y_T-x_T+x_t]+k^{-}[x_T-x_t]\\ =-k^{+}\left(x_t[x_t+y_T-x_T+K_D]-K_D{x}_T\right).\end{array}$

We will apply the excess buffer approximation to this system to obtain a linear ODE, and in subsequent sections derive considerably simplified stochastic versions of the kinetic equations.

In the present work, we have used the physico-chemical properties of Ca²⁺ ions and buffers as published elsewhere and as given in Table 1, along with references. We provide the variable simulation parameters such as reaction volumes in the corresponding results sections.

TABLE 1

Table 1. Physiological and artificial calcium buffers used in simulations.

2.2. The Excess Buffer Approximation (EBA)

In biological systems, at rest, the free calcium concentration is approximately 100 nM, whereas the total buffer concentration lies in the micromolar range (Smith et al., 1998; Jiang et al., 1999; Uttenweiler et al., 2002; Novo et al., 2003; von Wegner and Fink, 2010). The free buffer concentration can therefore often be assumed to be constant, i.e., y_t ≈ ỹ. In biological systems, the total buffer concentration is usually large enough to assume a constant free buffer concentration [B], i.e., only a small fraction of the total buffer molecules is bound to Ca²⁺. Thus, under the EBA assumption (Heinemann et al., 1986; Smith et al., 2001; Fall et al., 2002), the free buffer concentration [B] is assumed to be identical to the equilibrium concentration, i.e., y_t = ỹ. In the original derivation of the EBA, calcium and buffer diffusion were considered close to a calcium channel pore. The approximation was considered valid for small ratios K_D/[B]_T, i.e., when the total buffer concentration is [B]_T large, or when dealing with high-affinity buffers. In the current work, additional calcium influx through a calcium channel is not considered, and therefore, all systems are far from buffer saturation. Using the approximation y_T = ỹ + z_t, valid under the EBA assumption, the simplified deterministic dynamics can be rewritten as:

$\begin{array}{ll}\begin{array}{l}{\dot{x}}_t=-k^{+}{x}_t\tilde{y}+k^{-}[x_T-x_t]\\ =-(k^{+}\tilde{y}+k^{-})x_t+k^{-}{x}_T.\end{array}\hfill & (4)\hfill \end{array}$

The free calcium concentration [Ca²⁺] follows an inhomogeneous ODE with constant coefficients

$\begin{array}{ll}\begin{array}{l}{K}_1=k^{+}\tilde{y}+k^{-}\\ K_2=k^{-}{x}_T\end{array}\hfill & (5)\hfill \end{array}$

and with initial condition [Ca²⁺]₀ = x₀:

$\begin{array}{ll}{\dot{x}}_t=-K_1{x}_t+K_2.\hfill & (6)\hfill \end{array}$

The equilibrium calcium concentration $\tilde{x}$ in terms of the new ODE coefficients is $\tilde{x}=\frac{K_2}{K_1}$. The straightforward solution to Equation 6 is

$\begin{array}{ll}{x}_t=(x_0-\tilde{x})\exp (-K_{1}t)+\tilde{x}.\hfill & (7)\hfill \end{array}$

Clearly, the EBA-simplified dynamics predicts a mono-exponential relaxation of the single buffer system to calcium elevations that do not violate the EBA-assumption, in other words, the Ca²⁺ peak must not alterate the free buffer concentration significantly.

2.3. Stochastic Dynamics

In this section, we introduce stochastic descriptions of the simplified calcium dynamics. We first review Gillespie's original algorithm for exact stochastic simulations without re-deriving all details. Extended presentations of the algorithm for general systems and for calcium dynamics in particular can be found in the literature (Gillespie, 1977, 2007; von Wegner and Fink, 2010; Wieder et al., 2011; von Wegner et al., 2012). Subsequently, we formulate the chemical Langevin equation for the model systems at hand, i.e., a single buffer system, a single buffer system with calcium diffusion, and a multi-buffer system without diffusion. Last, we introduce the Ornstein-Uhlenbeck process as the generic Gaussian, exponentially correlated colored noise process.

2.3.1. Exact stochastic simulation algorithm—Gillespie's algorithm (SSA)

Systems of chemical reactions can be represented by multivariate Markov processes (Gardiner, 2009). In the case of chemical reaction systems, each variable describes the time-dependent copy number of exactly one molecular species. The Gillespie algorithm generates exact sample paths of the Markov process alternating between two sampling steps. First, a sample of the state-dependent waiting time distribution until the next reaction event is generated. In a second step, exactly one reaction is selected from the list of all possible reactions according to their relative reaction propensities (Gillespie, 2000, 2007). We will restrict the description of the algorithm to the notions that will be used in the following sections. In order to correctly capture reaction probabilities in small volumes, the deterministic reaction rate constants k⁺, k⁻ have to be transformed into the corresponding stochastic rate constants c⁺, c⁻. For mono- and bimolecular reactions, the transformations are elementary and only involve a scaling by the reaction volume V for the bimolecular reaction: $c^{+}=\frac{k^{+}}{V}$, c⁻ = k⁻, using the notation used in Equation 1. As we will exclusively deal with systems of mono- and bimolecular reactions, all non-vanishing stoichiometric coefficients v_ji are either −1 or +1. The coefficient v_ji reflects the change in the copy number of reactant i due to the reaction indexed by j. Diffusion is implemented as derived in Elf et al. (2003) and as implemented for calcium microdomains in von Wegner and Fink (2010). As the waiting time to the next reaction event is a random variable, Gillespie's algorithm generates non-equidistant sample paths that are linearly interpolated prior to further processing. We chose an interpolation interval of dt = 0.01 ms that is identical to the integration time step dt used for the chemical Langevin equation and the Ornstein-Uhlenbeck process.

2.3.2. The chemical Langevin equation (CLE)

The CLE approach combines the deterministic dynamics as described by reaction rates with the stochastic description of the Gillespie algorithm. Mathematically, it is based on the approximation of a Poisson-distributed random variable, representing the number of reaction events in a short interval of time, by a normally distributed random variable, thus achieving the connection to the general Langevin equation (Gillespie, 2000). The dynamics are driven by a stochastic process dW_t. The terms dW_t represent the stochastic increments of a Brownian motion (also called a standard Wiener process) (Gardiner, 2009). The increments are normally distributed as W_t+dt − W_t ~ (0, dt) and are pairwise uncorrelated. Furthermore, the initial condition is W₀ = 0. Algorithmically, such increments are computed using normally distributed pseudo-random numbers with mean μ = 0 and variance σ² = dt. Equivalently, pseudo-random numbers drawn from a standard normal distribution ((0, 1)) can be multiplied by the factor $\sqrt{dt}$ to obtain the same process. Explicit examples for systems as discussed here are given in von Wegner et al. (2012); Higham (2001).

For the single buffer system, Equation 1, the chemical Langevin equation reads:

$\begin{array}{ll}\begin{array}{l}d{X}_t=-c^{+}{X}_t{Y}_{t}dt+c^{-}{Z}_{t}dt-\sqrt{c^{+}{X}_t{Y}_t}d{W}_t^{(1)}\\ +\sqrt{c^{-}{Z}_t}d{W}_t^{(2)}\end{array}\hfill & (8)\hfill \end{array}$

where dW⁽¹⁾_t, dW⁽²⁾_t are mutually independent, standard Wiener processes.

Second, consider a single buffer system with calcium diffusion, i.e., a system where calcium ions can diffuse into and out of the simulation volume. We will use a “constant pool” assumption, meaning that the simulation voxel is surrounded by a volume containing all molecular species at equilibrium concentrations. Thus, Ca²⁺ ions diffusing into and out of the simulation voxel change the simulations voxel's calcium concentration only, leaving the influx rate of calcium ions constant over time while the efflux rate varies. As diffusion is a monomolecular reaction, the deterministic and stochastic rate constants k_d and c_d are identical and given by $k_d=c_d=\frac{D}{L^2}$ (Elf et al., 2003; von Wegner and Fink, 2010). Here, D is the diffusion constant of calcium ions and the voxel volume is V = L³.

$\begin{array}{ll}\begin{array}{l}{\text{Ca}}^{2+}+\text{B}\stackrel{k^{+}}{\underset{k^{-}}{\rightleftharpoons }}\text{CaB}\\ {\text{Ca}}^{2+}\stackrel{k_d}{\underset{k_d}{\rightleftharpoons }}\varnothing \end{array}\hfill & (9)\hfill \end{array}$

The chemical Langevin equation for this system is:

$\begin{array}{ll}\begin{array}{l}d{X}_t=-c^{+}{X}_t{Y}_{t}dt+c^{-}{Z}_{t}dt-c_d{X}_{t}dt+c_d\tilde{X}dt\\ -\sqrt{c^{+}{X}_t{Y}_t}d{W}_t^{(1)}+\sqrt{c^{-}{Z}_t}d{W}_t^{(2)}-\sqrt{c_d{X}_t}d{W}_t^{(3)}\\ +\sqrt{c_d\tilde{X}}d{W}_t^{(4)}.\end{array}\hfill & (10)\hfill \end{array}$

Finally, we consider a multiple buffer system containing calcium ions and M different calcium buffers B_j, j = 1, …, M:

$\begin{array}{ll}\begin{array}{l}{\text{Ca}}^{2+}+{\text{B}}_1\stackrel{k_1^{+}}{\underset{k_1^{-}}{\rightleftharpoons }}{\text{CaB}}_1\\ ⋮\\ {\text{Ca}}^{2+}+{\text{B}}_M\stackrel{k_M^{+}}{\underset{k_M^{-}}{\rightleftharpoons }}{\text{CaB}}_M\end{array}\hfill & (11)\hfill \end{array}$

The corresponding chemical Langevin equation reads:

$\begin{array}{ll}\begin{array}{l}d{X}_t=-{\displaystyle \sum _{j=1}^M{c}_j^{+}}{X}_t{Y}_t^{(j)}dt+{\displaystyle \sum _{j=1}^M{c}_j^{-}}{Z}_t^{(j)}dt\\ -{\displaystyle \sum _{j=1}^M\sqrt{c_j^{+}{X}_t{Y}_t^{(j)}}}d{W}_t^{(j,1)}\text{}+{\displaystyle \sum _{j=1}^M\sqrt{c_j^{-}{Z}_t^{(j)}}}d{W}_t^{(j,2)}\end{array}\hfill & (12)\hfill \end{array}$

where any combination dW^(j,1)_t, dW^(j,2)_t is a pair of mutually independent Wiener processes.

2.3.3. The Ornstein-Uhlenbeck process (OUP)

The Ornstein-Uhlenbeck process uses a minimum of parameters to yield a stochastic process (X_t)_t≥0, X_t ∈ ℝ, t ∈ ℝ_≥0, which is stationary, Gaussian and exponentially autocorrelated in time. In one dimension, it is fully described by an initial condition X₀ and three parameters, the mean μ, the autocorrelation time τ, and the volatility σ. In the context of general stochastic processes, the term volatility is used to indicate a time-dependent parameter σ_t, leading to time-varying values of the variance of the process. Whereas σ is constant only for the more simple examples of stochastic processes, e.g., Brownian motion or the Ornstein-Uhlenbeck process, in general σ_t is itself a stochastic process. To take this feature into account from the beginning, we will use the term noise volatility. The variance of the Ornstein-Uhlenbeck process is constant and is given by $Var(X_t)=\frac{{\sigma }^2\tau }{2}$ (Gardiner, 2009). Both, buffered calcium dynamics and the OU-process are mean-reverting processes, i.e., deviations from the mean value μ drive the system back toward the mean. The autocorrelation time τ, or “noise color,” of the process quantifies how fast the system responds to deviations from the equilibrium value. Large τ values indicate strongly colored noise, i.e., a stochastic process with a slowly decaying exponential autocorrelation function. Small values of τ reflect short-lived fluctuations and the resulting noise approximates Gaussian white noise with τ approaching zero (Gillespie, 1996). The following stochastic differential equation describes the Ornstein-Uhlenbeck process:

$\begin{array}{ll}d{X}_t=-\frac{1}{\tau }(X_t-\mu )dt+\sigma d{W}_t.\hfill & (13)\hfill \end{array}$

In the following, we will try to recast the chemical Langevin equations presented above into forms similar to the Ornstein-Uhlenbeck process, where possible.

All simulations in the manuscript were coded and run in python 2.7.6 and can be obtained from the corresponding author by mail request. For larger simulation runs, we recommend the use of compiled code, e.g., Cython or C/C++.

3. Results

3.1. Approximate Stochastic Kinetics

In this section, the chemical Langevin equation is formulated for different systems, in particular for (i) single buffer systems without diffusion, (ii) single-buffer systems with calcium diffusion, and (iii) systems with several buffers but no diffusion. Application of the excess buffer approximation, Y_t = Ỹ leads to approximations of the Langevin equation. These approximations have a functional form similar to the Ornstein-Uhlenbeck process, however, the process parameters contain non-stationary terms. Substituting the equilibrium values of relevant molecular concentrations, non-stationarities can be eliminated and still, good approximations to exact simulation results are obtained. Furthermore, approximations of the CLE lead to estimates of the autocorrelation time τ and for the volatility σ of the process.

3.1.1. A single buffer system

For the single buffer system described by Equation 1, we substitute the EBA assumption Y_t = Ỹ and Z_t = X_T − X_t in Equation 8. Next, we make use of the fact that the sum of two independent, Gaussian random variables with means μ₁, μ₂ and variances σ²₁, σ²₂ yields another Gaussian random variable: (μ₁, σ²₁) + (μ₂, σ²₂) ~ (μ₁ + μ₂, σ²₁ + σ²₂) (Gillespie, 2000). Applying this formula to the mutually uncorrelated Brownian motion terms dW^(j)_t of the chemical Langevin equation, we obtain a single Brownian motion. The simplified Langevin equation of the single buffer system is then given by:

$d{X}_t=-(c^{+}\tilde{Y}+c^{-})X_{t}dt+c^{-}{X}_{T}dt+\sqrt{c^{+}{X}_t\tilde{Y}+c^{-}{Z}_t}d{W}_t.$

Using the definitions of the constants K₁ and K₂ (Equation 5), the analogous expressions in the stochastic framework are given by C₁ = c⁺Ỹ + c⁻ and C₂ = c⁻X_T. We obtain:

$\begin{array}{ll}d{X}_t=-C_1(X_t-\frac{C_2}{C_1})dt+\sqrt{c^{+}{X}_t\tilde{Y}+c^{-}{Z}_t}d{W}_t.\hfill & (14)\hfill \end{array}$

Now, using the fact that the equilibrium calcium concentration [Ca²⁺]_eq is related to $\tilde{X}$ by ${[{\text{Ca}}^{2+}]}_{eq}V=\tilde{X}=\frac{C_2}{C_1}$ and defining $\tau =\frac{1}{C_1}$ and the time-dependent volatility ${\sigma }_t=\sqrt{c^{+}{X}_t\tilde{Y}+c^{-}{Z}_t}$, we get:

$\begin{array}{ll}d{X}_t=-\frac{1}{\tau }(X_t-\tilde{X})dt+{\sigma }_{t}d{W}_t.\hfill & (15)\hfill \end{array}$

The last expression is almost identical to an Ornstein-Uhlenbeck process, Equation 13, except from the non-stationary term σ_t.

3.1.2. Single buffer with calcium diffusion

Consider a simple reaction system containing calcium ions, Ca²⁺, and a calcium buffer B. For the single buffer system with calcium diffusion, the CLE is transformed as:

$\begin{array}{ll}\begin{array}{l}d{X}_t=-c^{+}{X}_t\tilde{Y}dt+c^{-}(X_T-X_t)dt-c_d(X_t-\tilde{X})dt\\ -\sqrt{c^{+}{X}_t\tilde{Y}}d{W}_t^{(1)}+\sqrt{c^{-}{Z}_t}d{W}_t^{(2)}-\sqrt{c_d{X}_t}d{W}_t^{(3)}\\ +\sqrt{c_d\tilde{X}}d{W}_t^{(4)}\\ =(-C_1{X}_t+C_2)dt-c_d(X_t-\tilde{X})dt\\ +\sqrt{c^{+}{X}_t\tilde{Y}+c^{-}{Z}_t+c_d{X}_t+c_d\tilde{X}}d{W}_t\end{array}\hfill & (16)\hfill \end{array}$

The term −c_d (X_t − $\tilde{X}$) represents diffusion into and out of the simulation volume under the constant pool assumption. The term −c_d X_t is a function of the non-stationary calcium concentration X_t and represents the diffusive flux out of the simulation volume. Assuming an invariant equilibrium calcium concentration (a constant pool) outside of the simulation volume, the diffusive flux into the simulation volume, c_d $\tilde{X}$, is constant. Using C_1,2 as defined above, we get:

$\begin{array}{ll}\begin{array}{l}d{X}_t=-(C_1+c_d)(X_t-\tilde{X})dt\\ +\sqrt{c^{+}{X}_t\tilde{Y}+c^{-}{Z}_t+c_d(X_t+\tilde{X})}d{W}_t.\end{array}\hfill & (17)\hfill \end{array}$

3.1.3. Multiple buffer systems

Consider a reaction system containing calcium ions and M different calcium buffers B_j, j = 1, …, M. For multiple buffer systems, the following substitutions are applied to the chemical Langevin equation (Equation 12):

$\begin{array}{ll}\begin{array}{l}{Y}_t^{(j)}={\tilde{Y}}^{(j)}\\ Z_t^{(j)}=X_T-X_t-{\displaystyle \sum _{\begin{array}{c}i=1\\ i\ne j\end{array}}^M{Z}_t^{(i)}}.\end{array}\hfill & (18)\hfill \end{array}$

This last substitution, Equation 18, may seem to be of little help as we substitute the single, time-dependent variable Z^(j)_t by a more complicated expression involving the time-dependent variables X_t as well as all other CaB concentrations Z⁽ⁱ⁾_t, i ≠ j; however, this representation leads to C_1,2 terms that can easily be interpreted as modifications to the single buffer system. First, the corresponding CLE in the excess buffer approximation regime reads:

$\begin{array}{ll}\begin{array}{l}d{X}_t=-\left({\displaystyle \sum _{j=1}^M{c}_j^{+}}{\tilde{Y}}^{(j)}\right)X_{t}dt+{\displaystyle \sum _{j=1}^M{c}_j^{-}}\left(X_T-X_t-{\displaystyle \sum _{\begin{array}{c}i=1\\ i\ne j\end{array}}^M{Z}_t^{(i)}}\right)dt\\ +\sqrt{{\displaystyle \sum _{j=1}^M{c}_j^{+}}{X}_t{\tilde{Y}}^{(j)}+c_j^{-}{Z}_t^{(j)}}d{W}_t^{(j)}.\end{array}\hfill & (19)\hfill \end{array}$

To simplify, we will introduce a more compact notation. We will reuse the variable names C_1,2—however, it is important to note that in the context of multiple buffers, C₂ is not a constant term. This represents another source of non-stationarity. Defining

$\begin{array}{l}{C}_1={\displaystyle \sum _{j=1}^M{c}_j^{+}}{\tilde{Y}}^{(j)}+c_j^{-}\\ C_2={\displaystyle \sum _{j=1}^M{c}_j^{-}}\left(X_T-{\displaystyle \sum _{\begin{array}{c}i=1\\ i\ne j\end{array}}^M{Z}_t^{(i)}}\right)\end{array}$

as well as $\tau =\frac{1}{C_1}$ and ${\mu }_t=\frac{C_2}{C_1}$, we can rewrite:

$\begin{array}{ll}d{X}_t=-\frac{1}{\tau }(X_t-{\mu }_t)dt+\sqrt{{\displaystyle \sum _{j=1}^M{c}_j^{+}}{X}_t{\tilde{Y}}^{(j)}+c_j^{-}{Z}_t^{(j)}}d{W}_t.\hfill & (20)\hfill \end{array}$

Again, a functional form similar to an Ornstein-Uhlenbeck process, but with remaining non-stationary terms is obtained.

3.1.4. Estimated autocorrelation times τ and noise volatilities σ

In order to approximate the solutions of Gillespie's exact simulation algorithm and the chemical Langevin equation by a simple Ornstein-Uhlenbeck process, the non-stationarity contained in the terms μ_t and σ_t has to be eliminated by approximation. We observe that non-stationarity is introduced by the time-dependent reactant concentrations, X_t and Z_t. In the following, we substitute these expressions with their equilibrium values $\tilde{X}$ and $\tilde{Z}$ to obtain time-stationary expressions.

In the case of a single buffer without diffusion, the mean value is $\mu =\frac{C_2}{C_1}$, and the estimated autocorrelation time in the excess buffer approximation is given by:

$\begin{array}{ll}\begin{array}{l}\tau ={(c^{+}\tilde{Y}+c^{-})}^{-1}\\ =\frac{1}{C_1}.\end{array}\hfill & (21)\hfill \end{array}$

In terms of the deterministic rate constants, the total buffer concentration [B]_T and the equilibrium calcium concentration $\left[\widetilde{{\text{Ca}}^{2+}}\right]$, a more explicit expression for τ is:

$\begin{array}{ll}\tau =\frac{1}{k^{+}{[\text{B}]}_T\left(1-{(1+\frac{K_D}{[\widetilde{{\text{Ca}}^{2+}}]})}^{-1}\right)+k^{-}}.\hfill & (22)\hfill \end{array}$

Using the same approach to eliminate the non-stationary expression σ_t in Equation 15 using equilibrium concentrations, we obtain the constant term:

$\begin{array}{ll}\sigma =\sqrt{c^{+}\tilde{X}\tilde{Y}+c^{-}\tilde{Z}}.\hfill & (23)\hfill \end{array}$

The resulting variance of the calcium noise in the Ornstein-Uhlenbeck representation is then given by $Var(X)={\sigma }^2\frac{\tau }{2}$.

In the case of a single buffer with calcium diffusion, we obtained stationary expressions for $\mu =\frac{C_2}{C_1}$ and $\tau =\frac{1}{C_1}$, the mean calcium concentration and the autocorrelation time, respectively. When substituting equilibrium reactant concentrations in the non-stationary expression for σ_t, Equation 17, we get the stationary value σ for the approximating Ornstein-Uhlenbeck process:

$\sigma =\sqrt{\left(c^{+}\tilde{Y}+2c_d\right)\tilde{X}+c^{-}\tilde{Z}}.$

In the case of multiple buffers without diffusion, both parameters, μ_t and σ_t contain non-stationary terms as shown in Equation 20. Applying the equilibrium approximation, we get:

$\mu =\frac{{\displaystyle {\sum }_{j=1}^M{c}_j^{-}}\left(X_T-{\displaystyle {\sum }_{\begin{array}{c}i=1\\ i\ne j\end{array}}^M{\tilde{Z}}^{(i)}}\right)}{{\displaystyle {\sum }_{j=1}^M{c}_j^{+}{\tilde{Y}}^{(j)}+c_j^{-}}}.$

and

$\sigma =\sqrt{{\displaystyle \sum _{j=1}^M{c}_j^{+}}\tilde{X}{\tilde{Y}}^{(j)}+c_j^{-}{\tilde{Z}}^{(j)}}.$

Finally, the estimated autocorrelation time is

$\tau =\frac{1}{{\displaystyle {\sum }_{j=1}^M{c}_j^{+}{\tilde{Y}}^{(j)}+c_j^{-}}}.$

The results for the single buffer system are summarized in Figure 1, where the expected noise autocorrelation time τ and noise volatility σ are plotted as functions of different parameters and a simulation volume of V = 0.125 fl. Figure 1A shows the dependency of τ on the buffer association rate k⁺, at constant [Ca²⁺]_eq = 0.1 μM. This relationship is computed for different total buffer concentrations [B]_T (see legend). The shape of the curves can be derived from Equation 22, showing that the autocorrelation time τ decreases with faster buffer kinetics (k⁺) and with increasing [B]_T. Figure 1B illustrates the dependency of τ on the equilibrium free calcium concentration [Ca²⁺]_eq and it is observed that τ is larger at elevated calcium levels. The total buffer concentration was constant at [B]_T = 10 μM, while the dissociation constant was set to K_D = 1 μM. Each curve represents a different calcium binding rate k⁺ (see legend). Another section of the parameter space is explored in Figure 1C, with [Ca²⁺]_eq as the independent variable again. The total buffer concentration is constant at [B]_T = 100 μM and the buffer association rate is held constant at k⁺ = 0.1 μM⁻¹ ms⁻¹. This time, each curve corresponds to a different buffer dissociation constant K_D (see legend). The lower set of figures (Figures 1D–F) are analogous in design to Figures 1A–C. Here, the same dependencies are shown with the noise volatility σ as the dependent variable. It is noted that even in the range of physiological buffer kinetics, the noise parameters τ and σ differ substantially.

FIGURE 1

Figure 1. Theoretical noise autocorrelation time τ and noise volatility σ for the single buffer system without diffusion. (A) Shows the dependency of the noise autocorrelation time τ on the calcium binding rate constant k⁺ and for constant [Ca²⁺]_eq = 0.1 μM. Each curve corresponds to a fixed total buffer concentration [B]_T as indicated in the legend. (B) Shows that τ increases with increasing equilibrium calcium concentration [Ca²⁺]_eq, and for a fixed total buffer concentration [B]_T = 10 μM and fixed dissociation constant K_D = 1 μM. Each curve corresponds to a different calcium binding rate constant k⁺ (see legend). (C) Shows that τ increases with increasing [Ca²⁺]_eq, and for a fixed total buffer concentration [B]_T = 100 μM and fixed k⁺ = 0.1 μM⁻¹ ms⁻¹. Each curve corresponds to a different buffer dissociation constant K_D (see legend). (D–F) In analogy to (A–C) these lower row of panels shows the corresponding dependencies of the noise volatility σ. Note that all graphs are shown in log-log coordinates.

3.2. Numerical Results

3.2.1. Single buffer simulations

In this section, numerical results for a realistic single buffer system without diffusion are presented. In Figure 2, we chose the calcium-calmodulin system Ca²⁺ + CaM ⇌ Ca²⁺ · CaM in a microdomain of V = 0.125 fl and fixed the equilibrium calcium concentration [Ca²⁺] at 0.1 μM for all simulations. On the left side of Figures 2A–D, the intrinsic fluctuations of the local calcium concentration are shown for different simulation methods. From top to bottom, we show the results for Gillespie's exact stochastic simulation algorithm (SSA, black), the chemical Langevin approach (CLE, red), the CLE in the excess buffer approximation (CLE-EBA, green) and the Ornstein-Uhlenbeck approach (OUP, blue). To characterize these sample paths statistically, we show three basic statistical parameters of these processes, i.e., the mean value μ, the process' standard deviation σ and the autocorrelation time τ. Observe that these parameters show a high similarity between the simulation methods tested. Figure 2E shows the relaxational response of the system to a perturbation, where the initial calcium concentration was set to [Ca²⁺]_{t = 0} = 1 μM. Again, the response is tested for different methods as indicated by color (see Figures 2A–D). Additionally, the deterministic decay curve as calculated from analytical considerations, Equation 7, is shown (dashed black curves in E). Visual inspection shows a similar shape of the decay curves for the methods tested. Furthermore, it is seen that exact and approximated stochastic results fluctuate around the analytic solution of the deterministic system. In both cases, the results of Gillespie's exact stochastic simulation algorithm (black curves) should be taken as the reference curve for stochastic systems, as the associated algorithm samples the underlying process exactly. At this point, we conclude that the stochastic approximations to Gillespie's exact algorithm do not seem to affect the dynamic and stochastic properties of the single-buffer system significantly. A more detailed and quantitative assessment of this observation is presented in the following paragraphs. The relationship between equilibrium fluctuations and the relaxational response of the system is further analyzed quantitatively in Figure 3. For the frequently occurring calcium buffers given in Table 1, we compare the noise autocorrelation time τ of equilibrium calcium fluctuations, as computed from Equation 22, and the time constant $\widehat{\tau }$ that describes the exponential relaxation from an initial perturbation at t = 0. The system under consideration has a volume of V = 1.0 fl. As seen from Equations 22 and Equation 14, both values should be identical. In a realistic scenario, the perturbation could be due to the opening of a calcium channel. The relaxation time constant $\widehat{\tau }$ was estimated from n = 10 samples for each simulation method (see legend), and with an initial value of [Ca²⁺]_{t = 0} = 1.0 μM. It is observed that the values estimated from all simulations match the expected identity τ = $\widehat{\tau }$ (black line) with good accuracy, underscoring the usefulness of the stochastic approximations. Finally, in Figure 4 the expected values of τ and σ for several naturally occurring and artificial buffer species are shown. Note that all values refer to single buffer systems (V = 1.0 fl). Figure 4A shows the autocorrelation time τ as a function of the buffer association rate k⁺, and as a function of the buffer dissociation constant K_D in B, respectively. The lower row of panels (C, D), σ is shown for the same range of parameters. It is observed that τ and σ cover a range of several orders of magnitude. Numerical simulations (n = 10) using the SSA (red), CLE (green), OUP (blue) are shown as colored dots around the theoretical values (large black dots). Simulation results from different algorithms again show a high degree of accuracy and similarity among each other.

FIGURE 2

Figure 2. Equilibrium fluctuations and relaxation from transients. All simulations represent the calcium-calmodulin single buffer system in a microdomain of V = 0.125 fl. The panels on the left side (A–D) show equilibrium fluctuations of the free calcium concentration [Ca²⁺] as obtained from different methods. (A) Gillespie's exact stochastic simulation algorithm (SSA, black), (B) chemical Langevin equation (CLE, red), (C) the CLE in the excess buffer approximation (CLE-EBA, green), (D) the Ornstein-Uhlenbeck process approximation (OUP, blue). The values of the basic statistical properties analyzed in this paper are given to the right of each panel (mean μ, standard deviation σ, estimated autocorrelation time τ). The data suggests that all methods lead to similar statistical properties. The right panel (E) shows the relaxational response of the system to an initial perturbation, [Ca²⁺]_{t = 0} = 1 μM and colors indicate the corresponding simulation method (see A–D). The analytical solution of the corresponding deterministic system is also shown (dashed black lines). Visually, the decay curves for different simulation methods are nearly identical. A detailed and quantitative analysis of autocorrelation times and decay time constants is given in the main text.

FIGURE 3

Figure 3. Fluctuation-dissipation relations in a single buffer system without diffusion: for the set of calcium buffers from Table 1, we compare the noise autocorrelation time τ of equilibrium calcium fluctuations (Equation 22) and the time constant $\widehat{\tau }$ that describes the exponential relaxation from a perturbation (Equation 7), e.g., due to the opening of a calcium channel. For each simulation method (SSA, red; CLE, green; OUP, blue), $\widehat{\tau }$ was estimated for n = 10 sample paths starting at [Ca²⁺]_{t = 0} = 1.0 μM. The proximity to the identity (black line) shows that τ ≈ $\widehat{\tau }$.

FIGURE 4

Figure 4. Noise autocorrelation times τ and volatilities σ for several physiological and artificial buffer species. The autocorrelation time τ is shown as a function of the calcium binding rate constant k⁺ (A), and as a function of the buffer dissociation constant K_D (B), respectively. In the lower row (C,D), the calcium noise volatility σ is shown for the same parameters. It is observed that τ and σ cover a range of several orders of magnitude. Numerical simulations (n = 10) using the SSA (red), CLE (green), OUP (blue) are shown as colored dots around the theoretical values (large black dots).

3.2.2. Multiple buffers and diffusion

In order to add complexity and simulate a more realistic situation, the effect of multiple buffers and calcium diffusion on the autocorrelation time τ is analyzed in Figure 5. In the given example, τ and σ are analyzed for (i) a single buffer system (CaM), (ii) two different systems containing two buffers each (CaM+EGTA, CaM+Fluo-4), and (iii) a single buffer(CaM) system with calcium diffusion. Comparing the kinetic constants of the different buffers as given in Table 1, we see that both, EGTA and Fluo-4, share a similar calcium affinity. However, the association and dissociation events to and from Fluo-4 occur much faster than for EGTA. The constants as found in the literature differ by an approximate factor of 100. Consequently, the autocorrelation time of the CaM+EGTA system is mainly determined by the fast CaM dynamics and is hardly changed by addition of the two slow EGTA reactions. The comparatively few binding and dissociation events between Ca²⁺ and EGTA that can occur in a given time interval accelerate the decorrelation of the calcium fluctuations only slightly. Therefore, the autocorrelation time of the combined buffer system is smaller, but almost identical to the single buffer system containing CaM only (τ_CaM = 0.516 ms, τ_CaM,EGTA = 0.483 ms). In the case of adding Fluo-4 however, we find major changes in the statistical properties of the resulting calcium signal. Table 1 shows that especially the dissociation constant of the Ca-Fluo-4 complex is significantly larger than that of the Ca-CaM complex. The faster pace of reactions leads to an earlier decorrelation of the signal as shown by the resulting autocorrelation times (τ_CaM,Fluo−4 = 0.044 ms). The magnitude of diffusion effects lies in-between the effects of EGTA and Fluo-4, in the present example. The reader is warned not to generalize the relative effect of buffers and diffusion too quickly. From our derivation, it is explicitly clear that the effects also depend on the absolute concentration of buffers and calcium, respectively. In our case, the calcium diffusion constant of 200 μm²/s at a mean calcium concentration of 100 nM leads to an autocorrelation time of τ_CaM,Diff. = 0.190. At the same time, the noise volatility σ is influenced in the opposite direction. Clearly, the number of possible parameter combinations grows rapidly with the number of parameters allowed. The current article does not aim to scan this parameter space exhaustively.

FIGURE 5

Figure 5. The effect of multiple buffers and calcium diffusion on the autocorrelation time τ. In the example, the calcium noise autocorrelation time τ and volatility σ are analyzed for different conditions: a single buffer system (CaM), two double buffer systems (CaM+EGTA, CaM+Fluo-4), and a single buffer (CaM) system with calcium diffusion. Is is observed that the presence of an additional buffer leads to a faster decorrelation of the fluctuating calcium concentration, i.e., smaller τ values. Analogously, including diffusion effects also decreases τ values.

3.2.3. A colored noise-induced transition

In order to illustrate a potential physiological role of autocorrelated calcium noise, we revisit a generic example for the study of colored noise effects on dynamical systems. In particular, we consider the one-dimensional dynamics of a single variable X_t moving in a bistable potential represented by a fourth-order polynomial $V(X)=-\frac{a}{2}{X}^2+\frac{b}{4}{X}^4$ (Haenggi and Jung, 1995). The dynamics is driven by a colored noise process ξ_t representing calcium noise, and integration is performed using a first-order Euler scheme with dt = 0.01 ms. The equation of motion is given by

$\begin{array}{ll}d{X}_t=-{\partial }_{X}V(X_t)dt+{\xi }_t.\hfill & (24)\hfill \end{array}$

The system can represent the generalized dynamics of an ion channel whose two functional states (open, closed) are represented by the local minima of the potential (Liebovitch and Czegledy, 1992). The system is an abstract representation of bistable dynamics and can also represent more involved biochemical systems whose common feature is bistability. Figure 6 shows the results of two simulation runs, both using the Ornstein-Uhlenbeck approximation of calcium noise as the “noisy drive” ξ_t of the target variable X_t. To assess the effects of noise autocorrelation time, we used two different driving signals, one with a small autocorrelation time of τ = 0.01 (Figures 6A,C), and another one using τ = 0.75 (Figures 6B,D). By comparison with our previous results, we see that both τ-values fall within the range of autocorrelation times obtained with natural calcium buffers. In both cases, the noisy drive ξ_t induces a motion of the variable X_t that shows bistable (or “on-off”) kinetics due to the double-well shape of the potential. Using the analogy of a calcium-regulated ion channel, the “on-off” kinetics of the variable X_t represents the continuous gating of the calcium sensitive ion channel and we assume that the closed-open transition is calcium-dependent. To show the influence of noise correlations on channel gating, we look at the empirical joint probability distribution P(ξ, X) of the noise variable ξ_t and the channel variable X_t. At τ = 0.01, the joint distribution is almost radially symmetric around the origin, indicating a lack of correlation between the noise ξ_t and the channel variable X_t. In other words, the joint distribution approximately factorizes as P(ξ, X) ~ P(ξ) × P(X). At τ = 0.75 however, the result is a heavily skewed joint distribution, showing a strong coupling of channel gating (X_t) to the correlated calcium noise (ξ_t). It could be argued that this effect only occurs for the Ornstein-Uhlenbeck approximation of calcium noise. We also tested other algorithms (Gillespie's algorithm, chemical Langevin equation) to simulate ξ_t and observed the same qualitative effects (data not shown). Please note that in the present example, the mean value of the noise was removed to comply with the symmetry of V around the origin. Furthermore, in the example shown only τ was varied, keeping σ constant. This leads to differing variances for the noise processes, however, when correcting for this, the described effects persist.

FIGURE 6

Figure 6. A colored noise-induced transition. (A,C) The continuous gating of a calcium-sensitive ion channel can be reduced to the motion of a single particle X_t in a double-well potential. X_t is the channel variable and ξ_t is the calcium noise variable driving the motion of X_t (see Equation 24). The shape of the potential then leads to characteristic time series with “on-off” kinetics as observed in (A) (τ = 0.01) and (C) (τ = 0.75). Here, we assume the closed-open transition of the channel to be calcium-dependent. Calcium noise ξ_t is implemented with the OUP method (other methods yield qualitatively identical results, not shown). (B,D) Correlations between the calcium noise ξ and the channel state variable X are analyzed via their empirical joint probability distribution P(ξ, X). For a short autocorrelation time of the noise, τ = 0.01, the joint distribution is almost symmetrical around the origin, indicating that the noise and the channel variables have a low correlation (B). At τ = 0.75 however, the joint distribution is clearly skewed. This indicates a strong positive correlation between both variables (D). Such correlations between microdomain calcium noise and the functional state of calcium-sensitive target molecules, controlled by the autocorrelation time of calcium noise, await experimental validation.

4. Discussion

In the present article, we analyze stochastic simulation strategies for intracellular calcium microdomains. In a previous article, we investigated the potential of Gillespie's exact stochastic simulation algorithm to simulate calcium microdomains containing buffers and calcium-regulated calcium channels (von Wegner and Fink, 2010). Here, we focus on the properties of calcium noise itself. Already in our previous work, we assessed the effects of including diffusion events in the simulation by analyzing the autocorrelation time and the power spectral density of calcium noise. However, while our previous investigations were focused on accuracy, therefore implementing Gillespie's exact algorithm, we here concentrate on approximations to the exact solution of the chemical master equation. In particular, we use the chemical Langevin equation (Gillespie, 2000) and apply the excess buffer approximation (Heinemann et al., 1986; Smith et al., 2001; Fall et al., 2002) to simplify our results. We find that close to the resting calcium concentration, the excess buffer approximation (EBA) is a valid assumption over the range of physiological buffer concentrations. Fluctuations around the equilibrium buffer concentrations are very small, yielding coefficients of variation smaller than 10⁻³ (Table 1). For more complex situations, including Ca²⁺ channels for instance, it has to be remembered that buffer saturation can occur close to the channel pore, thereby invalidating the EBA assumption (Smith, 1996). We are aware of only one very recent publication using the chemical Langevin equation to analyze stochastic fluctuations of the calcium concentration in a microdomain and the influence of calcium buffers (Weinberg and Smith, 2014). For whole cell models, the chemical Langevin equation has been used more often (Li et al., 2005; Manninen et al., 2006). Other stochastic approaches for microdomains are based on explicit molecular random walk simulations (Franks et al., 2002; Shahrezaei and Delaney, 2004; Keller et al., 2008) and on the Fokker-Planck representation of the chemical master equation. The latter has been used to simulate calcium microdomains as found in the dyadic cleft of cardiac myocytes (Winslow et al., 2006; Hake and Lines, 2008).

The main focus of the present work are the autocorrelation properties of calcium noise, i.e., the noise “color.” Examples from many scientific areas show that noise color alone can determine a systems behavior and can even eventuate transitions between distinct stable states (Haenggi and Jung, 1995). However, so far there is no experimental evidence showing a functional role of calcium noise properties within a microdomain reaction network. To obtain quantitative results, we formulate the chemical Langevin equations of simplified model systems containing calcium ions, buffers and diffusion events. Combining noise terms and using the excess buffer approximation, we obtain surrogate processes that can be characterized by few parameters (μ, τ, σ), similar to an Ornstein-Uhlenbeck process. In some expressions, non-stationary terms were further substituted by their expected values at chemical equilibrium and thus, stationary expressions were obtained. With these expressions, the influence of buffer kinetic parameters and the calcium diffusion constant were quantified and illustrated for a section of the physiologically plausible parameter space. To summarize, we are able to characterize microdomain calcium noise as a colored noise process, noise color being characterized by physico-chemical system properties. In terms of computational speed, the chemical Langevin equation and the approximations presented here represent a major improvement in performance. Using the Ornstein-Uhlenbeck type approximations, the implementation of calcium noise is elementary and integration of realistic calcium noise into other simulation frameworks reduces to adding one further equation. In the context of systems described by ordinary or stochastic differential equations, adding calcium noise does not change the simulation framework conceptually as would be necessary when using Gillespie's algorithm.

To discuss situations in which the statistical properties of local calcium noise may be of physiological relevance, we chose a well studied example of a noise-induced transition in a generic bistable system (Figure 6, Haenggi and Jung, 1995). Bistability is a dynamical feature frequently encountered in biological systems, especially in intracellular signaling networks. Prominent examples are given by the genetic toggle switch (Tian and Burrage, 2006), phosphorylation-dephosphorylation networks (Kholodenko, 2006) and abstract representations of ion channels (Liebovitch and Czegledy, 1992). Each simulation method to generate calcium noise, namely Gillespie's algorithm, the chemical Langevin equation, and the Ornstein-Uhlenbeck representation, yield results analogous to Figure 6. The simulations show a coupling of the bistable variable to the noise process, or in terms of a more biological interpretation, a possible coupling of a calcium-regulated ion channel to the local calcium noise surrounding the channel. Research on more specific examples including realistic calcium channels and complex buffer situations are currently under way. Previously published complex models of cellular signaling networks driven by calcium implemented calcium noise as a random flux with a fixed rate (Bhalla, 2004). As our results suggest a possible functional role of calcium noise properties on calcium-sensitive reaction networks, the calcium noise implementation should be chosen carefully. Our results show how system parameters can be used to calculate realistic noise parameters.

In contrast to a full stochastic description of the system (Gillespie's algorithm), the Langevin method approximates the number of reactions in a fixed time interval as an appropiately scaled, normally distributed random variable with identical mean and variance. The principal limitation of the method becomes apparent with small reaction volumes, where the choice of a normally distributed variable can lead to physically meaningless negative molecule counts. While the solution to this problem would be to use the stochastically and physically correct Poisson distribution, the Langevin approximation has deeper theoretical connections and provides formulae for spectral properties and the cross-correlation structure of multi-variable systems (Gillespie, 2000; Gardiner, 2009; Weinberg and Smith, 2014). However, the validity of the Langevin approach has to be validated for each new model. A yet unexplored, but possibly interesting alternative could arise from the stochastic Cox-Ingersoll-Ross process (Goeing-Jaeschke and Yor, 2003). The process is also mean-reverting and its volatility σ_t scales as $\sqrt{X_t}$, similar to Equation 15.

Finally, the numerical predictions presented should receive further attention by experimentalists. An interesting perspective in this context is the possibility of performing measurements of localized calcium fluctuations by means of a fluorescent calcium indicator dye and a stationary laser, namely a fluorescent correlation spectroscopy (FCS) setup. To the best of our knowledge, Ca²⁺-FCS measurements have not been performed yet but the numerical results indicate that such a setup could be feasible. Such measurements would allow to quantify the magnitude and correlation structure of calcium flucutations in different cellular compartments and would allow conclusions about the buffer composition. In summary, as the attention of researchers shifts toward the stochastic nature of calcium dynamics, novel insights on the regulation of microdomain signaling networks are expected.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Bhalla, U. S. (2004). Signaling in small subcellular volumes. i. stochastic and diffusion effects on individual pathways. Biophys. J. 87, 733–744. doi: 10.1529/biophysj.104.040469

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bhattacharya, S., Bunick, C. G., and Chazin, W. J. (2004). Target selectivity in EF-hand calcium binding proteins. Biochim. Biophys. Acta 1742, 69–79. doi: 10.1016/j.bbamcr.2004.09.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Cheng, H., and Lederer, W. J. (2008). Calcium sparks. Physiol. Rev. 88, 1491–1545. doi: 10.1152/physrev.00030.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Choi, T., Maurya, M. R., Tartakovsky, D. M., and Subramaniam, S. (2010). Stochastic hybrid modeling of intracellular calcium dynamics. J. Chem. Phys. 133, 165101. doi: 10.1063/1.3496996

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Elf, J., Doncic, A., and Ehrenberg, M. (2003). Mesoscopic reaction-diffusion in intracellular signaling. Proc. SPIE 5110, 114–125. doi: 10.1117/12.497009

CrossRef Full Text | Google Scholar

Fall, C. P., Marland, E. S., Wagner, J. M., and Tyson, J. J. (eds.). (2002). Computational Cell Biology. New York, NY: Springer.

Franks, K. M., Bartol, T. M. Jr., and Sejnowski, T. J. (2002). A monte carlo model reveals independent signaling at central glutamatergic synapses. Biophys. J. 83, 2333–2348. doi: 10.1016/S0006-3495(02)75248-X

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gardiner, C. (2009). Stochastic Methods: A Handbook for the Natural and Social Sciences (Springer Series in Synergetics). Heidelberg: Springer.

Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361. doi: 10.1021/j100540a008

CrossRef Full Text | Google Scholar

Gillespie, D. T. (1996). Exact numerical simulation of the Ornstein-Uhlenbeck process and its integral. Phys. Rev. E 54, 2084–2091. doi: 10.1103/PhysRevE.54.2084

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gillespie, D. T. (2000). The chemical Langevin equation. J. Chem. Phys. 113, 297–306. doi: 10.1063/1.481811

CrossRef Full Text | Google Scholar

Gillespie, D. T. (2007). Stochastic simulation of chemical kinetics. Ann. Rev. Phys. Chem. 58, 35–55. doi: 10.1146/annurev.physchem.58.032806.104637

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Goeing-Jaeschke, A., and Yor, M. (2003). A survey and some generalizations of bessel processes. Bernoulli 9, 313–349. doi: 10.3150/bj/1068128980

CrossRef Full Text | Google Scholar

Haenggi, P., and Jung, P. (1995). Colored noise in dynamical systems. Adv. Chem. Phys. LXXXIX, 239–326.

Hake, J., and Lines, G. T. (2008). Stochastic binding of ca2+ ions in the dyadic cleft; continuous versus random walk description of diffusion. Biophys. J. 94, 4184–4201. doi: 10.1529/biophysj.106.103523

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Heinemann, U., Klee, M., Neher, E., and Singer, W. (eds.). (1986). Calcium Electrogenesis and Neuronal Functioning. Experimental Brain Research Series, Vol. 14. Berlin: Springer. doi: 10.1007/978-3-642-70744-5

CrossRef Full Text

Higham, D. J. (2001). An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev. 43, 525–546. doi: 10.1137/S0036144500378302

CrossRef Full Text | Google Scholar

Jiang, Y. H., Klein, M. G., and Schneider, M. F. (1999). Numerical simulation of Ca²⁺ “sparks” in skeletal muscle. Biophys. J. 77, 2333–2357. doi: 10.1016/S0006-3495(99)77072-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Keller, D. X., Franks, K. M., Bartol, T. M. Jr., and Sejnowski, T. J. (2008). Calmodulin activation by calcium transients in the postsynaptic density of dendritic spines. PLoS ONE 3:e2045. doi: 10.1371/journal.pone.0002045

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kholodenko, B. N. (2006). Cell-signalling dynamics in time and space. Nat. Rev. Mol. Cell Biol. 7, 165–176. doi: 10.1038/nrm1838

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Li, H., Hou, Z., and Xin, H. (2005). Internal noise stochastic resonance for intracellular calcium oscillations in a cell system. Phys. Rev. E Stat. Nonlin. Soft Matt. Phys. 71(6 Pt 1):061916. doi: 10.1103/PhysRevE.71.061916

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Liebovitch, L. S., and Czegledy, F. P. (1992). A model of ion channel kinetics based on deterministic, chaotic motion in a potential with two local minima. Ann. Biomed. Eng. 20, 517–531. doi: 10.1007/BF02368171

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Manninen, T., Linne, M.-L., and Ruohonen, K. (2006). Developing itô stochastic differential equation models for neuronal signal transduction pathways. Comput. Biol. Chem. 30, 280–291. doi: 10.1016/j.compbiolchem.2006.04.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Novo, D., DiFranco, M., and Vergara, J. L. (2003). Comparison between the predictions of diffusion-reaction models and localized Ca²⁺ transients in amphibian skeletal muscle fibers. Biophys. J. 85, 1080–1097. doi: 10.1016/S0006-3495(03)74546-9

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Shahrezaei, V., and Delaney, K. R. (2004). Consequences of molecular-level ca2+ channel and synaptic vesicle colocalization for the ca2+ microdomain and neurotransmitter exocytosis: a monte carlo study. Biophys. J. 87, 2352–2364. doi: 10.1529/biophysj.104.043380

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Smith, G. D. (1996). Analytical steady-state solution to the rapid buffering approximation near an open ca2+ channel. Biophys. J. 71, 3064–3072. doi: 10.1016/S0006-3495(96)79500-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Smith, G. D., Dai, L., Miura, R. M., and Sherman, A. (2001). Asymptotic analysis of buffered calcium diffusion near a point source. SIAM J. Appl. Math. 61, 1816–1838. doi: 10.1137/S0036139900368996

CrossRef Full Text | Google Scholar

Smith, G. D., Keizer, J. E., Stern, M. D., Lederer, W. J., and Cheng, H. (1998). A simple numerical model of calcium spark formation and detection in cardiac myocytes. Biophys. J. 75, 15–32. doi: 10.1016/S0006-3495(98)77491-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Tian, T., and Burrage, K. (2006). Stochastic models for regulatory networks of the genetic toggle switch. Proc. Natl. Acad. Sci. U.S.A. 103, 8372–8377. doi: 10.1073/pnas.0507818103

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Uttenweiler, D., Kirsch, W. G., Schulzke, E., Both, M., and Fink, R. H. A. (2002). Model-based analysis of elementary Ca²⁺ release events in skinned mammalian skeletal muscle fibres. Eur. Biophys. J. 31, 331–340. doi: 10.1007/s00249-002-0223-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

von Wegner, F., and Fink, R. H. A. (2010). Stochastic simulation of calcium microdomains in the vicinity of an L-type calcium channel. Eur. Biophys. J. 39, 1079–1088. doi: 10.1007/s00249-009-0504-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

von Wegner, F., Wieder, N., and Fink, R. H. A. (2012). Simulation strategies for calcium microdomains and calcium-regulated calcium channels. Adv. Exp. Med. Biol. 740, 553–567. doi: 10.1007/978-94-007-2888-2-25

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Weinberg, S. H., and Smith, G. D. (2014). The influence of Ca(2+) buffers on free [Ca(2+)] fluctuations and the effective volume of Ca(2+) microdomains. Biophys. J. 106, 2693–2709. doi: 10.1016/j.bpj.2014.04.045

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wieder, N., Fink, R. H. A., and Wegner, F. v. (2011). Exact and approximate stochastic simulation of intracellular calcium dynamics. J. Biomed. Biotechnol. 2011, 572492. doi: 10.1155/2011/572492

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Winslow, R. L., Tanskanen, A., Chen, M., and Greenstein, J. L. (2006). Multiscale modeling of calcium signaling in the cardiac dyad. Ann. N.Y. Acad. Sci. 1080, 362–375. doi: 10.1196/annals.1380.027

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zhang, J., Hou, Z., and Xin, H. (2004). System-size biresonance for intracellular calcium signaling. Chemphyschem 5, 1041–1045. doi: 10.1002/cphc.200400089

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zhu, C.-L., Jia, Y., Liu, Q., Yang, L.-J., and Zhan, X. (2007). A mesoscopic stochastic mechanism of cytosolic calcium oscillations. Biophys. Chem. 125, 201–212. doi: 10.1016/j.bpc.2006.08.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar