## Abstract

Ca^{2+}-dependent signaling is often localized in spatially restricted
microdomains and may involve only 1 to 100 Ca^{2+} ions. Fluctuations in
the microdomain Ca^{2+} concentration (Ca^{2+}) can arise from a
wide range of elementary processes, including diffusion, Ca^{2+} influx,
and association/dissociation with Ca^{2+} binding proteins or buffers.
However, it is unclear to what extent these fluctuations alter
Ca^{2+}-dependent signaling. We construct Markov models of a general
Ca^{2+}-dependent signaling cascade and Ca^{2+}-triggered
synaptic vesicle release. We compare the hitting (release) time distribution and
statistics for models that account for [Ca^{2+}] fluctuations with the
corresponding models that neglect these fluctuations. In general, when
Ca^{2+} fluctuations are much faster than the characteristic time
for the signaling event, the hitting time distributions and statistics for the
models with and without Ca^{2+} fluctuation are similar. However, when
the timescale of Ca^{2+} 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 Ca^{2+} ions is small, a consequence of a long-tailed hitting time distribution. In a
model of Ca^{2+}-triggered synaptic vesicle release, we demonstrate the
conditions for which [Ca^{2+}] 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
Ca^{2+} fluctuations are an important aspect of microdomain
Ca^{2+} signaling and further suggesting that Ca^{2+} 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, Ca^{2+}-dependent signaling regulates important
physiologically subcellular processes. In neurons, Ca^{2+} signaling often
occurs in spatially localized subspaces, or Ca^{2+} microdomains, due to
spatially confined or restricted Ca^{2+} 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, Ca^{2+} signaling often
occurs over a wide range of timescales, on the millisecond timescale for SVR (Smith, 2004), to many hours, for example, for
Ca^{2+}-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 Ca^{2+} microdomains, the relevant
system size is often small, such that the likelihood of zero Ca^{2+} ions is
nonnegligible and significant: intracellular [Ca^{2+}] 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 Ca^{2+} ions.

Many previous studies have focused on the influence of stochasticity in
Ca^{2+}-driven processes, such as gating of Ca^{2+}-activated
channels (Williams et al., 2011; Gaur
& Rudy, 2011; Cannell, Kong, Imtiaz,
& Laver, 2013). Recently, several
modeling studies have demonstrated that [Ca^{2+}] fluctuations may
significantly influence the properties of Ca^{2+}-dependent processes,
specifically Ca^{2+}-regulated Ca^{2+} channels (Wieder, Fink, &
von Wegner, 2011, 2015). Flegg, Rüdiger, and Erban (2013) found that [Ca^{2+}] fluctuations reduce
the interpuff time in a cluster of Ca^{2+}-activated Ca^{2+} channels. Tanskanen, Greenstein, Chen, Sun, and Winslow (2007) showed that a continuous model of [Ca^{2+}]
overestimates channel open probability in cardiac myocytes, compared with a model
that accounts for individual Ca^{2+} ions and protein geometry. Weinberg and
Smith (2012) similarly showed in a minimal
Ca^{2+} microdomain model that [Ca^{2+}] fluctuations typically
reduced steady-state channel open probability in Ca^{2+}-activated channels.
Hake and Lines (2008) argued that a
random-walk model of Ca^{2+} signaling in the cardiac dyad can be well
approximated by a deterministic representation of Ca^{2+} 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 [Ca^{2+}] fluctuations over a longer
timescale is a general feature of Ca^{2+}-dependent signaling. More
specifically, in neurons, do presynaptic terminal [Ca^{2+}] fluctuations
alter release timing in a model of Ca^{2+}-triggered SVR? If so, at what
timescales does this occur, and how does this depend on the properties of the
Ca^{2+}-dependent signaling process?

During peak Ca^{2+} entry in the presynaptic terminal, SVR probability is
expected to be maximal, while when free [Ca^{2+}] is low and highly
buffered, SVR probability is expected to be low. However, Weinberg and Smith (2014) recently showed that Ca^{2+} binding proteins or buffers, which may be expected to reduce [Ca^{2+}]
fluctuations, in fact did not reduce and may enhance the magnitude of
[Ca^{2+}] fluctuations. Further, this enhancement of [Ca^{2+}]
fluctuations increased as average [Ca^{2+}] increased. Thus, it is not clear
to what extent [Ca^{2+}] fluctuations may influence SVR timing and how this
may depend on Ca^{2+} microdomain and SVR properties.

