LETTER

Communicated by Nicolas Wieder

Microdomain [Ca2+ ] Fluctuations Alter Temporal Dynamics in Models of Ca2+ -Dependent Signaling Cascades and Synaptic Vesicle Release Seth H. Weinberg Virginia Modeling, Analysis and Simulation Center, Old Dominion University, Suffolk, Virginia 23435, U.S.A.

Ca2+ -dependent signaling is often localized in spatially restricted microdomains and may involve only 1 to 100 Ca2+ ions. Fluctuations in the microdomain Ca2+ concentration ([Ca2+ ]) can arise from a wide range of elementary processes, including diffusion, Ca2+ influx, and association/dissociation with Ca2+ binding proteins or buffers. However, it is unclear to what extent these fluctuations alter Ca2+ -dependent signaling. We construct Markov models of a general Ca2+ -dependent signaling cascade and Ca2+ -triggered synaptic vesicle release. We compare the hitting (release) time distribution and statistics for models that account for [Ca2+ ] fluctuations with the corresponding models that neglect these fluctuations. In general, when [Ca2+ ] fluctuations are much faster than the characteristic time for the signaling event, the hitting time distributions and statistics for the models with and without [Ca2+ ] fluctuation are similar. However, when the timescale of [Ca2+ ] fluctuations is on the same order as the signaling cascade or slower, the hitting time mean and variability are typically increased, in particular when the average number of microdomain Ca2+ ions is small, a consequence of a long-tailed hitting time distribution. In a model of Ca2+ -triggered synaptic vesicle release, we demonstrate the conditions for which [Ca2+ ] fluctuations do and do not alter the distribution, mean, and variability of release timing. We find that both the release time mean and variability can be increased, demonstrating that [Ca2+ ] fluctuations are an important aspect of microdomain Ca2+ signaling and further suggesting that [Ca2+ ] fluctuations in the presynaptic terminal may contribute to variability in synaptic vesicle release and thus variability in neuronal spiking. 1 Introduction In many cell types, Ca2+ -dependent signaling regulates important physiologically subcellular processes. In neurons, Ca2+ signaling often occurs in spatially localized subspaces, or Ca2+ microdomains, due to spatially confined or restricted Ca2+ diffusion (Korkotian & Segal, 2006; Biess, Korkotian, & Holcman, 2011) and triggers many key physiological functions, including Neural Computation 28, 493–524 (2016) doi:10.1162/NECO_a_00811

c 2016 Massachusetts Institute of Technology 

494

S. Weinberg

synaptic vesicle release (SVR) (Berridge, 1998; Berridge, Bootman, & Roderick, 2003; Bollmann & Sakmann, 2005). Furthermore, Ca2+ signaling often occurs over a wide range of timescales, on the millisecond timescale for SVR (Smith, 2004), to many hours, for example, for Ca2+ -calmodulin-dependent protein kinase gene transcription regulation (Soderling, 1999; Carrasco & Hidalgo, 2006). In a biochemical reaction network, it is well known that the number of molecules in the system will randomly fluctuate due to the random timing of individual reactions (Keizer, 1987). Fluctuations in molecular concentration are large in amplitude when the system size is small, while for sufficiently large system size, concentration fluctuations become negligible (Keizer, 1987). For bimolecular reactions, it has been shown that the expected value of the stationary distribution for the stochastic system does not agree with the corresponding deterministic steady state (McQuarrie, Jachimowski, & Russell, 1964; Darvey, Ninham, & Staff, 1966). Morita (1988) studied concentration fluctuations in the reversible bimolecular association reaction, in which the fluctuating species was externally driven by gaussian white noise and had sufficiently large concentration, such that the likelihood of negative concentration was considered negligible. He showed that concentration fluctuations also alter reaction dynamics, in addition to the steady state, compared with the deterministic model. In physiological conditions associated with Ca2+ microdomains, the relevant system size is often small, such that the likelihood of zero Ca2+ ions is nonnegligible and significant: intracellular [Ca2+ ] in cells at rest is typically very low, near 100 nM, and microdomain volumes are between 0.001 and 1 μm3 (Harris, Jensen, & Tsao, 1992), which correspond to 0.06 to 60 Ca2+ ions. Many previous studies have focused on the influence of stochasticity in Ca2+ -driven processes, such as gating of Ca2+ -activated channels (Williams et al., 2011; Gaur & Rudy, 2011; Cannell, Kong, Imtiaz, & Laver, 2013). Recently, several modeling studies have demonstrated that [Ca2+ ] fluctuations may significantly influence the properties of Ca2+ -dependent processes, specifically Ca2+ -regulated Ca2+ channels (Wieder, Fink, & von Wegner, ¨ 2011, 2015). Flegg, Rudiger, and Erban (2013) found that [Ca2+ ] fluctuations reduce the interpuff time in a cluster of Ca2+ -activated Ca2+ channels. Tanskanen, Greenstein, Chen, Sun, and Winslow (2007) showed that a continuous model of [Ca2+ ] overestimates channel open probability in cardiac myocytes, compared with a model that accounts for individual Ca2+ ions and protein geometry. Weinberg and Smith (2012) similarly showed in a minimal Ca2+ microdomain model that [Ca2+ ] fluctuations typically reduced steady-state channel open probability in Ca2+ -activated channels. Hake and Lines (2008) argued that a random-walk model of Ca2+ signaling in the cardiac dyad can be well approximated by a deterministic representation of Ca2+ diffusion, coupled with stochastic channel gating, and that large fluctuations due to diffusion are smoothed at the longer timescale

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

495

of channel gating. In this study, we sought to determine if smoothing of fast [Ca2+ ] fluctuations over a longer timescale is a general feature of Ca2+ dependent signaling. More specifically, in neurons, do presynaptic terminal [Ca2+ ] fluctuations alter release timing in a model of Ca2+ -triggered SVR? If so, at what timescales does this occur, and how does this depend on the properties of the Ca2+ -dependent signaling process? During peak Ca2+ entry in the presynaptic terminal, SVR probability is expected to be maximal, while when free [Ca2+ ] is low and highly buffered, SVR probability is expected to be low. However, Weinberg and Smith (2014) recently showed that Ca2+ binding proteins or buffers, which may be expected to reduce [Ca2+ ] fluctuations, in fact did not reduce and may enhance the magnitude of [Ca2+ ] fluctuations. Further, this enhancement of [Ca2+ ] fluctuations increased as average [Ca2+ ] increased. Thus, it is not clear to what extent [Ca2+ ] fluctuations may influence SVR timing and how this may depend on Ca2+ microdomain and SVR properties. In our prior work, we showed that the relative size (coefficient of variation) of [Ca2+ ] fluctuations, driven solely by passive exchange (diffusion) (Schuss, Singer, & Holcman, 2007), is inversely proportional with the square √ root of the average number of microdomain Ca2+ ions—that is, 1/ cV, where c is the average [Ca2+ ] and V is the microdomain volume—and, as noted above, Ca2+ buffers may in fact enhance the size of fluctuations, as the Ca2+ -buffer association and dissociation reactions are an additional source of intrinsic noise (Weinberg & Smith, 2012, 2014). Further, the timescale of microdomain [Ca2+ ] fluctuations, specifically the autocorrelation time, can range over several orders of magnitude, from approximately 10−2 to 102 ms, depending on the properties of buffers (e.g., concentration and binding affinity, and microdomain geometry; Weinberg & Smith, 2014; von Wegner, Wieder, & Fink, 2014). In this study, we investigate how the average [Ca2+ ], microdomain volume, and [Ca2+ ] fluctuations timescale influence the dynamics of Ca2+ -dependent signaling, in particular, multistep signaling cascades. Do [Ca2+ ] fluctuations alter Ca2+ -dependent signaling when fluctuations are much faster, slower, or the same timescale as the Ca2+ dependent signaling event? Answering this question is complicated by the fact that the Ca2+ -dependent signaling event itself is a source of intrinsic [Ca2+ ] fluctuations. Since [Ca2+ ] fluctuations are of significant amplitude due to the small system size, we ask a related (and, in some sense, equivalent) question: Does a model that neglects [Ca2+ ] fluctuations under- or overestimate the temporal dynamics of Ca2+ -dependent signaling? Previously, Holcman and Schuss (2005) described the two-dimensional Markov chain to investigate concentration fluctuations driven by biochemical reactions and diffusion. Dao Duc and Holcman (2010) used a similar two-dimensional Markov chain description to characterize the mean first passage time for a molecule to sequentially bind an immobile target site, studying the situation in which the mobile effector molecule (comparable

496

S. Weinberg

to Ca2+ in our study) can and cannot escape the microdomain (i.e., passive exchange out of the domain). Our study extends this analysis by also considering passive exchange into the domain from the bulk, an important source of intrinsic noise. In this study, we first investigate a general signaling cascade and how [Ca2+ ] fluctuations timescale and microdomain volume may influence the timing of Ca2+ -dependent signaling. We derive analytical approximations for the hitting time statistics in the limit of a small microdomain volume. Neher and others demonstrated that SVR can be modeled as a multistep signaling cascade, the sequential Ca2+ binding steps preceding a final Ca2+ independent step (Heidelberger, Heinemann, Neher, & Matthews, 1994; Schneggenburger & Neher, 2000; Neher & Sakaba, 2008), and following the approach applied to the general signaling cascade, we analyze the influence of [Ca2+ ] fluctuations on Ca2+ -triggered SVR. Finally, we simulate SVR triggered by neuron spiking and characterize the variability in SVR timing. 2 Methods 2.1 Calculation of Hitting Times in Signal Cascades. In any biochemical signaling cascade, the time for a transition to a downstream state, a consequence of a biochemical reaction, is inherently random. The final stage of a signaling cascade often represents a physiologically important event, such as triggering neurotransmitter release (Bollmann & Sakmann, 2005) or complete phosphorylation of a protein with multiple sites (Yang, MacLellan, Han, Weiss, & Qu, 2004). We can represent a Ca2+ -dependent signaling cascade, with or without [Ca2+ ] fluctuations present, that is, a model that either accounts for or neglects [Ca2+ ] fluctuations, by a discretestate continuous-time Markov chain (CTMC). If we define M(t) to denote the state of the CTMC that takes integer values between 1 and N + 1 (the total number of states), then the probability of a transition from state i to state j during a small time interval dt is given by qi j dt = Pr{M(t + dt) = j|M(t) = i}, (i = j),

