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.
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.1 Calculation of Hitting Times in Signal Cascades
2.1.1 In the Absence of Microdomain [Ca2+] Fluctuations
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).
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.1 General Signaling Cascade
3.1.1 Single Irreversible Reaction
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).
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.
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.
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).
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
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 .
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.
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
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.
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.
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.
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.
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.
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.
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
A.2 Model Details and Monte Carlo Simulation Methods for Synaptic Vesicle Release in the Spiking Neuron
I acknowledge and thank Old Dominion University High Performance Computing for use of the Turing cluster to perform numerical calculations and simulations.