In our prior work, we showed that the relative size (coefficient of variation) of
[Ca^{2+}] fluctuations, driven solely by passive exchange (diffusion)
(Schuss, Singer, & Holcman, 2007), is
inversely proportional with the square root of the average number of microdomain
Ca^{2+} ions—that is, , where *c* is the average [Ca^{2+}] and *V* is the microdomain volume—and, as noted above,
Ca^{2+} buffers may in fact enhance the size of fluctuations, as the
Ca^{2+}-buffer association and dissociation reactions are an additional
source of intrinsic noise (Weinberg & Smith, 2012, 2014). Further, the
timescale of microdomain [Ca^{2+}] 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 [Ca^{2+}], microdomain volume, and [Ca^{2+}]
fluctuations timescale influence the dynamics of Ca^{2+}-dependent
signaling, in particular, multistep signaling cascades. Do [Ca^{2+}]
fluctuations alter Ca^{2+}-dependent signaling when fluctuations are much
faster, slower, or the same timescale as the Ca^{2+}-dependent signaling
event? Answering this question is complicated by the fact that the
Ca^{2+}-dependent signaling event itself is a source of intrinsic
[Ca^{2+}] fluctuations. Since [Ca^{2+}] 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 [Ca^{2+}]
fluctuations under- or overestimate the temporal dynamics of
Ca^{2+}-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 Ca^{2+} 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
[Ca^{2+}] fluctuations timescale and microdomain volume may influence
the timing of Ca^{2+}-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 Ca^{2+} binding steps preceding a final
Ca^{2+}-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 [Ca^{2+}] fluctuations on
Ca^{2+}-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

^{2+}-dependent signaling cascade, with or without [Ca

^{2+}] fluctuations present, that is, a model that either accounts for or neglects [Ca

^{2+}] 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 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 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).

*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), where is a matrix exponential and is the standard Heaviside function. The expected (mean) hitting time is given by and in general, the

*q*th moment of is given by where is a shorthand for the matrix inverse, raised to the

*q*th 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 Ca

^{2+}-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 [Ca

^{2+}]. 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 [Ca^{2+}] Fluctuations

*S _{n}* 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
[Ca

^{2+}] fluctuations are present.

#### 2.1.2 In the Presence of Microdomain [Ca^{2+}] Fluctuations

As discussed in section 1, the timescale
of [Ca^{2+}] fluctuations depends on many factors, such as
microdomain geometry and the presence and properties of Ca^{2+} buffers. Since the size and timescale of [Ca^{2+}] 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 [Ca^{2+}], microdomain volume, and timescale
separately. To investigate the influence of [Ca^{2+}] fluctuations
in general, without the introduction of several buffer- or other
noise-related parameters, we consider [Ca^{2+}] fluctuations that
arise as a consequence of passive exchange (diffusion) between the
microdomain and a bulk compartment, characterized by a single parameter . If Ca^{2+} passive exchange occurs with rate , with units of time^{−1}, then the
timescale of [Ca^{2+}] fluctuations , specifically the autocorrelation time, is simply given by (Weinberg & Smith, 2014).

^{2+}] fluctuations on cascade signaling, we consider a system of reactions including passive exchange (diffusion) of Ca

^{2+}(Schuss et al., 2007),

^{2+}-dependence, which arises due to a Ca

^{2+}-association reaction, that is, where

*k*is an association rate constant, with units of . For appropriate comparison with the hitting time calculation in the absence of [Ca

_{i}^{2+}] fluctuations, we define and . While Ca

^{2+}-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 [Ca

^{2+}].

^{2+}] fluctuations and Ca

^{2+}-dependent state transitions by , where ranges over all possible discrete values for the number of microdomain Ca