(2.1)

where matrix Q = (qi j ) is known as the infinitesimal generator matrix. Conservation  of probability implies that the diagonal entries of Q are given by qii = − j=i qi j . For a CTMC with a single absorbing state, N + 1, Q can be written in the form   T u Q= , (2.2) 0 0 where T is an N × N matrix of transition rates between transient states, u = −Te collects the rates into the absorbing state, and e and 0 are a column

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

497

vector of ones and a row vector of zeros, respectively, of length N (Breuer & Baum, 2005; Bladt, 2005; Lamar, Kemper, & Smith, 2011). The hitting time for state N + 1, or time until absorption of an absorbing CTMC, τ is well known to have a phase-type distribution (Breuer & Baum, 2005; Bladt, 2005; Lamar et al., 2011), written as τ ∼ PH(ζ, T ), where ζ = (ζi ) is a row vector denoting the initial probability of state i such that i ζi = 1. The probability density function for τ , f (x)dx = Pr{x < τ ≤ x + dx}, that is, the hitting time distribution, is given by Breuer and Baum (2005), Bladt (2005), and Lamar et al. (2011), f (x) = ζ exp(xT )uH (x),

(2.3)

where exp(xT ) is a matrix exponential and H (x) is the standard Heaviside function. The expected (mean) hitting time τ is given by  μ=

∞

x f (x)dx = −ζ T−1 e,

(2.4)

0

and in general, the qth moment of f (x) is given by  E[τ q ] =

∞

xq f (x)dx = (−1)q q!ζ T−q e,

(2.5)

0

where T−q is a shorthand for the matrix inverse, raised to the qth power. The hitting time τ variability is characterized by the coefficient of variation (CV), the hitting time standard deviation divided by the mean, cv = σ /μ, where the hitting time variance σ 2 = E[τ 2 ] − μ2 . In the following analysis, a Ca2+ -dependent signaling cascade, either a generic cascade or specifically a model of SVR, is represented as an absorbing CTMC, with the final cascade stage representing the absorbing state. Our analysis focuses on calculating the hitting time for this absorbing state, assuming an initial state that corresponds to the first cascade stage and a specified [Ca2+ ]. In general, analytical expressions for the hitting time probability density function and moments can be calculated using equations 2.3 to 2.5. However, these expressions involve calculation of the matrix exponential and inverse of T (often a very large matrix); as such, these expressions are complex and in most cases do not provide insight. Therefore, numerical calculations are performed and presented. In the limit of a small microdomain volume, we will derive expressions for the hitting time mean and variance. 2.1.1 In the Absence of Microdomain [Ca2+ ] Fluctuations. We present a general approach to analyze the influence of [Ca2+ ] fluctuations on a multistep Ca2+ -dependent signaling cascades. The most general signaling cascade

498

S. Weinberg

consists of n reversible steps, with transition rates λ1 , . . . , λn , λ−1 , . . . , λ−n , with units of time−1 , given by λ1

λ2

λn−1

λ−1

λ−2

λ−(n−1)

S1  S2  · · ·



λn

Sn  Sn+1 . λ−n

(2.6)

Sn represents possible system states, which, for example, could represent states of a protein with multiple binding sites or states of an ion channel with multiple closed, inactivated, and open states. For any generic signaling cascade, the hitting time distribution, expected value, and variability can be calculated using equations 2.3 to 2.5 and provide a reference for comparison with numerical calculations of these measures when [Ca2+ ] fluctuations are present. 2.1.2 In the Presence of Microdomain [Ca2+ ] Fluctuations. As discussed in section 1, the timescale of [Ca2+ ] fluctuations depends on many factors, such as microdomain geometry and the presence and properties of Ca2+ buffers. Since the size and timescale of [Ca2+ ] fluctuations are complex functions of buffer properties and inherently linked via these dependencies (Weinberg & Smith, 2014; von Wegner et al., 2014), we decouple these characteristics in this study by varying the average [Ca2+ ], microdomain volume, and timescale separately. To investigate the influence of [Ca2+ ] fluctuations in general, without the introduction of several buffer- or other noise-related parameters, we consider [Ca2+ ] fluctuations that arise as a consequence of passive exchange (diffusion) between the microdomain and a bulk compartment, characterized by a single parameter τˆe . If Ca2+ passive exchange occurs with rate kˆ e , with units of time−1 , then the timescale of [Ca2+ ] fluctuations τˆe , specifically the autocorrelation time, is simply given by τˆe = 1/kˆ e (Weinberg & Smith, 2014). To investigate the influence of intrinsic [Ca2+ ] fluctuations on cascade signaling, we consider a system of reactions including passive exchange (diffusion) of Ca2+ (Schuss et al., 2007), ∅

kˆ e c∞  kˆ e

Ca2+ ,

(2.7a)

where c∞ is the average [Ca2+ ] in the microdomain and bulk, and further modify the reaction scheme in equation 2.6 to account for Ca2+ -dependence, which arises due to a Ca2+ -association reaction, that is, S1 + Ca2+

k1  k−1

S2

(2.7b)

S2 + Ca2+

k2  k−2

S3

(2.7c)

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

.. .

.. .

499

.. .

Sn + Ca2+

kn  k−n

Sn+1 ,

(2.7d)

where ki is an association rate constant, with units of concentration−1 · time−1 . For appropriate comparison with the hitting time calculation in the absence of [Ca2+ ] fluctuations, we define ki = λi /c∞ and k−i = λ−i . While Ca2+ -dependence may arise via a more complicated biochemical reaction scheme, we consider the bimolecular reaction that has the simplest linear dependence in the large system size limit, that is, when governed by mass-action kinetics, the rates of the reactions in equations 2.7b to 2.7d are proportional to [Ca2+ ]. If we define the state of the two-dimensional CTMC that accounts for [Ca2+ ] fluctuations and Ca2+ -dependent state transitions by M = (C, S), where C ∈ {0, 1, . . . , ∞} ranges over all possible discrete values for the number of microdomain Ca2+ ions, S ∈ {1, . . . , n + 1} represents the system state of a signal cascade, and n is the number of cascade steps, we can write the state-transition diagram for the CTMC model that is governed by equation 2.7 in a microdomain of volume V: (0, 1)

α1  β

(1, 1)

α1  2β

α1  3β

2k¯ 1 k−1

k¯ 1 k−1 (0, 2)

(2, 1)

α2  β

(1, 2)

(3, 1)

···

α2  3β

···

.. .

..

αn  3β

···

3k¯ 1 k−1 α2  2β

(2, 2)

k¯ 2 k−2 .. .

2k¯ 2 k−2 .. .

2k¯ n−1 k−(n−1)

3k¯ n−1 k−(n−1)

(1, n)

α1  4β

αn  2β

(2, n)

k¯ n k−n

2k¯ n k−n

(∗, n + 1)

(∗, n + 1)

.

,

··· (2.8)

where k¯ i = ki /V = λi /(c∞V ), k−i = λ−i , α1 = α2 = · · · = αn = kˆ e c∞V, and β = kˆ e . Note that the transition rate between cascade stages k¯ i is inversely

500

S. Weinberg

proportional to the microdomain volume V, while the influx rate due to passive exchange αi is proportional to V. In the transition diagram, transitions between states within a row arise from passive exchange with the bulk, equation 2.7a, and transitions between rows arise from a change from system state i to state i + 1, equations 2.7b to 2.7d. Each ordered pair M = (C, S) represents a state of two-dimensional CTMC. Since we are interested in the time until the system is absorbed into system state Sn+1 , we can lump the bottom row into a single CTMC state. Further, while in theory, the CTMC represented by equation 2.8 has an infinite state space, since C is unbounded, in practice, we can limit the range for C ∈ {0, 1, . . . , Cm }, where Cm is a maximum number of Ca2+ ions. This value can be determined by incrementally increasing Cm until there are insignificant changes in numerical calculations in hitting distribution and moments. We found that a reasonable value for Cm = max( 2c∞V, 50), where · is the ceiling function, yielding a Q matrix of dimension n(Cm + 1) + 1. After enumerating the N = n(Cm + 1) transient system states, such that state 1 corresponds to (0,1), 2 to (1,1), . . . , and the absorbing state N + 1 corresponds to (∗, n + 1), we can write the infinitesimal generator matrix Q corresponding to the finite transition diagram, equation 2.8, in the form given by equation 2.2. We define the initial probability vector ζ such that there is agreement between models with and without [Ca2+ ] fluctuations present (see equation 2.6 and kˆ e = 0 in equation 2.8) for an appropriate comparison. Since the number of Ca2+ ions must be an integer and in general c∞V is not an integer, the initial probability is distributed appropriately between states π = c∞V and π − 1, such that the initial expected number of Ca2+ ions is equal to c∞V. The timescales associated with [Ca2+ ] fluctuations and Ca2+ -dependent signaling cascades can range over several orders of magnitude. Here, we focus on the influence of the relative timescale between [Ca2+ ] fluctuations and signaling; thus, we define λi = 1 (in arbitrary units of time−1 ) such that the passive exchange rate kˆ e = λi /τe , where τe is a unitless [Ca2+ ] fluctuations timescale parameter (equivalently, we define τe = λi /kˆ e ). Importantly, while τe defines the [Ca2+ ] fluctuations timescale in a microdomain in which passive exchange is the only elementary process (Weinberg & Smith, 2014), here, since we explicitly account for Ca2+ binding in the Ca2+ -triggered signaling, the association and dissociation reactions are an additional source of intrinsic noise that also drive [Ca2+ ] fluctuations. In the subsequent analysis, hitting time distributions and moments are numerically calculated using equations 2.3 to 2.5. We calculate the hitting time distribution, equation 2.3, and mean μ, ˆ equation 2.4, normalized by the mean of the model without [Ca2+ ] fluctuations. We also calculate the (not normalized) hitting variability cv to illustrate how model parameters alter the hitting time variability and the normalized hitting time variability, cˆv , to illustrate to what extent [Ca2+ ] fluctuations alter hitting time variability.

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

