Abstract

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 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 (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, Rüdiger, 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 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, , 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 to  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 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 discrete-state continuous-time Markov chain (CTMC). If we define to denote the state of the CTMC that takes integer values between 1 and (the total number of states), then the probability of a transition from state i to state j during a small time interval is given by
formula
2.1
where matrix is known as the infinitesimal generator matrix. Conservation of probability implies that the diagonal entries of are given by . For a CTMC with a single absorbing state, , can be written in the form
formula
2.2
where is an matrix of transition rates between transient states, collects the rates into the absorbing state, and and are a column 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 , 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 , where is a row vector denoting the initial probability of state i such that . The probability density function for , , that is, the hitting time distribution, is given by Breuer and Baum (2005), Bladt (2005), and Lamar et al. (2011),
formula
2.3
where is a matrix exponential and is the standard Heaviside function. The expected (mean) hitting time is given by
formula
2.4
and in general, the qth moment of is given by
formula
2.5
where 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, , where the hitting time variance . 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 (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 consists of n reversible steps, with transition rates , with units of time−1, given by
formula
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 . If Ca2+ passive exchange occurs with rate , with units of time−1, then the timescale of [Ca2+] fluctuations , specifically the autocorrelation time, is simply given by (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),
formula
2.7a
where is the average [ 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,
formula
2.7b
formula
2.7c
formula
2.7d
where ki is an association rate constant, with units of . For appropriate comparison with the hitting time calculation in the absence of [Ca2+] fluctuations, we define and . 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 , where ranges over all possible discrete values for the number of microdomain Ca2+ ions, 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:
formula
2.8
where , , , and . Note that the transition rate between cascade stages is inversely proportional to the microdomain volume V, while the influx rate due to passive exchange 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 , equations 2.7b to 2.7d. Each ordered pair represents a state of two-dimensional CTMC. Since we are interested in the time until the system is absorbed into system state , 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 , where is a maximum number of Ca2+ ions. This value can be determined by incrementally increasing until there are insignificant changes in numerical calculations in hitting distribution and moments. We found that a reasonable value for , where is the ceiling function, yielding a matrix of dimension .

After enumerating the transient system states, such that state 1 corresponds to (0,1), 2 to (1,1),  , and the absorbing state corresponds to , we can write the infinitesimal generator matrix 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 in equation 2.8) for an appropriate comparison. Since the number of Ca2+ ions must be an integer and in general is not an integer, the initial probability is distributed appropriately between states and , such that the initial expected number of Ca2+ ions is equal to . 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 such that the passive exchange rate , where is a unitless [Ca2+] fluctuations timescale parameter (equivalently, we define ). Importantly, while 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, , to illustrate to what extent [Ca2+] fluctuations alter hitting time variability.

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 and (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. Rüdiger, 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, Rüdiger, Warnecke, & Falcke, 2012; Rückl et al., 2015). Flegg, Rüdiger, 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:
formula
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, , with an expected hitting time . The hitting time variance , and thus the relative hitting time variability .

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 and a resting [Ca2+] of (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 signaling reaction (, red), the hitting time distribution is nearly identical to that as when [Ca2+] fluctuations are absent. However, when [Ca2+] fluctuations are much slower than (, 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 are on the same timescale (, green).

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 (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) , .

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 (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) , .

In a smaller microdomain of volume, (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 . 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 reaction may be relatively long.

In Figure 2, we plot the hitting time mean (see equation 2.4) and coefficient of variation , normalized by the mean and CV, respectively, of 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 . For volumes of to 1 , and are biphasic functions of . When is either very small or large, and 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 ), even though the fluctuating [Ca2+] is typically above or below its expected value of , 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+], . In contrast, when [Ca2+] fluctuations are much slower than the signaling timescale (large ), [Ca2+] fluctuations are negligibly influential, since [Ca2+] is effectively constant on the timescale of Ca2+-triggered signaling, and again hitting time is minimally altered.

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 for different microdomain volumes V. (B) Hitting time mean and CV are shown as a function of V for different values of . 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: .

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 for different microdomain volumes V. (B) Hitting time mean and CV are shown as a function of V for different values of . 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: .

However, for near 1—the timescale of [Ca2+] fluctuations and are similar— and are increased. In this case, when the expected number of Ca2+ ions, , 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 are increased to a lesser extent.

When microdomain volume is smaller ( to ), both and increase as 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 are shown as a function volume V. increases several orders of magnitude when , which corresponds to an average of 3 Ca2+ ions, and 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 and microdomain volume V (see Figure 3). As described above, simulations demonstrate that when [Ca2+] fluctuations are slow () and microdomain volume is larger (see Figure 3A), [Ca2+] is effectively constant until the transition from S1 to S2. For a smaller microdomain volume (; see Figure 3B), there are long time periods of zero microdomain Ca2+ ions, resulting in an increase in hitting time.

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 . Simulations terminate when the system transition from state 1 to state 2 (i.e., the reaction ) occurs. Parameters: Microdomain volume V = 0.1 (A) and 0.01 (B) , .

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 . Simulations terminate when the system transition from state 1 to state 2 (i.e., the reaction ) occurs. Parameters: Microdomain volume V = 0.1 (A) and 0.01 (B) , .

The parameter regime for which the dependence of and on transitions from biphasic to monotonically increasing is on the order of to . 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 (). Thus, we find that when the number of cascade steps , and are monotonically increasing functions of , and slow [Ca2+] fluctuations increase and (see Figure S3). For the remainder of the study, we consider microdomain volumes at the small end of this range, , 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 typically increased to lesser extents (not shown).

In the appendix, we derive analytical approximations for the hitting time mean and variability, and , respectively, in the limit of a small microdomain volume V, specifically , for the single-step signaling cascade, relating hitting time statistics to key model parameters (, V, ),
formula
3.2
and , where
formula
3.3
From equation 3.2, it is clear that in this limit, , that is, the expected hitting time is always increased. Further, since , also always increases as increases.

We plot and in Figure 2 (thin dashed black lines, bottom of panel A and panel B) for values of V for which , 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 . Further, it is straightforward to show that in agreement with numerical calculations. In this small microdomain volume limit, we also find that , V, , and are related by .

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,
formula
3.4
and for all i, for which without [Ca2+] fluctuations, the hitting time distribution is well known, the hypoexponential distribution, (Breuer & Baum, 2005), with mean and variance . When transition rates between all states are equal, that is, , with mean , , and . 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 step signaling cascade, for different values of . In the absence of fluctuations (dashed line), as noted above, the hitting time distribution is hypoexponential. When [Ca2+] fluctuations are fast (small ), the hitting time distribution is similar to the distribution for the model without fluctuations. As increases, the distribution becomes right-shifted with a long tail, leading to an increase in the normalized expected hitting time .

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 (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: , , , n = 4.

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 (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: , , , n = 4.

In Figures 5A and 5B, we plot (left), cv (middle), and (right) as functions of and the cascade length n. In general, we find that when [Ca2+] fluctuations are fast, and are near 1, as in the single-step cascade. However, when [Ca2+] fluctuations occur on a similar timescale or slower (i.e., ), increases, more so as n increases, approaching an asymptotic value around n = 5. cv decreases as n decreases, while is increased for large 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 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.

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 are shown as functions of for different cascade lengths n. (B) , cv, and are shown as a function of n for different values of . Note that the black line in the top row () 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: , , .

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 are shown as functions of for different cascade lengths n. (B) , cv, and are shown as a function of n for different values of . Note that the black line in the top row () 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: , , .

For the general, irreversible cascade, we also derive the small microdomain volume approximation, shown in the appendix:
formula
3.5
and , where
formula
3.6
and equations 3.2 and 3.3 are equivalent to equations 3.5 and 3.6 for n = 1. As before, , and as the number of cascade steps increases,
formula
3.7
which approaches 1 for both large and small . For intermediate near 1, this limit is less than 1, although is still near 1 for values of . We plot and in Figure 5B (thin dashed black lines), and as in the single-step cascade, we find nice agreement with numerical calculations, with closer agreement for larger .

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 , 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 of the general signaling cascade studied above. Here, we consider a well-characterized model of Ca2+-dependent SVR (Bollmann, 2000), where the release time is equivalent to the hitting time for the terminal step in the cascade:
formula
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 ) and dissociation rate b (with units of ), related by dissociation constant (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 ), respectively, and a terminal step that represents the triggering of synaptic vesicle release from the activated state, with rate , 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 , such that . Note that physiological values of on the order of ms (Weinberg & Smith, 2014; von Wegner et al., 2014) correspond to values on the order of . As in the previous sections, we present analysis for ranging from because there are negligible differences between measures for less than , that is, faster [Ca2+] fluctuations. We also present analysis for slower [Ca2+]fluctuations, that is, values an order of magnitude greater than the upper limit of this range to , as these values may be appropriate in pathological settings that slow or impair diffusion.

3.2.1  Numerical Calculations of Vesicle Release Timing for Constant Average

Using the approach presented above, we plot the distributions for the release timing for different values of . 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 increases (cyan, blue), the distributions become more long tailed and right shifted, and increases, as in the previous general multistep example.

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 (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 are shown as functions of for different values of average [Ca2+], c. (C) , cv, and are shown as functions of c for different values of . Parameters (Bollmann, 2000): , , , , , , , , . In panel A, .

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 (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 are shown as functions of for different values of average [Ca2+], c. (C) , cv, and are shown as functions of c for different values of . Parameters (Bollmann, 2000): , , , , , , , , . In panel A, .

We next analyze the influence of the [Ca2+] fluctuations timescale parameter on the release time for varying levels of the average [Ca2+], c. For a given [Ca2+] level, monotonically increases as increases, and for some values of c, approaches an asymptotic value, while is only slightly increased for some parameters (see Figure 6B). For a given value of , is a biphasic function of c, such that for c near resting levels and very elevated near 100 , is not increased and near 1 (see Figure 6C), while for intermediate concentrations, between approximately 1 and 10 , is greater than 1 and generally increases as increases. 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 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 , with increased up to a factor of 2 and over a similar range for c values.

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 and Ca2+ dissociation constant Kd. Note that the plots in the third row () are the same as in Figure 6C. Parameters: , . Other parameters as in Figure 6.

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 and Ca2+ dissociation constant Kd. Note that the plots in the third row () are the same as in Figure 6C. Parameters: , . Other parameters as in Figure 6.

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
formula
3.9

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 , respectively.

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 and number of Ca2+ binding sites . Note that the plots in the fourth row () are the same as in Figure 6C. Parameters as in Figure 6.

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 and number of Ca2+ binding sites . Note that the plots in the fourth row () are the same as in Figure 6C. Parameters as in Figure 6.

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 differences of over the entire parameter space for 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 when is large (see Figure 9). Although the passive exchange rate is constant in these simulations, varying c clearly alters the timescale of [Ca2+] fluctuations, such that fluctuations are faster for smaller c. For large , both the association and dissociation rate, and b, respectively, are much faster than passive exchange rate , 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 , such that and the association rate is faster than dissociation, , 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 [ is effectively constant, and in general are not increased.

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: . Other parameters as in Figure 6.

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: . Other parameters as in Figure 6.

In contrast, for , such that the dissociation rate is faster, SVR occurs relatively slowly, [Ca2+] fluctuations are relatively fast on the timescale of the signaling cascade, and and 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 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 to a larger extent, compared with both large c and Kd and relatively small fluctuations. As decreases, [Ca2+] fluctuations are driven to a larger extent by passive exchange, which leads to faster [Ca2+] fluctuations, which in turn results in and 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 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.

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 , 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 : model 1 (color), 2 (gray), and 3 (black). (C) Top: SVR time mean (left) and CV (right) are shown as functions of 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 . SVR parameters as in Figure 6.

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 , 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 : model 1 (color), 2 (gray), and 3 (black). (C) Top: SVR time mean (left) and CV (right) are shown as functions of 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 . SVR parameters as in Figure 6.

In Figure 10A, [Ca2+] is shown as a function of time for different values of . For small , Ca2+ transients are shorter due to faster passive exchange, such that SVR events (solid lines below each trace) typically occur during neuron spikes. For larger , 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 . For small (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 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 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 in Figure 10C. For all three models, the expected SVR time decreases as increases. For all values of , 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 , largest typically near for model 1, and for most values of is largest when accounting for [Ca2+] fluctuations. For most values of , 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 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 ), 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 [ 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 ), 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 . 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.

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 is much longer than the characteristic time associated with diffusion across (within) the domain . Schuss et al. (2007) showed that the timescale for microdomain escape via a small opening (“escape window”) of radius r, , where D is the diffusion coefficient. For three-dimensional diffusion, , where L is a length characteristic of the microdomain and if we assume , then the ratio of these timescales 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 . 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 ) and 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 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 , then , and the size of the state space (i.e., the dimension of the matrix) for a signaling cascade of n = 10 steps, , is large, such that is a (sparse) matrix with approximately elements. For an exceedingly large state space, numerical calculations will be infeasible, and Monte Carlo simulations will be necessary to analyze hitting time statistics.

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 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, , also decreases, such that for sufficiently small V, the likelihood of more than one Ca2+ ion present in the domain becomes negligible: . Therefore, we can define the maximum value for C in the CTMC described by equation 2.8, , which simplifies considerably. For the case of the single irreversible reaction ,
formula
A.1
and . We then use equations 2.4 and 2.5 to analytically derive equations 3.2 and 3.3, respectively.
For the general, irreversible cascade , we can also derive the small microdomain volume approximation. Again assuming , we can define , simplifying , such that the matrix, , is given by
formula
A.2
, 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:
formula
A.3a
formula
A.3b
where , and the ionic currents are given by , , , and The equations governing the dynamics of the gating variables and n are given by
formula
where v is the shifted transmembrane potential, in which the resting potential mV has been subtracted. The instantaneous activation gating of ICa is given by
formula
where and mV. Other parameter values are given as follows: maximum ionic current conductances (): , , ; and reversal potentials, in which Vrest has been subtracted (mV): , , , . In all simulations, . Since increasing [Ca2+] fluctuations parameter decreases the passive exchange rate, for a given Ca2+ influx, via ICa, increasing results in a larger [Ca2+]. Therefore, the value of gCa is decreased as is increased, such that the peak [Ca2+] in the deterministic model is fixed and approximately 70 . For , and , and 0.48 , 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,
formula
A.5
where 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 , where radius is calculated from the volume V. is determined by numerical integration of equation A.3. In the stochastic SVR model, 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 .

Acknowledgments

I acknowledge and thank Old Dominion University High Performance Computing for use of the Turing cluster to perform numerical calculations and simulations.

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 IP3-mediated 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
.
Flegg
,
M. B.
,
Rüdiger
,
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.
,
Rüdiger
,
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
.
Neher
,
E.
, &
Sakaba
,
T.
(
2008
).
Multiple roles of calcium ions in the regulation of neurotransmitter release
.
Neuron Review
,
59
(
6
),
861
872
.
Rückl
,
M.
,
Parker
,
I.
,
Marchant
,
J. S.
,
Nagaiah
,
C.
,
Johenning
,
F. W.
, &
Rüdiger
,
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 calcium-regulated 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
.

Supplementary data