^{2+}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*: 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 Ca^{2+} 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 [Ca^{2+}] fluctuations present (see equation 2.6 and in equation 2.8) for an appropriate comparison. Since the number of
Ca^{2+} 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 Ca^{2+} ions is equal to . The timescales associated with [Ca^{2+}]
fluctuations and Ca^{2+}-dependent signaling cascades can range over
several orders of magnitude. Here, we focus on the influence of the relative
timescale between [Ca^{2+}] fluctuations and signaling; thus, we
define such that the passive exchange rate , where is a unitless [Ca^{2+}] fluctuations timescale
parameter (equivalently, we define ). Importantly, while defines the [Ca^{2+}] fluctuations timescale in a
microdomain in which passive exchange is the only elementary process
(Weinberg & Smith, 2014), here,
since we explicitly account for Ca^{2+} binding in the
Ca^{2+}-triggered signaling, the association and dissociation
reactions are an additional source of intrinsic noise that also drive
[Ca^{2+}] 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
[Ca^{2+}] fluctuations. We also calculate the (not normalized)
hitting variability *c _{v}* to illustrate how model parameters alter the hitting time
variability and the normalized hitting time variability, , to illustrate to what extent [Ca

^{2+}] 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 Ca^{2+} microdomain and does not consider spatial aspects of diffusion and
Ca^{2+} 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 Ca^{2+} signaling. Rüdiger, Falcke, and colleagues used a hybrid stochastic
spatial modeling approach, accounting for stochastic aspects of
Ca^{2+} 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 Ca^{2+} ion motion and showed that
neglecting [Ca^{2+}] fluctuations altered channel opening time.
Modchang et al. (2010)
investigated SVR in a stochastic spatial model and found that fluctuations
in Ca^{2+} channel gating were the most significant noise source
influencing SVR probability, but also that Ca^{2+} diffusion was a
significant noise source when peak [Ca^{2+}] and release probability
were low.

Although our study does not fully account for spatial aspects of
[Ca^{2+}] 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

*S*

_{2}molecule assuming one initial

*S*

_{1}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 [Ca^{2+}] 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 [Ca^{2+}] of (corresponding to an average of 6 Ca^{2+} ions), the hitting time distributions for
both fast and slow [Ca^{2+}] fluctuations also appear exponential
(see Figure 1A, left). When shown on a
logarithmic scale (right), for [Ca^{2+}] fluctuations much faster
than the timescale of the Ca^{2+}-dependent signaling reaction (, red), the hitting time distribution is nearly identical
to that as when [Ca^{2+}] fluctuations are absent. However, when
[Ca^{2+}] 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 Ca^{2+} ions), long hitting times are more
likely when [Ca^{2+}] 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 Ca^{2+} ions, zero Ca^{2+} ions
may be present in the microdomain, and because of the slow timescale of
passive exchange with the bulk, the time until one Ca^{2+} 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 [Ca^{2+}] fluctuations (which for this single-step cascade
are both equal to 1), as functions of the microdomain volume *V* and [Ca^{2+}] fluctuations timescale
parameter . For volumes of to 1 , and are biphasic functions of . When is either very small or large, and are near 1—[Ca^{2+}] fluctuations do not
alter the hitting time—which we explain as follows. When
[Ca^{2+}] fluctuations are much faster than the signaling
timescale (small ), even though the fluctuating [Ca^{2+}] is
typically above or below its expected value of , the duration of these events is very brief due to passive
exchange with the bulk. [Ca^{2+}] fluctuations are essentially
smoothed at the longer timescale, as proposed by Hake and Lines (2008), such that the timing of
Ca^{2+}-triggered signaling depends on only the average
[Ca^{2+}], . In contrast, when [Ca^{2+}] fluctuations are much
slower than the signaling timescale (large ), [Ca^{2+}] fluctuations are negligibly
influential, since [Ca^{2+}] is effectively constant on the
timescale of Ca^{2+}-triggered signaling, and again hitting time is
minimally altered.

However, for near 1—the timescale of [Ca^{2+}]
fluctuations and are similar— and are increased. In this case, when the expected number of
Ca^{2+} ions, , is small and fluctuations drive [Ca^{2+}] above
or below this value, the likelihood of zero microdomain Ca^{2+} ions
is nonnegligible, such that the duration of zero microdomain Ca^{2+} 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 Ca^{2+} ions increases, reducing
the likelihood of zero Ca^{2+} 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 Ca^{2+} 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 Ca^{2+} ions, and is largest for microdomain volumes near this value. In
Figure S2, we show that [Ca^{2+}] fluctuations driven by
Ca^{2+}-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
[Ca^{2+}] fluctuations for different values of and microdomain volume *V* (see Figure 3). As described above, simulations
demonstrate that when [Ca^{2+}] fluctuations are slow () and microdomain volume is larger (see Figure 3A), [Ca^{2+}] is effectively constant until the
transition from *S*_{1} to *S*_{2}. For a smaller microdomain volume (; see Figure 3B),
there are long time periods of zero microdomain Ca^{2+} 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 Ca^{2+} binding steps, each of which
decreases the number of Ca^{2+} ions, increase the likelihood of
zero Ca^{2+} 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 [Ca^{2+}] 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
[Ca^{2+}] fluctuations may alter Ca^{2+}-dependent
signaling. For larger microdomain volumes, we find similar qualitative
dependence on parameters, with and typically increased to lesser extents (not shown).