501

It is worth noting that calculations of the hitting time distribution and moments can be determined by Monte Carlo simulations (Gillespie, 1977), shown in the online supporting material. However, the approach used in our study enables direct numerical calculations of these same measures without simulation. As discussed below, for many parameters the hitting time distribution has a long tail, such that rare events (relative to an exponential distribution) occur with an increased likelihood. Accurate estimation of the long-tailed hitting time distribution and moments requires a large number of simulations, between 103 and 104 (see Figure S1), such that analysis of the wide range of cascade parameters and timescales considered in this study would be prohibitive. We also note that our analysis assumes a well-mixed Ca2+ microdomain and does not consider spatial aspects of diffusion and Ca2+ binding. The validity of this compartment-based modeling assumption will depend on microdomain geometry and the passive exchange rate (see also section 4.2). Prior modeling studies have accounted for spatial aspects of Ca2+ signaling. ¨ Rudiger, Falcke, and colleagues used a hybrid stochastic spatial modeling approach, accounting for stochastic aspects of Ca2+ channel gating and specific channel locations but neglecting fluctuations arising from diffusion ¨ ¨ (Nagaiah, Rudiger, Warnecke, & Falcke, 2012; Ruckl et al., 2015). Flegg, ¨ Rudiger, and Erban (2013) recently compared this hybrid model with a fully stochastic model with Brownian Ca2+ ion motion and showed that neglecting [Ca2+ ] fluctuations altered channel opening time. Modchang et al. (2010) investigated SVR in a stochastic spatial model and found that fluctuations in Ca2+ channel gating were the most significant noise source influencing SVR probability, but also that Ca2+ diffusion was a significant noise source when peak [Ca2+ ] and release probability were low. Although our study does not fully account for spatial aspects of [Ca2+ ] fluctuations, our compartment-based formulation enables numerical calculations for hitting time distributions, rather than computationally intensive Monte Carlo simulations, and thus facilitates a wide range of parameter studies. Additionally, our numerical approach can be used to predict parameter regimes of interest that can be investigated in a fully stochastic spatial simulation (Flegg et al., 2012). Further, Nadkarni, Bartol, Sejnowski, and Levine (2010) showed in a stochastic spatial model that the timescales associated with SVR are independent of the detailed synaptic geometry, suggesting that a compartment-based model may be sufficient to reproduce stochastic SVR timing dynamics. 3 Results 3.1 General Signaling Cascade 3.1.1 Single Irreversible Reaction. We first consider the simplest possible general signaling cascade, a single irreversible reaction:

502

S. Weinberg

Figure 1: [Ca2+ ] fluctuations can induce long-tailed hitting time distributions. The hitting time distribution for the single-step signaling cascade (see equation 2.8) is shown for different values of the [Ca2+ ] fluctuations timescale parameter τe (solid color lines). The distribution in the absence of [Ca2+ ] fluctuations is shown for comparison (black dashed line). The distribution is shown on a linear (left) and logarithmic (right) scale. Parameters: Microdomain volume V = 0.1 (A) and 0.01 (B) μm3 , c∞ = 0.1 μM.

λ1

S1 → S2 .

(3.1)

In the absence of fluctuations, the hitting time τ , in this case, the time for formation of one S2 molecule assuming one initial S1 molecule, is well characterized: the distribution is exponentially distributed, τ ∼ Exp(λ1 ), −2 2 with an expected hitting time μ = λ−1 1 . The hitting time variance σ = λ1 , and thus the relative hitting time variability cv = 1. In Figure 1, we plot the hitting time distribution, equation 2.3, for the CTMC that includes [Ca2+ ] fluctuations with various timescales (solid line), compared with the model that does not include fluctuations (dashed black line). In a microdomain of volume V = 0.1 μm3 and a resting [Ca2+ ] of c∞ = 0.1 μM (corresponding to an average of ∼6 Ca2+ ions), the hitting time distributions for both fast and slow [Ca2+ ] fluctuations also appear exponential (see Figure 1A, left). When shown on a logarithmic scale (right), for [Ca2+ ] fluctuations much faster than the timescale of the Ca2+ -dependent

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

503

Figure 2: [Ca2+ ] fluctuations can increase hitting time and variability. (A) The hitting time mean (left) and coefficient of variation (CV, right) are shown as a function of [Ca2+ ] fluctuations timescale parameter τe for different microdomain volumes V. (B) Hitting time mean and CV are shown as a function of V for different values of τe . Thin dashed black lines in panel A (bottom) and panel B are the small microdomain volume approximations, calculated using equations 3.2 and 3.3. Parameter: c∞ = 0.1 μM.

signaling reaction S1 → S2 (τe = 10−2 , red), the hitting time distribution is nearly identical to that as when [Ca2+ ] fluctuations are absent. However, when [Ca2+ ] fluctuations are much slower than S1 → S2 (τe = 102 , blue), the probability of long hitting times increases slightly, that is, the distribution is long-tailed. The long tail is more prominent when fluctuations and S1 → S2 are on the same timescale (τe = 1, green). In a smaller microdomain of volume, V = 0.01 μm3 (corresponding to average of ∼0.6 Ca2+ ions), long hitting times are more likely when [Ca2+ ] fluctuations are on the same timescale or slower than S1 → S2 . This can be explained by the fact that due to the near-zero expected number of Ca2+ ions, zero Ca2+ ions may be present in the microdomain, and because of the slow timescale of passive exchange with the bulk, the time until one Ca2+ ion reenters the microdomain and subsequently triggers the S1 → S2 reaction may be relatively long. In Figure 2, we plot the hitting time mean μˆ (see equation 2.4) and coefficient of variation cˆv , normalized by the mean and CV, respectively, of

504

S. Weinberg

the model without [Ca2+ ] fluctuations (which for this single-step cascade are both equal to 1), as functions of the microdomain volume V and [Ca2+ ] fluctuations timescale parameter τe . For volumes of V = 0.1 μm3 to 1 μm3 , μˆ and cˆv are biphasic functions of τe . When τe is either very small or large, μˆ and cˆv are near 1—[Ca2+ ] fluctuations do not alter the hitting time— which we explain as follows. When [Ca2+ ] fluctuations are much faster than the signaling timescale (small τe ), even though the fluctuating [Ca2+ ] is typically above or below its expected value of c∞ , the duration of these events is very brief due to passive exchange with the bulk. [Ca2+ ] fluctuations are essentially smoothed at the longer timescale, as proposed by Hake and Lines (2008), such that the timing of Ca2+ -triggered signaling depends on only the average [Ca2+ ], c∞ . In contrast, when [Ca2+ ] fluctuations are much slower than the signaling timescale (large τe ), [Ca2+ ] fluctuations are negligibly influential, since [Ca2+ ] is effectively constant on the timescale of Ca2+ -triggered signaling, and again hitting time is minimally altered. However, for τe near 1—the timescale of [Ca2+ ] fluctuations and S1 → S2 are similar—μˆ and cˆv are increased. In this case, when the expected number of Ca2+ ions, c∞V, is small and fluctuations drive [Ca2+ ] above or below this value, the likelihood of zero microdomain Ca2+ ions is nonnegligible, such that the duration of zero microdomain Ca2+ ion events is significant at the timescale of the signaling event. Therefore, longer hitting times become more likely, as seen in the long tail in the hitting time distribution, which results in the expected hitting time and variability increasing. As microdomain volume V increases, the expected number of Ca2+ ions increases, reducing the likelihood of zero Ca2+ ions such that μˆ and cˆv are increased to a lesser extent. When microdomain volume is smaller (V = 10−3 to 10−2 μm3 ), both μˆ and cˆv increase as τe increases, a consequence of the near-zero expected number of Ca2+ ions in the microdomain, as discussed above. Note the different scales in the top and bottom plots in Figure 2A. In Figure 2B, μˆ and cˆv are shown as a function volume V. μˆ increases several orders of magnitude when V < 0.05 μm3 , which corresponds to an average of ∼3 Ca2+ ions, and cˆv is largest for microdomain volumes near this value. In Figure S2, we show that [Ca2+ ] fluctuations driven by Ca2+ -buffer association and dissociation yield results qualitatively similar compared with fluctuations driven by passive exchange. Gillespie-type Monte Carlo simulations (Gillespie, 1977) illustrate the timescales of [Ca2+ ] fluctuations for different values of τe and microdomain volume V (see Figure 3). As described above, simulations demonstrate that when [Ca2+ ] fluctuations are slow (τe = 102 ) and microdomain volume V = 0.1 μm3 is larger (see Figure 3A), [Ca2+ ] is effectively constant until the transition from S1 to S2 . For a smaller microdomain volume (V = 0.01 μm3 ; see Figure 3B), there are long time periods of zero microdomain Ca2+ ions, resulting in an increase in hitting time.

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

505

Figure 3: Monte Carlo simulations of microdomain [Ca2+ ] and system state. Representative Monte Carlo simulations of system state S and [Ca2+ ] are shown for different values of the [Ca2+ ] fluctuations timescale parameter τe . Simulations terminate when the system transition from state 1 to state 2 (i.e., the reaction S1 → S2 ) occurs. Parameters: Microdomain volume V = 0.1 (A) and 0.01 (B) μm3 , c∞ = 0.1 μM.

The parameter regime for which the dependence of μˆ and cˆv on τe transitions from biphasic to monotonically increasing is on the order of V = 0.01 to 0.1 μm3 . However, in multistep signaling cascades (discussed in the next section), multiple Ca2+ binding steps, each of which decreases the number of Ca2+ ions, increase the likelihood of zero Ca2+ ions present in the microdomain, including volumes at the large end of this range (V = 0.1 μm3 ). Thus, we find that when the number of cascade steps n ≥ 4 − 5, μˆ and cˆv are monotonically increasing functions of τe , and slow [Ca2+ ] fluctuations increase μˆ and cˆv (see Figure S3). For the remainder of the study, we consider microdomain volumes at the small end of this range, V = 0.01 μm3 , to calculate the largest extent for which [Ca2+ ] fluctuations may alter Ca2+ dependent signaling. For larger microdomain volumes, we find similar qualitative dependence on parameters, with μˆ and cˆv typically increased to lesser extents (not shown). In the appendix, we derive analytical approximations for the hitting time mean and variability, μˆ small and cˆsmall , respectively, in the limit of a small v microdomain volume V, specifically c∞V  1, for the single-step signaling

506

S. Weinberg

cascade, relating hitting time statistics (μ, ˆ cˆv ) to key model parameters (c∞ , V, τe ),  μˆ small = 1 + c∞V + τe

 1 −1 , c∞V

(3.2)

= σˆ small /μˆ small , where and cˆsmall v  σˆ

2,small

= 1 + c∞V (c∞V + 2) +

τe2

 2 1 −1 . + τe c∞V (c∞V )2

(3.3)

From equation 3.2, it is clear that in this limit, μˆ small > 1, that is, the expected hitting time is always increased. Further, since c∞V < 1, μˆ small also always increases as τe increases. We plot μˆ small and cˆsmall in Figure 2 (thin dashed black lines, bottom of v panel A and panel B) for values of V for which c∞V < 1, and we find that the small microdomain limit approximations agree nicely with the numerical calculations, with closer agreement for smaller V and larger values of τe . Further, it is straightforward to show that limV→0 cˆsmall = 1, in agreement v with numerical calculations. In this small microdomain volume limit, we also find that c∞ , V, τe , and μˆ small are related by limV→0 (c∞V μˆ small ) = τe . 3.1.2 Irreversible Signaling Cascade. We briefly investigate the influence of [Ca2+ ] fluctuations in the general, multistep Ca2+ -dependent signaling cascades. For simplicity, we consider an irreversible cascade, λ1

λ2

λn−1

λn

S1 → S2 → · · · → Sn → Sn+1 ,

(3.4)

λ−i = 0 and λi = 1 for all i, for which without [Ca2+ ] fluctuations, the hitting time distribution is well known, the hypoexponential distribution,  τ ∼ Hypo(λ1 , . . . , λn ) (Breuer & Baum, 2005), with mean μ = ni λ−1 and i n −2 2 variance σ = i λi . When transition rates between all states are equal, that is, λ = λ√1 = · · · = λn , τ ∼ Erlang(n, λ) with mean μ = n/λ, σ 2 = n/λ2 , and cv = 1/ n. For more general signaling cascades, including the SVR model, for which the hitting time distributions and moments cannot be described by simple analytical expressions, we use equations 2.3 to 2.5 to calculate these measures. In Figure 4, we plot the hitting time distribution for an n = 4 step signaling cascade, for different values of τe . In the absence of fluctuations (dashed line), as noted above, the hitting time distribution is hypoexponential. When [Ca2+ ] fluctuations are fast (small τe ), the hitting time distribution is similar to the distribution for the model without fluctuations. As τe increases, the

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

507

Figure 4: [Ca2+ ] fluctuations alter hitting time in an irreversible signaling cascade. (A) The hitting time distribution on a linear (left) and logarithmic (right) scale for different values of [Ca2+ ] fluctuations timescale parameter τe (solid color lines) or no fluctuations (dashed black line) is shown with time normalized to the hitting time mean in the absence of fluctuations μ. Parameters: λ1 = · · · = λn = 1, V = 0.01 μm3 , c∞ = 0.1 μM, n = 4.

Figure 5: [Ca2+ ] fluctuations alter hitting time in an irreversible signaling cascade. (A) The normalized expected hitting time μ, ˆ hitting time CV cv , and normalized hitting time CV cˆv are shown as functions of τe for different cascade lengths n. (B) μ, ˆ cv , and cˆv are shown as a function of n for different values of τe . Note that the black line in the top row (n = 1) is the same as the green line in Figure 2A. Thin dashed black lines in panel B are the small microdomain volume approximations, calculated using equations 3.5 and 3.6. Parameters: λ1 = · · · = λn = 1, V = 0.01 μm3 , c∞ = 0.1 μM.

distribution becomes right-shifted with a long tail, leading to an increase in the normalized expected hitting time μ. ˆ In Figures 5A and 5B, we plot μˆ (left), cv (middle), and cˆv (right) as functions of τe and the cascade length n. In general, we find that when [Ca2+ ] fluctuations are fast, μˆ and cˆv are near 1, as in the single-step cascade. However, when [Ca2+ ] fluctuations occur on a similar timescale or slower (i.e., τe ≥ 1), μˆ increases, more so as n increases, approaching an asymptotic

508

S. Weinberg

value around n = 5. cv decreases as n decreases, while cˆv is increased for large τe and approaches 1 as n increases. As discussed above, the right-shift and long tail in the hitting time distribution, leading to an increase in μˆ and cˆv for slow [Ca2+ ] fluctuations, are a consequence of the increased likelihood of zero microdomain Ca2+ ions, which results from multiple Ca2+ binding steps in the signaling cascade. For the general, irreversible cascade, we also derive the small microdomain volume approximation, shown in the appendix:  μˆ small = 1 + c∞V + τe

1 1 − c∞V n

 (3.5)

and cˆsmall = σˆ 2,small /μˆ small , where v  σˆ

2,small

= 1 + c∞V (c∞V + 2) +

τe2

 2 1 1 − + , τe c∞V (c∞V )2 n

(3.6)

and equations 3.2 and 3.3 are equivalent to equations 3.5 and 3.6 for n = 1. = 1, and as the number of cascade steps increases, As before, limV→0 cˆsmall v  lim cˆsmall n→∞ v

=

(c∞V )2 (c∞V + 1)2 + τe (2c∞V + τe ) , c∞V (c∞V + 1) + τe

(3.7)

which approaches 1 for both large and small τe . For intermediate τe near 1, this limit is less than 1, although cˆsmall is still near 1 for values of n ≤ 20. v We plot μˆ small and cˆsmall in Figure 5B (thin dashed black lines), and as in the v single-step cascade, we find nice agreement with numerical calculations, with closer agreement for larger τe . Our general approach to analyze the influence of the [Ca2+ ] fluctuations timescale can be applied to any model represented as a signaling cascade and can be utilized to investigate specific cascade properties of interest. In the online supplement, we analyze in more detail the influence of several cascade properties, include reversibility (see Figure S4), cooperativity (see Figure S5), and a fast or slow terminal cascade step (see Figure S6) and found that changes in hitting time mean and variability, as a function of [Ca2+ ] fluctuations timescale parameter τe , were highly dependent on these cascade properties and with complex dependence. 3.2 Synaptic Vesicle Release. Following our analysis showing that [Ca2+ ] fluctuations alter timing in a general signaling cascade, we next analyze the influence of [Ca2+ ] fluctuations on SVR time. Experiments have shown that SVR occurs following sequential Ca2+ binding steps preceding a final Ca2+ -independent step (Heidelberger et al., 1994; Schneggenburger & Neher, 2000), suggesting that SVR can be modeled as a specific case

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

509

of the general signaling cascade studied above. Here, we consider a wellcharacterized model of Ca2+ -dependent SVR (Bollmann, 2000), where the release time is equivalent to the hitting time for the terminal step in the cascade: 5ac

4ac

3ac

2ac

ac

γ

b

2b

3b

4b

5b

δ

pν