*V*, specifically , for the single-step signaling cascade, relating hitting time statistics to key model parameters (,

*V*, ), and , where 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

^{2+}] fluctuations in the general, multistep Ca

^{2+}-dependent signaling cascades. For simplicity, we consider an irreversible cascade, and for all

*i*, for which without [Ca

^{2+}] 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
[Ca^{2+}] 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), *c _{v}* (middle), and (right) as functions of and the cascade length

*n*. In general, we find that when [Ca

^{2+}] fluctuations are fast, and are near 1, as in the single-step cascade. However, when [Ca

^{2+}] fluctuations occur on a similar timescale or slower (i.e., ), increases, more so as

*n*increases, approaching an asymptotic value around

*n*= 5.

*c*decreases as

_{v}*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 [Ca

^{2+}] fluctuations, are a consequence of the increased likelihood of zero microdomain Ca

^{2+}ions, which results from multiple Ca

^{2+}binding steps in the signaling cascade.

*n*= 1. As before, , and as the number of cascade steps increases, 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 [Ca^{2+}]
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 [Ca^{2+}] fluctuations timescale
parameter , were highly dependent on these cascade properties and
with complex dependence.

### 3.2 Synaptic Vesicle Release

^{2+}] fluctuations alter timing in a general signaling cascade, we next analyze the influence of [Ca

^{2+}] fluctuations on SVR time. Experiments have shown that SVR occurs following sequential Ca

^{2+}binding steps preceding a final Ca

^{2+}-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 Ca

^{2+}-dependent SVR (Bollmann, 2000), where the release time is equivalent to the hitting time for the terminal step in the cascade: In this signaling cascade, the first five steps represent Ca

^{2+}sequentially and reversibly binding to five independent binding sites on Ca

^{2+}-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 Ca

^{2+}-bound sensor can then undergo a reversible Ca

^{2+}-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 [Ca

^{2+}],

*c*, we define the passive exchange rate relative to Ca

^{2+}binding at rest conditions and the [Ca

^{2+}] 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 [Ca

^{2+}] fluctuations. We also present analysis for slower [Ca

^{2+}]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 [Ca^{2+}] fluctuations
(dashed line) and fast [Ca^{2+}] fluctuations (red, magenta, green),
the distribution is approximately exponential (see Figure 6A). As [Ca^{2+}] 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 [Ca^{2+}] fluctuations timescale
parameter on the release time for varying levels of the average
[Ca^{2+}], *c*. For a given [Ca^{2+}]
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 [Ca^{2+}] values, as is .

We next consider the influence of varying several SVR model parameters. We
find that for all values of dissociation constant *K _{d}*, increases as increases (see Figure 7), as in Figure 6. Slow
[Ca

^{2+}] 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

*K*increases (decreased binding affinity). is increased to a larger extent for small

_{d}*K*, that is, [Ca

_{d}^{2+}] 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

*c*and , with increased up to a factor of 2 and over a similar range for

_{v}*c*values.