X  XCa1  XCa2  XCa3  XCa4  XCa5  XCa∗5 → F. (3.8) In this signaling cascade, the first five steps represent Ca2+ sequentially and reversibly binding to five independent binding sites on Ca2+ -sensor molecule X, with association rate constant a (with units of concentration−1 × time−1 ) and dissociation rate b (with units of time−1 ), related by dissociation constant Kd = b/a (units of concentration). The fully Ca2+ -bound sensor can then undergo a reversible Ca2+ -independent activation step, with forward and backward rates γ and δ (both with units of time−1 ), respectively, and a terminal step that represents the triggering of synaptic vesicle release from the activated state, with rate pν, where p is a release rate per available vesicle and ν is the number of available vesicles for release. For reasonable comparisons of calculations for different values of the average [Ca2+ ], c, we define the passive exchange rate relative to Ca2+ binding at rest conditions and the [Ca2+ ] fluctuations parameter τe , such that kˆ e = ac∞ /τe . Note that physiological values of τˆe = 1/kˆ e on the order of 10−2 − 102 ms (Weinberg & Smith, 2014; von Wegner et al., 2014) correspond to τe values on the order of 10−4 − 10. As in the previous sections, we present analysis for τe ranging from 10−2 − 102 because there are negligible differences between measures for τe less than 10−2 , that is, faster [Ca2+ ] fluctuations. We also present analysis for slower [Ca2+ ]fluctuations, that is, τe values an order of magnitude greater than the upper limit of this range to τe = 102 , as these τe values may be appropriate in pathological settings that slow or impair diffusion. 3.2.1 Numerical Calculations of Vesicle Release Timing for Constant Average [Ca2+ ]. Using the approach presented above, we plot the distributions for the release timing for different values of τe . In both the absence of [Ca2+ ] fluctuations (dashed line) and fast [Ca2+ ] fluctuations (red, magenta, green), the distribution is approximately exponential (see Figure 6A). As [Ca2+ ] fluctuations become slower and τe increases (cyan, blue), the distributions become more long tailed and right shifted, and μˆ increases, as in the previous general multistep example. We next analyze the influence of the [Ca2+ ] fluctuations timescale parameter τe on the release time for varying levels of the average [Ca2+ ], c. For a given [Ca2+ ] level, μˆ monotonically increases as τe increases, and for some values of c, μˆ approaches an asymptotic value, while cˆv is only slightly increased for some parameters (see Figure 6B). For a given value of τe , μˆ

510

S. Weinberg

Figure 6: [Ca2+ ] fluctuations alter synaptic vesicle release timing. (A) SVR time distribution on a linear (left) and logarithmic (right) scale for different values of [Ca2+ ] fluctuations timescale parameter τe (solid color lines) or no fluctuations (dashed black line) are shown with time normalized to the release time mean in the absence of fluctuations μ. Note that the red and magenta lines overlap with the green line. (B) The normalized expected release time μ, ˆ release time CV cv , and normalized release time CV cˆv are shown as functions of τe for ˆ cv , and cˆv are shown as functions different values of average [Ca2+ ], c. (C) μ, of c for different values of τe . Parameters (Bollmann, 2000): a = 0.3 μM−1 ms−1 , b = 3 ms−1 , δ = 8 ms−1 , γ = 30 ms−1 , p = 40 ms−1 vesicles−1 , ν = 800 releasable vesicles, V = 0.01 μm3 , c∞ = 0.1 μM, kˆ e = ac∞ /τe . In panel A, c = 1 μM.

is a biphasic function of c, such that for c near resting levels c∞ and very elevated near 100 μM, μˆ is not increased and near 1 (see Figure 6C), while for intermediate concentrations, between approximately 1 and 10 μM, μˆ is greater than 1 and generally increases as τe increases. cˆv exhibits a more complex triphasic dependence on c, but in general is also increased over a similar range of [Ca2+ ] values, as is μ. ˆ We next consider the influence of varying several SVR model parameters. We find that for all values of dissociation constant Kd , μˆ increases as τe

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

511

Figure 7: Ca2+ binding affinity alters the influence of [Ca2+ ] fluctuations on synaptic vesicle release timing. The normalized expected release time μ, ˆ release time CV cv , and normalized release time CV are shown as functions of average [Ca2+ ], c, for different values of [Ca2+ ] fluctuations timescale parameter τe and Ca2+ dissociation constant Kd . Note that the plots in the third row (Kd = 10 μM) are the same as in Figure 6C. Parameters: a = 0.3 μM−1 ms−1 , b = aKd . Other parameters as in Figure 6.

increases (see Figure 7), as in Figure 6. Slow [Ca2+ ] fluctuations increase the expected SVR time, and this increase is most prominent over a range of c values, which is shifted to higher values of c as Kd increases (decreased binding affinity). μˆ is increased to a larger extent for small Kd , that is, [Ca2+ ] fluctuations increase expected SVR time to a larger extent when binding affinity is high, consistent with our findings in a general multistate reversible cascade (see Figure S4B). This increase in expected SVR time can be quite large, between one and two orders of magnitude. We find a similar response for cv and cˆv , with cˆv increased up to a factor of 2 and over a similar range for c values. We next consider the influence of varying the number of Ca2+ binding sites, , in the sensor molecule X, that is, analyzing the variable-length SVR signaling cascade given by ac

( −1)ac

b

2b

X  XCa1



···

2ac



( −1)b

ac

γ

b

δ

pν

XCa −1  XCa  XCa∗ → F.

(3.9)

512

S. Weinberg

Figure 8: The number of Ca2+ binding sites alters the influence of [Ca2+ ] fluctuations on synaptic vesicle release timing. The normalized expected release time μ, ˆ release time CV cv , and normalized release time CV are shown as functions of average [Ca2+ ], c, for different values of [Ca2+ ] fluctuations timescale parameter τe and number of Ca2+ binding sites . Note that the plots in the fourth row ( = 5) are the same as in Figure 6C. Parameters as in Figure 6.

Interestingly, we find that increasing or decreasing the number of binding sites has a similar effect as that of increasing or decreasing Kd (see Figure 8). Decreasing shifts the range of c values and increase the extent for which fluctuations increase the expected SVR time and variability, μˆ and cˆv , respectively. We also vary the number of available vesicles, ν, and find that increasing or decreasing ν by an order of magnitude (80–8000) from the baseline value of 800 negligibly influenced release time statistics (μˆ and cˆv differences of < 0.4% over the entire parameter space for τe and c investigated in Figure 6), as this terminal step is still much faster than the preceding cascade steps in this parameter regime (not shown). We show Monte Carlo simulations of SVR to illustrate the dependence on the average [Ca2+ ], dissociation constant Kd , and [Ca2+ ] fluctuations timescale parameter τe when τe = 102 is large (see Figure 9). Although the passive exchange rate kˆ e is constant in these simulations, varying c clearly alters the timescale of [Ca2+ ] fluctuations, such that fluctuations are faster for smaller c. For large τe , both the association and dissociation rate, ac and b,

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

513

Figure 9: Monte Carlo simulation of synaptic vesicle release for a constant [Ca2+ ]. Representative Monte Carlo simulations of the system state (top) and [Ca2+ ] (bottom) for different values of the average [Ca2+ ], c. Parameters: τe = 102 . Other parameters as in Figure 6.

respectively, are much faster than passive exchange rate kˆ e , such that [Ca2+ ] fluctuations are primarily driven by Ca2+ binding to and unbinding from the sensor molecule, not passive exchange. For large values of c near 100 μM, such that c > Kd = 10 μM and the association rate is faster than dissociation, ac > b, SVR occurs relatively quickly, such that [Ca2+ ] is effectively constant over the timescale of the signaling cascade and [Ca2+ ] fluctuations are relatively slow on this timescale. As in Figure 3, when [Ca2+ ] is effectively constant, μˆ and cˆv in general are not increased. In contrast, for c = c∞ < Kd , such that the dissociation rate b > ac is faster, SVR occurs relatively slowly, [Ca2+ ] fluctuations are relatively fast on the timescale of the signaling cascade, and μˆ and cˆv are also not increased. However, for intermediate values of c near Kd , the association and dissociation reactions are on a comparable timescale, driving [Ca2+ ] fluctuations of comparable timescale as the signaling cascade; as a consequence, μˆ and cˆv are increased. When Kd is shifted, the range of values of c for which the release time is increased is similarly shifted (see Figure 7). When c is both near Kd and small, the magnitude of [Ca2+ ] fluctuations is relatively large and increases μˆ and cˆv to a larger extent, compared with both large c and Kd and relatively small fluctuations. As τe decreases, [Ca2+ ] fluctuations are driven to a larger extent by passive exchange, which leads to faster [Ca2+ ] fluctuations, which in turn results in μˆ and cˆv increasing to a smaller extent. 3.2.2 Monte Carlo Simulation of Vesicle Release Timing in the Spiking Neuron. Finally, we perform Monte Carlo simulations of SVR in a spiking neuron to demonstrate the significance of [Ca2+ ] fluctuations in the setting of a

514

S. Weinberg

Figure 10: Monte Carlo simulation of synaptic vesicle release in a spiking neuron. (A) Transmembrane potential (Vm ), Ca2+ current (ICa ), and [Ca2+ ], for different values of [Ca2+ ] fluctuations timescale parameter τe , are shown as a function of time. Solid lines below each [Ca2+ ] trace denote individual SVR events. Insets illustrate fluctuations in [Ca2+ ] that occur between neuron spikes. (B) Histogram for SVR time for the three models described in the text, for different values of τe : model 1 (color), 2 (gray), and 3 (black). (C) Top: SVR time mean (left) and CV (right) are shown as functions of τe for the three models: model 1 (thick black), 2 (gray), and 3 (black). Bottom: SVR time mean and CV for model 1, normalized by the respective value in models 2 (gray) and 3 (thin black), are shown as functions of τe . SVR parameters as in Figure 6.

time-varying Ca2+ influx. We simulate the Hodgkin-Huxley neuron, augmented with a high-threshold Ca2+ channel current. A constant current stimulus is applied, such that the neuron spikes repetitively at approximately 70 Hz, and the Ca2+ current triggers Ca2+ influx and subsequent Ca2+ -triggered SVR (see Figure 10A, top). To illustrate the significance of [Ca2+ ] fluctuations, we consider three model formulations: (1) the stochastic SVR/Ca2+ model, in which both SVR and Ca2+ influx are stochastic; (2) the stochastic SVR model, in which [Ca2+ ] is deterministic and governed by the Ca2+ current and mass-action kinetics; and (3) the deterministic model, in which both SVR and [Ca2+ ] are governed by mass-action kinetics. Further details of the model and simulations methods are provided in the appendix. In Figure 10A, [Ca2+ ] is shown as a function of time for different values of τe . For small τe , Ca2+ transients are shorter due to faster passive exchange, such that SVR events (solid lines below each trace) typically occur during

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

515

neuron spikes. For larger τe , Ca2+ transients are longer, such that SVR events also occur between spikes, although at reduced frequency. Insets illustrate [Ca2+ ] fluctuations between spikes. In Figure 10B, histograms of SVR time are shown for different values of τe . For small τe (model 1; red), the distribution is bimodal, with the larger peak corresponding to the time between spikes and the second peak corresponding to the timing between SVR events during a single spike (i.e., intraspike SVR). As τe increases and SVR events occur between spikes, this first peak is reduced in magnitude until it is no longer present and the second peak is left-shifted as the time between SVR events is reduced (model 1; green, blue). We compare these distributions with SVR timing calculated from the models in which [Ca2+ ] fluctuations are not present. In a model of deterministic SVR (model 3, black), the shift from a bimodal to unimodal distribution as τe increases is reproduced. However, the spread of each peak is greatly reduced, demonstrating that the deterministic model does not fully reproduce SVR time variability. Further, the peak corresponding to intraspike SVR events is right-shifted. In a model of stochastic SVR and deterministic Ca2+ influx (model 2, gray), the SVR time distribution is left-shifted relative to model 1 (i.e., neglecting [Ca2+ ] fluctuations underestimates SVR timing). We plot the SVR time mean and CV as functions of τe in Figure 10C. For all three models, the expected SVR time decreases as τe increases. For all values of τe , expected SVR time in model 1 (thick black) is greater than model 2 (gray) and less than model 3 (thin black). Importantly, we find that SVR time mean is underestimated when neglecting [Ca2+ ] fluctuations (see Figure 10C, bottom left, gray) and overestimated when neglecting both SVR and [Ca2+ ] fluctuations (black). SVR time CV is also a biphasic function of τe , largest typically near τe = 1 for model 1, and for most values of τe is largest when accounting for [Ca2+ ] fluctuations. For most values of τe , SVR time CV is also underestimated when neglecting [Ca2+ ] fluctuations (see Figure 10C, bottom right), demonstrating that [Ca2+ ] fluctuations modulate SVR release time and potentially subsequently spike timing and variability. 4 Discussion and Conclusion 4.1 Summary of Main Findings. In this study, we demonstrate that [Ca2+ ] fluctuations can greatly influence the hitting time for a general Ca2+ dependent signaling cascade and specifically in a model of Ca2+ -triggered SVR. In a general signaling cascade, small microdomain volume V and multiple Ca2+ binding steps can lead to an increased likelihood of zero microdomain Ca2+ ions, such that [Ca2+ ] fluctuations on the same timescale or slower than a Ca2+ -dependent signaling cascade greatly increase the hitting time mean and variability due to an increased likelihood of very long hitting times (i.e., a long-tailed distribution). However, when V or [Ca2+ ] is sufficiently large, such that the likelihood of zero Ca2+ ions is

516

S. Weinberg

small, slow [Ca2+ ] fluctuations result in an effectively constant [Ca2+ ] and minimal influence on hitting times. Similarly, when [Ca2+ ] fluctuations due to passive exchange are faster (small τe ), even when V is small, [Ca2+ ] fluctuations also negligibly alter hitting times, since the Ca2+ -dependent signaling effectively responds to the average [Ca2+ ]. In an SVR model, [Ca2+ ] fluctuations may be driven by both passive exchange and Ca2+ association and dissociation reactions. When the resulting fluctuations are on a comparable timescale, SVR time mean and variability may be greatly increased (see Figure 6). The influence of [Ca2+ ] fluctuations highly depends on Ca2+ binding affinity Kd (see Figure 7) and the number of Ca2+ binding sites (see Figure 8) but negligibly depends on the number of vesicles available for release. Monte Carlo simulations similarly demonstrate that [Ca2+ ] fluctuations can greatly increase SVR time mean and variability in a spiking neuron (see Figure 10). Importantly, these points highlight that both passive exchange and Ca2+ dependent reactions are significant intrinsic sources of [Ca2+ ] fluctuations that can alter the temporal dynamics of Ca2+ -dependent signaling cascades, including SVR and, further, that a physiological model that neglects the influence of [Ca2+ ] fluctuations may be dramatically underestimating the timing of these Ca2+ -mediated signaling events. 4.2 Comparison with Prior Simulation Studies. In this study, we demonstrate in both general and physiological settings that [Ca2+ ] fluctuations may or may not influence the timing of Ca2+ -dependent signaling. In many parameter regimes, our study agrees with the findings of Hake and Lines (2008) that [Ca2+ ] fluctuations are smoothed by signaling events occurring on a longer timescale; when [Ca2+ ] fluctuations are much faster than the signaling event (small τe ), the hitting time is not increased or decreased. However, we find that this is not always the case, specifically in the setting of a very small microdomain volume V (see Figures 1 and 2). The reason for the disagreement between our results and the analysis of [Ca2+ ] fluctuations in the cardiac dyad by Hake and Lines (2008) may arise from the specific details of the models. The stochastic simulations they performed are three-dimensional random walks in a small volume of ∼5 × 10−4 μm3 . For a subspace of such small volume, we would expect [Ca2+ ] fluctuations to greatly influence Ca2+ -dependent signaling hitting times. However, in their study, the dyad geometry is modeled as an ideal cylinder, with the top and bottom representing membranes and the side representing the region for escape to the bulk. However, the dyad structure may be more appropriately described as a labyrinth (Cannell et al., 2013), in which the sarcoplasmic reticulum membrane defines a complex geometry in the dyad. In our study, while we do not specify a specific geometry, we treat the Ca2+ microdomain as a spatially restricted (well-mixed) subcellular compartment.

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

517

In general, spatially localized Ca2+ signaling can arise as a consequence of complex membrane configurations and other obstructions to diffusion (e.g., molecular crowding). This compartmental modeling approach is appropriate when the characteristic time for microdomain escape τˆx is much longer than the characteristic time associated with diffusion across (within) the domain τˆdi f f . Schuss et al. (2007) showed that the timescale for microdomain escape via a small opening (“escape window”) of radius r, τˆx = V/(4rD), where D is the diffusion coefficient. For three-dimensional diffusion, τˆdi f f = L2 /(6D), where L is a length characteristic of the microdomain and if we assume V = L3 , then the ratio of these timescales τˆx /τˆdi f f = 32 L/r becomes large as the escape radius decreases. The validity of the compartmental modeling approach will highly depend on subspace geometry. Tanskanen et al. (2007) showed that [Ca2+ ] fluctuations are influential when accounting for the complex geometry of the cardiac dyad, in agreement with our conclusion that [Ca2+ ] fluctuations are significant in spatially restricted microdomains and specifically the presynaptic terminal. Collectively, these studies and our results suggest that accounting for both [Ca2+ ] fluctuations and microdomain geometry is important. 4.3 Physiological Significance. In this study, we demonstrate that [Ca2+ ] fluctuations can significantly alter the timing of SVR and further that models of SVR neglecting the influence of [Ca2+ ] fluctuations may underestimate the timing between release events. Much work has investigated sources of neuronal spike timing and variability over a wide range of spatial and temporal scales, including the role of stochastic ion channel gating (Anwar, Hepburn, Nedelescu, Chen, & De Schutter, 2013) and noisy spike thresholds (Coombes, Thul, Laudanski, Palmer, & Sumner, 2011). Our study suggests that stochasticity in SVR, and specifically modulation by [Ca2+ ] fluctuations, may also play a significant role in altering the timing and variability of neural firing, especially within larger models of synaptically coupled neurons or networks (Anderson, Kudela, Weinberg, Bergey, & Franaszczuk, 2009; Weinberg, 2015). Interestingly, synaptic release timing is increased to the greatest extent when the average [Ca2+ ] is slightly elevated above resting levels, near 1 to 10 μM. Dynamics at such moderately elevated Ca2+ levels may be significant physiologically, as Ca2+ -dependent signaling may trigger feedback mechanisms that either suppress or further increase Ca2+ influx (e.g., Ca2+ -activated or Ca2+ -inactivated Ca2+ channels). In many physiological settings, Ca2+ binding proteins or buffers can play a significant role in the spatiotemporal dynamics of Ca2+ signaling. For example, buffers can alter the decay of residual Ca2+ following Ca2+ channel closure and influence release dynamics of a Ca2+ release cluster (Dargan, Schwaller, & Parker, 2004). Buffers may play the additional role of modulating the timescale of [Ca2+ ] fluctuations (Weinberg & Smith, 2014), and in some settings, driving fast [Ca2+ ] fluctuations (i.e., small τe ) and

518

S. Weinberg