Interestingly, we find that increasing or decreasing the number of binding
sites has a similar effect as that of increasing or decreasing *K _{d}* (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 [Ca^{2+}], dissociation constant *K _{d}*, and [Ca

^{2+}] 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 [Ca

^{2+}] 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 [Ca

^{2+}] fluctuations are primarily driven by Ca

^{2+}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 [Ca

^{2+}] is effectively constant over the timescale of the signaling cascade and [Ca

^{2+}] 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, [Ca^{2+}]
fluctuations are relatively fast on the timescale of the signaling cascade,
and and are also not increased. However, for intermediate values
of *c* near *K _{d}*, the association and dissociation reactions are on a comparable
timescale, driving [Ca

^{2+}] fluctuations of comparable timescale as the signaling cascade; as a consequence, and are increased.

When *K _{d}* 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

*K*and small, the magnitude of [Ca

_{d}^{2+}] fluctuations is relatively large and increases and to a larger extent, compared with both large

*c*and

*K*and relatively small fluctuations. As decreases, [Ca

_{d}^{2+}] fluctuations are driven to a larger extent by passive exchange, which leads to faster [Ca

^{2+}] 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 [Ca^{2+}] fluctuations in the
setting of a time-varying Ca^{2+} influx. We simulate the
Hodgkin-Huxley neuron, augmented with a high-threshold Ca^{2+} channel current. A constant current stimulus is applied, such that the
neuron spikes repetitively at approximately 70 Hz, and the
Ca^{2+} current triggers Ca^{2+} influx and subsequent
Ca^{2+}-triggered SVR (see Figure 10A, top). To illustrate the significance of [Ca^{2+}]
fluctuations, we consider three model formulations: (1) the stochastic
SVR/Ca^{2+} model, in which both SVR and Ca^{2+} influx
are stochastic; (2) the stochastic SVR model, in which [Ca^{2+}] is
deterministic and governed by the Ca^{2+} current and mass-action
kinetics; and (3) the deterministic model, in which both SVR and
[Ca^{2+}] are governed by mass-action kinetics. Further details
of the model and simulations methods are provided in the appendix.

In Figure 10A, [Ca^{2+}] is
shown as a function of time for different values of . For small , Ca^{2+} transients are shorter due to faster
passive exchange, such that SVR events (solid lines below each trace)
typically occur during neuron spikes. For larger , Ca^{2+} transients are longer, such that SVR
events also occur between spikes, although at reduced frequency. Insets
illustrate [Ca^{2+}] 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 [Ca^{2+}] 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 Ca^{2+} influx (model 2, gray), the SVR time
distribution is left-shifted relative to model 1 (i.e., neglecting
[Ca^{2+}] 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 [Ca^{2+}]
fluctuations (see Figure 10C, bottom
left, gray) and overestimated when neglecting both SVR and [Ca^{2+}]
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 [Ca^{2+}]
fluctuations. For most values of , SVR time CV is also underestimated when neglecting
[Ca^{2+}] fluctuations (see Figure 10C, bottom right), demonstrating that
[Ca^{2+}] 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 [Ca^{2+}] fluctuations can greatly
influence the hitting time for a general Ca^{2+}-dependent signaling
cascade and specifically in a model of Ca^{2+}-triggered SVR. In a
general signaling cascade, small microdomain volume *V* and
multiple Ca^{2+} binding steps can lead to an increased likelihood of
zero microdomain Ca^{2+} ions, such that [Ca^{2+}] fluctuations
on the same timescale or slower than a Ca^{2+}-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 [Ca^{2+}] is
sufficiently large, such that the likelihood of zero Ca^{2+} ions is
small, slow [Ca^{2+}] fluctuations result in an effectively constant
[Ca^{2+}] and minimal influence on hitting times. Similarly, when
[Ca^{2+}] fluctuations due to passive exchange are faster (small ), even when *V* is small, [Ca^{2+}]
fluctuations also negligibly alter hitting times, since the
Ca^{2+}-dependent signaling effectively responds to the average
[Ca^{2+}].

In an SVR model, [Ca^{2+}] fluctuations may be driven by both passive
exchange and Ca^{2+} 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 Ca^{2+} binding
affinity *K _{d}* (see Figure 7) and the number of
Ca

^{2+}binding sites (see Figure 8) but negligibly depends on the number of vesicles available for release. Monte Carlo simulations similarly demonstrate that [Ca

^{2+}] 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
Ca^{2+}-dependent reactions are significant intrinsic sources of
[Ca^{2+}] fluctuations that can alter the temporal dynamics of
Ca^{2+}-dependent signaling cascades, including SVR and, further,
that a physiological model that neglects the influence of [Ca^{2+}]
fluctuations may be dramatically underestimating the timing of these
Ca^{2+}-mediated signaling events.

### 4.2 Comparison with Prior Simulation Studies

In this study, we demonstrate in both general and physiological settings that
[Ca^{2+}] fluctuations may or may not influence the timing of
Ca^{2+}-dependent signaling. In many parameter regimes, our study
agrees with the findings of Hake and Lines (2008) that [Ca^{2+}] fluctuations are smoothed by signaling
events occurring on a longer timescale; when [Ca^{2+}] 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
[Ca^{2+}] 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
[Ca^{2+}] fluctuations to greatly influence
Ca^{2+}-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
Ca^{2+} microdomain as a spatially restricted (well-mixed)
subcellular compartment.

In general, spatially localized Ca^{2+} 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 [Ca^{2+}] fluctuations are influential when accounting for the
complex geometry of the cardiac dyad, in agreement with our conclusion that
[Ca^{2+}] fluctuations are significant in spatially restricted
microdomains and specifically the presynaptic terminal. Collectively, these
studies and our results suggest that accounting for both [Ca^{2+}]
fluctuations and microdomain geometry is important.

### 4.3 Physiological Significance

In this study, we demonstrate that [Ca^{2+}] fluctuations can
significantly alter the timing of SVR and further that models of SVR neglecting
the influence of [Ca^{2+}] 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 [Ca^{2+}]
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 [Ca^{2+}] is
slightly elevated above resting levels, near 1 to 10 . Dynamics at such moderately elevated Ca^{2+} levels
may be significant physiologically, as Ca^{2+}-dependent signaling may
trigger feedback mechanisms that either suppress or further increase
Ca^{2+} influx (e.g., Ca^{2+}-activated or
Ca^{2+}-inactivated Ca^{2+} channels).

In many physiological settings, Ca^{2+} binding proteins or buffers can
play a significant role in the spatiotemporal dynamics of Ca^{2+} signaling. For example, buffers can alter the decay of residual Ca^{2+} following Ca^{2+} channel closure and influence release dynamics of a
Ca^{2+} release cluster (Dargan, Schwaller, & Parker, 2004). Buffers may play the additional role
of modulating the timescale of [Ca^{2+}] fluctuations (Weinberg &
Smith, 2014), and in some settings,
driving fast [Ca^{2+}] fluctuations (i.e., small ) and reducing the timing variability in
Ca^{2+}-dependent signaling. Of course, the influence of buffers may be
complex, as Ca^{2+} binding proteins such as calmodulin can both trigger
or modulate Ca^{2+}-dependent processes and drive [Ca^{2+}]
fluctuations.

While we investigate Ca^{2+}-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 [Ca^{2+}] fluctuation requires
many more iterations of the Gillespie stochastic simulation algorithm (SSA),
compared with slow [Ca^{2+}] 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
[Ca^{2+}] 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 [Ca^{2+}] fluctuations is influenced by diffusion,
microdomain geometry, and buffer properties. In this study, we demonstrated that
both passive exchange and Ca^{2+}-dependent association and dissociation
are potentially significant sources of [Ca^{2+}] fluctuations that can
alter hitting times, and we focus on the timescale of [Ca^{2+}]
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 [Ca^{2+}] fluctuations (i.e.,
buffers), competing Ca^{2+}-binding reactions, and other intrinsic or
extrinsic sources, occurring on potentially multiple timescales, are significant
in the context of Ca^{2+}-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

*V*for the single-step signaling cascade. As

*V*decreases, the expected number of microdomain Ca

^{2+}ions, , also decreases, such that for sufficiently small

*V*, the likelihood of more than one Ca

^{2+}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 , and . We then use equations 2.4 and 2.5 to analytically derive equations 3.2 and 3.3, respectively.

#### A.2 Model Details and Monte Carlo Simulation Methods for Synaptic Vesicle Release in the Spiking Neuron

^{2+}-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

*I*, given by the following system of differential equations: where , and the ionic currents are given by , , , and The equations governing the dynamics of the gating variables and

_{Ca}*n*are given by where

*v*is the shifted transmembrane potential, in which the resting potential mV has been subtracted. The instantaneous activation gating of

*I*is given by where and mV. Other parameter values are given as follows: maximum ionic current conductances (): , , ; and reversal potentials, in which

_{Ca}*V*has been subtracted (mV): , , , . In all simulations, . Since increasing [Ca

_{rest}^{2+}] fluctuations parameter decreases the passive exchange rate, for a given Ca

^{2+}influx, via

*I*, increasing results in a larger [Ca

_{Ca}^{2+}]. Therefore, the value of

*g*is decreased as is increased, such that the peak [Ca

_{Ca}^{2+}] in the deterministic model is fixed and approximately 70 . For , and , and 0.48 , respectively.

^{2+}model, we augment the set of reactions given in equation 2.7, with the zeroth-order reaction, where is the Ca

^{2+}flux, in units of (ions time

^{−1}), that can be calculated by multiplying

*I*by the microdomain surface area

_{Ca}*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

^{2+}signal

^{2+}release by stochastic RyR gating within a 3D model of the cardiac dyad and importance of induction decay for CICR termination

^{2+}signals in Xenopus oocytes by Ca

^{2+}-binding proteins

^{2+}ions in the dyadic cleft: Continuous versus random walk description of diffusion

^{2+}buffers on free [Ca

^{2+}] fluctuations and the effective volume of Ca

^{2+}microdomains

^{2+}fluctuations on IP3R gating