reducing the timing variability in Ca2+ -dependent signaling. Of course, the influence of buffers may be complex, as Ca2+ binding proteins such as calmodulin can both trigger or modulate Ca2+ -dependent processes and drive [Ca2+ ] fluctuations. While we investigate Ca2+ -dependent signaling cascades in this study, our approach is broadly applicable to essentially all signaling cascades mediated by an effector present in small copy number in spatially restricted domains or subcellular regions, such as protein kinase cascades (Herskowitz, 1995) and transcription factors with multiple sequential targets (Alon, 2006). Holcman, Dao Duc, Jones, Byrne, and Burrage (2014) have extended the general two-dimensional Markov chain approach to investigate the influence of the small number of messenger RNA molecules on the mean first passage time in the setting of post-transcriptional nuclear regulation. Further, while we focus on signaling cascades, that is, reaction networks in which signaling events occur in a sequential order, due to their prevalence in physiological signaling, our approach is applicable to all biochemical reaction networks. It is straightforward to enumerate the state-space for any biochemical reaction network, define the transition between these states (Higham, 2008), and characterize the timing for any specific state of physiological interest. However, as the size of reaction network increases, numerical calculations may be limited by the size of the state-space, discussed in the next section. 4.4 Limitations. The hitting time for the final stage in a signaling cascade is random, and the approach used in this study enables the direct numerical calculation of hitting time distribution, mean, and variability, without Monte Carlo simulation, for a wide range of parameters and cascade properties. Estimation of the hitting time via simulation in the presence of fast [Ca2+ ] fluctuation requires many more iterations of the Gillespie stochastic simulation algorithm (SSA), compared with slow [Ca2+ ] fluctuations, due to many more passive exchange events. In contrast, the size of the Q matrix used in equations 2.3 to 2.5 is independent of the fluctuations timescale, although convergence in numerical algorithms will be parameter dependent. This facilitated a wide range of parameter studies. However, numerical calculations of a matrix exponential, inversion, and multiplication (see equations 2.3 to 2.5) become computationally expensive as the size of the state-space increases. For example, if the average [Ca2+ ] is 10 μM and microdomain volume V = 0.1 μm3 , then Cm = 1204, and the size of the state space (i.e., the dimension of the Q matrix) for a signaling cascade of n = 10 steps, N = n(Cm + 1) = 12,050, is large, such that Q is a (sparse) matrix with approximately 1.45 · 108 elements. For an exceedingly large state space, numerical calculations will be infeasible, and Monte Carlo simulations will be necessary to analyze hitting time statistics.

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

519

For simplicity in this study, the fluctuation timescale is governed solely by passive exchange (which can be defined by a single parameter). We (Weinberg & Smith, 2014) and others (von Wegner et al., 2014) have shown that the timescale of [Ca2+ ] fluctuations is influenced by diffusion, microdomain geometry, and buffer properties. In this study, we demonstrated that both passive exchange and Ca2+ -dependent association and dissociation are potentially significant sources of [Ca2+ ] fluctuations that can alter hitting times, and we focus on the timescale of [Ca2+ ] fluctuations driven by passive exchange, which enabled the investigation of a wide range of model parameters, for both the general signaling cascade and the synaptic vesicle release model. However, further analysis is necessary to understand whether additional sources of [Ca2+ ] fluctuations (i.e., buffers), competing Ca2+ -binding reactions, and other intrinsic or extrinsic sources, occurring on potentially multiple timescales, are significant in the context of Ca2+ -dependent signaling. In the online supplement, we show that there was minimal difference when fluctuations are driven by immobile buffer (see Figure S2). However, accounting for a mobile buffer would increase the size of the state space by potentially many orders of magnitude, and thus this analysis would face the computational challenges discussed above. Appendix: Additional Model Details, Methods, and Approximations A.1 Small Volume Approximation. We derive analytical approximations for μˆ and cˆv in the limit of a small microdomain volume V for the single-step signaling cascade. As V decreases, the expected number of microdomain Ca2+ ions, c∞V, also decreases, such that for sufficiently small V, the likelihood of more than one Ca2+ ion present in the domain becomes negligible: c∞V  1. Therefore, we can define the maximum value for C in the CTMC described by equation 2.8, Cm = 1, which simplifies Q considerably. For the case of the single irreversible reaction S 1 → S2 ,  T=

−c∞V/τe

c∞V/τe

1/τe

−1/(c∞V ) − 1/τe

 ,

(A.1)

and ζ = (1 − c∞V, c∞V ). We then use equations 2.4 and 2.5 to analytically derive equations 3.2 and 3.3, respectively. For the general, irreversible cascade S1 → S2 → · · · → Sn → Sn+1 , we can also derive the small microdomain volume approximation. Again assuming c∞V  1, we can define Cm = 1, simplifying Q, such that the 2n × 2n matrix, T, is given by

520

S. Weinberg

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ T=⎜ ⎜ ⎜ ⎜ ⎝

−c∞V/τe 1/τe 0 0 .. . 0 0 ··· ··· ··· ··· .. . 0 0

c∞V/τe −1/(c∞V ) − 1/τe 0 0 .. . 0 0

0 1/(c∞V ) −c∞V/τe 1/τe .. . 0 0 ⎞

0 0 0 0

0 0 0 0

1/(c∞V ) −c∞V/τe 1/τe

0 c∞V/τe −1/(c∞V ) − 1/τe

0 0 c∞V/τe −1/(c∞V ) − 1/τe .. . ··· ···

⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠

(A.2)

ζ = (1 − c∞V, c∞V, 0, . . . , 0), and then, as before, use equations 2.4 and 2.5 to analytically derive equations 3.5 and 3.6, respectively. A.2 Model Details and Monte Carlo Simulation Methods for Synaptic Vesicle Release in the Spiking Neuron. To simulate the stochastic timing of Ca2+ -triggered synaptic vesicle release, we incorporate the model of vesicle release, equation 3.8, into the classical Hodgkin-Huxley (HH) neuron model (Hodgkin & Huxley, 1952), augmented with a high-threshold calcium current ICa , given by the following system of differential equations: v˙ = (Iapp − INa − IK − IL − ICa )/C,

(A.3a)

x˙ = αx (1 − x) − βx x,

(A.3b)

where x ∈ {m, h, n}, and the ionic currents are given by INa = gNa m3 h(v − ENa ), IK = gK n4 (v − EK ), IL = gL (v − EL ), and ICa = gCa d∞ (v − ECa ). The equations governing the dynamics of the gating variables m, h, and n are given by αm =

2.5 − 0.1v , exp(2.5 − 0.1v) − 1

αh = 0.07 exp(−v/20), αn =

0.1(1 − 0.1v) , exp(1 − 0.1v) − 1

βm = 4 exp(−v/18),

βh =

1 , exp(3 − 0.1v) + 1

βn = 0.125 exp(−v/80),

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

521

where v is the shifted transmembrane potential, in which the resting potential Vrest = −80 mV has been subtracted. The instantaneous activation gating of ICa is given by d∞ =

   v − V1 1 1 + tanh , 2 V2

where V1 = 78.8 and V2 = 18 mV. Other parameter values are given as follows: maximum ionic current conductances (mS/cm2 ): gNa = 120, gK = 36, gL = 0.3; and reversal potentials, in which Vrest has been subtracted (mV): ENa = 115, EK = −12, EL = 10.6, ECa = 200. In all simulations, Iapp = 10 μA/cm2 . Since increasing [Ca2+ ] fluctuations parameter τe decreases the passive exchange rate, for a given Ca2+ influx, via ICa , increasing τe results in a larger [Ca2+ ]. Therefore, the value of gCa is decreased as τe is increased, such that the peak [Ca2+ ] in the deterministic model is fixed and approximately 70 μM. For τe = 10−2 , 1, and 102 , gCa = 2.2, 0.59, and 0.48 mS/cm2 , respectively. Monte Carlo simulations are performed using a modified version of the Gillespie algorithm that accounts for time-varying propensity functions, described by Alfonsi, Cances, Turinici, Ventura, and Huisinga (2005). In the stochastic SVR/Ca2+ model, we augment the set of reactions given in equation 2.7, with the zeroth-order reaction, JCa (t)

∅ → Ca2+ ,

(A.5)

where JCa (t) = AICa (t)/(2e) is the Ca2+ flux, in units of (ions · time−1 ), that can be calculated by multiplying ICa by the microdomain surface area A, and dividing by the elementary charge e and a factor of 2 for the divalent ion. For simplicity, we treat the presynaptic microdomain as a sphere, such that A = 4πr2d , where radius rd = [3V/(4π )]1/3 is calculated from the volume V. ICa (t) is determined by numerical integration of equation A.3. In the stochastic SVR model, [Ca2+ ](t) is determined by numerical integration of equation A.3, coupled with the corresponding mass action kinetic equations governed by the passive exchange and binding reactions in equation 2.7, and then Monte Carlo simulations are performed using the Alfonsi method with the corresponding [Ca2+ ](t). Acknowledgments I acknowledge and thank Old Dominion University High Performance Computing for use of the Turing cluster to perform numerical calculations and simulations.

522

S. Weinberg

References Alfonsi, A., Cances, E., Turinici, G., Ventura, B. D., & Huisinga, W. (2005). Adaptive simulation of hybrid stochastic and deterministic models for biochemical systems. ESAIM: Proceedings, 14, 1–13. Alon, U. (2006). An introduction to systems biology: Design principles of biological circuits. Boca Raton, FL: CRC Press. Anderson, W. S., Kudela, P., Weinberg, S., Bergey, G. K., & Franaszczuk, P. J. (2009). Phase-dependent stimulation effects on bursting activity in a neural network cortical simulation. Epilepsy Res., 84(1), 42–55. Anwar, H., Hepburn, I., Nedelescu, H., Chen, W., & De Schutter, E. (2013). Stochastic calcium mechanisms cause dendritic calcium spike variability. J. Neurosci., 33(40), 15848–15867. Berridge, M. J. (1998). Neuronal calcium signaling. Neuron, 21(1), 13–26. Berridge, M. J., Bootman, M. D., & Roderick, H. L. (2003). Calcium signalling: Dynamics, homeostasis and remodelling. Nat. Rev. Mol. Cell Biol., 4(7), 517–529. Biess, A., Korkotian, E., & Holcman, D. (2011). Barriers to diffusion in dendrites and estimation of calcium spread following synaptic inputs. PLOS Comput. Biol., 7(10), e1002182–14. Bladt, M. (2005). A review on phase-type distributions and their use in risk theory. Astin Bulletin, 35(1), 145–161. Bollmann, J. H. (2000). Calcium sensitivity of glutamate release in a calyx-type terminal. Science, 289(5481), 953–957. Bollmann, J. H., & Sakmann, B. (2005). Control of synaptic strength and timing by the release-site Ca2+ signal. Nat. Neurosci., 8(4), 426–434. Breuer, L., & Baum, D. (2005). An introduction to queueing theory and matrix-analytic methods. New York: Springer. Cannell, M. B., Kong, C.H.T., Imtiaz, M. S., & Laver, D. R. (2013). Control of sarcoplasmic reticulum Ca2+ release by stochastic RyR gating within a 3D model of the cardiac dyad and importance of induction decay for CICR termination. Biophys. J., 104(10), 2149–2159. Carrasco, M. A., & Hidalgo, C. (2006). Calcium microdomains and gene expression in neurons and skeletal muscle cells. Cell Calcium, 40(5–6), 575–583. Coombes, S., Thul, R., Laudanski, J., Palmer, A. R., & Sumner, C. J. (2011). Neuronal spike-train responses in the presence of threshold noise. Frontiers in Life Science, 5(3–4), 91–105. Dao Duc, K., & Holcman, D. (2010). Threshold activation for stochastic chemical reactions in microdomains. Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 81(4 Pt. 1), 041107. Dargan, S. L., Schwaller, B., & Parker, I. (2004). Spatiotemporal patterning of IP3mediated Ca2+ signals in Xenopus oocytes by Ca2+ -binding proteins. J. Physiol., (London), 556(Pt. 2), 447–461. Darvey, I. G., Ninham, B. W., & Staff, P. J. (1966). Stochastic models for second-order chemical reaction kinetics: The equilibrium state. J. Chem. Phys., 45(6), (1966) 2145– 2155. Flegg, M. B., Chapman, S. J., & Erban, R. (2012). The two-regime method for optimizing stochastic reaction-diffusion simulations. J. R. Soc. Interface, 9(70), 859–868.

Microdomain [Ca2+ ] Fluctuations Alter Synaptic Vesicle Release

523

¨ Flegg, M. B., Rudiger, S., & Erban, R. (2013). Diffusive spatio-temporal noise in a first-passage time model for intracellular calcium release. J. Chem. Phys., 138(15), 154103. Gaur, N., & Rudy, Y. (2011). Multiscale modeling of calcium cycling in cardiac ventricular myocyte: Macroscopic consequences of microscopic dyadic function. Biophys. J., 100(12), 2904–2912. Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. Journal Phys. Chem., 81(25), 2340–2361. Hake, J., & Lines, G. T. (2008). Stochastic binding of Ca2+ ions in the dyadic cleft: Continuous versus random walk description of diffusion. Biophys. J., 94(11), 4184– 4201. Harris, K. M., Jensen, F. E., & Tsao, B. (1992). Three-dimensional structure of dendritic spines and synapses in rat hippocampus (CA1) at postnatal day 15 and adult ages: Implications for the maturation of synaptic physiology and long-term potentiation. J. Neurosci., 12(7), 2685–2705. Heidelberger, R., Heinemann, C., Neher, E., & Matthews, G. (1994). Calcium dependence of the rate of exocytosis in a synaptic terminal. Nature, 371(6), 513–515. Herskowitz, I. (1995). MAP kinase pathways in yeast: For mating and more. Cell, 80(2), 187–197. Higham, D. J. (2008). Modeling and simulating chemical reactions. SIAM Rev., 50(2), 347–368. Hodgkin, A. L., & Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol., 117(4), 500–544. Holcman, D., Dao Duc, K., Jones, A., Byrne, H., & Burrage, K. (2014). Posttranscriptional regulation in the nucleus and cytoplasm: Study of mean time to threshold (MTT) and narrow escape problem. J. Math. Biol., 70(4), 805–828. Holcman, D., & Schuss, Z. (2005). Stochastic chemical reactions in microdomains. J. Chem. Phys., 122(11), 114710–114716. Keizer, J. (1987). Statistical thermodynamics of nonequilibrium processes. New York: Springer. Korkotian, E., & Segal, M. (2006). Spatially confined diffusion of calcium in dendrites of hippocampal neurons revealed by flash photolysis of caged calcium. Cell Calcium, 40(5–6), 441–449. Lamar, M. D., Kemper, P., & Smith, G. D. (2011). Reduction of calcium release site models via moment fitting of phase-type distributions. Phys. Biol., 8(2), 026015. McQuarrie, D. A., Jachimowski, C. J., & Russell, M. E. (1964). Kinetics of small systems, II. J. Chem. Phys., 40(10), 2914–2919. Modchang, C., Nadkarni, S., Bartol, T. M., Triampo, W., Sejnowski, T. J., Levine, H., & Rappel, W.-J. (2010). A comparison of deterministic and stochastic simulations of neuronal vesicle release models. Phys. Biol., 7(2), 026008. Morita, A. (1988). Externally driven fluctuation of concentration and dynamics of second order chemical reaction. J. Chem. Phys., 88(12), 7481–7484. Nadkarni, S., Bartol, T. M., Sejnowski, T. J., & Levine, H. (2010). Modelling vesicular release at hippocampal synapses. PLoS Comput. Biol., 6(11), e1000983. ¨ Nagaiah, C., Rudiger, S., Warnecke, G., & Falcke, M. (2012). Adaptive space and time numerical simulation of reaction-diffusion models for intracellular calcium dynamics. Applied Mathematics and Computation, 218, 10194–10210.

524

S. Weinberg

Neher, E., & Sakaba, T. (2008). Multiple roles of calcium ions in the regulation of neurotransmitter release. Neuron Review, 59(6), 861–872. ¨ ¨ Ruckl, M., Parker, I., Marchant, J. S., Nagaiah, C., Johenning, F. W., & Rudiger, S. (2015). Modulation of elementary calcium release mediates a transition from puffs to waves in an IP3R cluster model. PLOS Comput. Biol., e1003965–12. Schneggenburger, R., & Neher, E. (2000). Intracellular calcium dependence of transmitter release rates at a fast central synapse. Nature, 406, 889–893. Schuss, Z., Singer, A., & Holcman, D. (2007). The narrow escape problem for diffusion in cellular microdomains. Proceedings of the National Academy of Sciences, 104(41), 16098–16103. Smith, G. D. (2004). Modeling intracellular calcium: Diffusion, dynamics, and domains. In G. N. Reeke, R. R. Poznanski, K. A. Lindsay, J. R. Rosenberg, & O. Sporns (Eds.), Modeling in the neurosciences (pp. 339–371). Boca Raton, FL: CRC Press. Soderling, T. R. (1999). The Ca-calmodulin-dependent protein kinase cascade. Trends Biochem. Sci., 24(6), 232–236. Tanskanen, A. J., Greenstein, J. L., Chen, A., Sun, S. X., & Winslow, R. L. (2007). Protein geometry and placement in the cardiac dyad influence macroscopic properties of calcium-induced calcium release. Biophys. J., 92(10), 3379–3396. von Wegner, F., Wieder, N., & Fink, R. H. A. (2014). Microdomain calcium fluctuations as a colored noise process. Front. Genet., 5, 376. Weinberg, S. H. (2015). Membrane capacitive memory alters spiking in neurons described by the fractional-order Hodgkin-Huxley model. PLoS One, 10(5), e0126629. Weinberg, S. H., & Smith, G. D. (2012). Discrete-state stochastic models of calciumregulated calcium influx and subspace dynamics are not well-approximated by ODEs that neglect concentration fluctuations. Computational and Mathematical Methods in Medicine, 2012, 1–17. Weinberg, S. H., & Smith, G. D. (2014). The influence of Ca2+ buffers on free [Ca2+ ] fluctuations and the effective volume of Ca2+ microdomains. Biophys. J., 106(12), 2693–2709. Wieder, N., Fink, R. H. A., & von Wegner, F. (2011). Exact and approximate stochastic simulation of intracellular calcium dynamics. Journal of Biomedicine and Biotechnology, 2011(5231), 1–5. Wieder, N., Fink, R. H. A., & von Wegner, F. (2015). Exact stochastic simulation of a calcium microdomain reveals the impact of Ca2+ fluctuations on IP3R gating. Biophysical Journal, 108(3), 557–567. Williams, G. S. B., Chikando, A. C., Tuan, H.-T.M., Sobie, E. A., Lederer, W. J., & Jafri, M. S. (2011). Dynamics of calcium sparks and calcium leak in the heart. Biophys. J., 101(6), 1287–1296. Yang, L., MacLellan, W. R., Han, Z., Weiss, J. N., & Qu, Z. (2004). Multisite phosphorylation and network dynamics of cyclin-dependent kinase signaling in the eukaryotic cell cycle. Biophys. J., 86(6), 3432–3443.

Received July 14, 2015; accepted November 2, 2015.