Models of stochastic Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ca}^{2+}$$\end{document} spiking

Complexity and limited knowledge render it impractical to write down the equations describing a cellular system completely. Cellular biophysics uses hypotheses-based modelling instead. How can we set up models with predictive power beyond the experimental examples used to develop them? The two textbook systems of cellular biophysics, Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ca}^{2+}$$\end{document} signalling and neuronal membrane potential dynamics, both face this question. Both systems also have a non-equilibrium feature in common: on different time scales and for different observables, they exhibit stochastic spiking, i.e., sequences of stereotypical events that are separated by statistically distributed intervals, the interspike intervals (ISI). Here we review recent progress on the description of Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ca}^{2+}$$\end{document} spikes in terms of blips, puffs and cellular Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ca}^{2+}$$\end{document} spikes and focus on stochastic models that can explain the statistics of the single ISIs, in particular its mean and variance and the cell-to-cell variability of these statistics. We also review models of the stochastic integrate-and-fire type and measures like the spike-train power spectrum or the serial correlation coefficient that are used to describe neuronal spike trains. These concepts from computational neuroscience might be applicable for understanding long-term memory effects in Ca2+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ca}^{2+}$$\end{document} spiking that extend beyond a single ISI, such as cumulative refractoriness.


Introduction
The program of theoretical physics for understanding a given system is to specify first principles to it and to solve the resulting equations. That program has been extremely successful and defined our idea of an exact and quantitative science. The predictive power of the first principles originates from the astonishing correspondence between experimental objects and mathematical structures. The mechanics of macroscopic objects corresponds to variational principles and differential equations, the behaviour of microscopic objects corresponds to operator theory in Hilbert spaces.
The biophysics and biochemistry of cells obey the first principles, too. But cells consist of many components and interactions. Specifying the fundamental equations of physics to a living cell is close to impracticable. The approach of theoretical biophysics is consequently what usually is called mathematical modelling. Instead of a derivation from first principles, a hypothesis on the components and interactions assumed to be most relevant for a specific process of interest defines the model equations. The assumptions need to be verified retrospectively by contrasting model predictions with experimental results. Modelling has to find the a e-mail: nicolai.friedhoff@mdc-berlin.de (corresponding author) balance between capturing all relevant components, manageable complexity and the purpose of the model. Within this balance and in particular since modelling lacks the certainty of first principles, it is fundamental to start the formulation of the model equations within the mathematical structures to which the system to be modelled corresponds to gain predictive power. Otherwise a model might capture the experiment used to develop it, but very likely fails in predictions beyond this specific setting.
Only a few cellular dynamical systems are currently characterized well enough for identifying the mathematical structure corresponding to them. Intracellular Ca 2+ dynamics is one of them. The Ca 2+ pathway translates extracellular signals into intracellular responses by increasing the cytosolic Ca 2+ concentration in a stimulus dependent pattern [7,32,94]. The concentration increase can be caused either by Ca 2+ entry from the extracellular medium through plasma membrane channels, or by Ca 2+ release from internal storage compartments. In the following, we will focus on inositol 1,4,5-trisphosphate (IP 3 )-induced Ca 2+ release from the endoplasmic reticulum (ER), which is the predominant Ca 2+ release mechanism in many cell types. IP 3 sensitizes Ca 2+ channels (IP 3 Rs) on the ER membrane for Ca 2+  its own release channel is called Ca 2+ -induced-Ca 2+release (CICR). Opening of an IP 3 R triggers a Ca 2+ flux into the cytosol due to the large concentration differences between the two compartments. CICR sometimes strongly multiplies channel opening to a global release and concentration spike. The released Ca 2+ is removed from the cytosol either by sarco-endoplasmic reticulum Ca 2+ ATPases (SERCAs) into the ER or by plasma membrane Ca 2+ ATPases into the extracellular space.
IP 3 Rs are spatially organized into clusters of up to about fifteen channels within an area with a diameter of 100-500 nm. These clusters are scattered across the ER membrane with distances of 1-7 µm [10,53,59,85,92,93]. CICR and Ca 2+ diffusion couple the state dynamics of the channels. Given that the diffusion length of free Ca 2+ is less than 2 µm due to the presence of Ca 2+ binding molecules in the cytoplasm and SERCAs, the coupling between channels in a cluster is much stronger than the coupling between adjacent clusters [96]. The structural hierarchy of IP 3 Rs from the single channel to clusters and cluster arrays on cell level shown in Fig. 1 is also reflected in the dynamic responses of the intracellular Ca 2+ concentration as revealed through fluorescence microscopy and simulations [10,62,97,108]. Openings of single IP 3 Rs (blips) may trigger collective openings of IP 3 Rs within a cluster (puffs). Ca 2+ diffusing from a puff site can then activate neighboring clusters, eventually leading to a global, i.e., cell wide, Ca 2+ spike [35,53,62,63]. Repetitive sequences of these Ca 2+ spikes encode information that is used to regulate many processes in various cell types [7,55,73].
Ca 2+ exerts also a negative feedback on the channel open probability, which acts on a slower time scale than the positive feedback, and has a higher Ca 2+ half maximum value than CICR [10,50,63,67,99,108]. This Ca 2+ -dependent negative feedback helps terminating puffs. Therefore, the puff probability immediately after a puff is smaller than the stationary value but typically not 0. Channel clusters recover within a few seconds to the stationary puff probability from this Ca 2+dependent inhibition [10,50,63,67,99,108].
The negative feedback terminating global release spikes causes an absolute refractory period T min as part of the interspike intervals (ISIs) lasting tens of seconds [71,100,107]. The molecular mechanism of this feedback is pathway and cell type specific and not always known. A negative feedback on the IP 3 concentration might be involved [5,69]. Hence, the negative feedback that determines the time scale of interspike intervals is different from the feedback contributing to interpuff intervals. It requires global (whole cell) release events.
Modelling of Ca 2+ signalling has relied heavily on ordinary differential equations in the last decades, established as the rate equations for the average fractions of IP 3 Rs in states corresponding to IP 3 R state schemes and spatially averaged Ca 2+ , IP 3 and buffer concentrations [86][87][88]104]. This approach neglects noise and fluctuations [89]. However, the experimental evidence both on puffs and sequences of global spikes An open channel entails Ca 2+ release into the cytosol due to the large concentration difference between the ER and the cytosol. Since channels are clustered, opening of a single channel, which is called a blip, leads to activation of other channels in the cluster, i.e., a puff (middle). The cluster corresponds to a region with Ca 2+ release with a radius R cl that is fixed by the number of open channels. The stochastic local events are orchestrated by diffusion and CICR into cell wide Ca 2+ waves, which form the spikes on cell level (top). (Figure from ref. [83].) demonstrated random behavior and, therefore, the relevance of higher moments. Additionally, most models do not distinguish between local and global processes and feedbacks. That entailed in the end dependencies of system characteristics like, e.g., the average interspike interval (period) on measurable parameter values which deviate from experimental observations or require parameter values not supported by measurements [88,104]. The purpose of most models is to simulate cellular behavior, and ordinary differential equations are very convenient to that end. Their derivation, however, has to take the large fluctuations into account, i.e., has to start from stochastic theory as the mathematical structure corresponding to Ca 2+ dynamics. We will illustrate with the Siekmann IP 3 R model, how this might be done.
An alternative to simulating cellular behavior by differential equations is to determine the distribution of cellular properties generated by the noise inherent to the system [38,54]. Such an approach would correspond more to the noisy character of cell dynamics, but will only take hold, if the analysis of experimental results engages into such a view on cellular behavior and measures distributions and/or their moments [54]. We will discuss a concept for calculating the first two moments of the interspike interval distribution.
Ca 2+ spikes and their statistical measures have also some similarity with sequences of neural action potentials, the famous neural spike trains. We will also briefly discuss how concepts from computational neuroscience, such as multidimensional integrate-and-fire models and spike train power spectra could be useful to model and analyze Ca 2+ spiking.  [8]. [Ca 2+ ] around or in clusters becomes very large fast during puffs, easily reaching tens of µM or more, Fig. 2 [6]. These are concentrations in the inhibitory regime (s. microdomains), such that Ca 2+ release also has a fast negative feedback component on clusters.
Recent experiments on puff behavior of all three isoforms of the IP 3 R shed light on their local dynamics in form of puff frequency, puff amplitudes, open channels per puff, rise and fall times, and duration [59]. The average puff duration (full duration at half-maximum, FDHM) is about 41 ms ± 3 ms for wild-type IP 3 Rs. While the opening of clusters is explained by CICR within clusters, possible closing mechanisms of single IP 3 Rs and clusters are still being discussed. Among possible puff termination mechanisms are stochastic attrition (there is always a probability for many channels to spontaneously close together within a short time window), local ER-depletion (the ER becomes devoid of Ca 2+ locally, not able to support local cluster Ca 2+ efflux), luminal activation (regulation by Ca 2+ or other molecule species on the ER-side of the IP 3 R), or coupled gating (single closing may trigger closing of the cluster due to coupled channel dynamics) [90]. High [Ca 2+ ] together with the biphasic Ca 2+ dependency is also assumed to be at least a contributing factor of puff termination [106].
Various single channel behaviours in the course of a puff have been measured in experiments. While a steep increase of the fluorescence signal measured directly at puff sites as a quick opening of coupled channels is common, the termination of puffs can be realized in numerous ways. Smooth decay, step-wise decay, or closing with infrequent re-opening or bursting re-openings are among the most occurring channel closing scenarios or puff shapes, respectively [106]. Sometimes multiple IP 3 Rs within one cluster close almost in nearsynchrony in experiment on some occasions, yielding the seldom occurring block puff [106]. This occurred more often compared to expectations based on sets of independently closing channels (stochastic attrition).
Observation of neighbouring open IP 3 Rs within clusters with either one or two open channels confirmed deviations from the behavior of pairs of independent channels. This overall behaviour cannot be explained by inhibitory fast high Ca 2+ (biphasic open probability) or local ER-depletion, suggesting an important but yet unknown channel coupling mechanism leading to coupled gating that renders puff duration and channelcoupled puff termination robust.
While regulation of IP 3 Rs by luminal Ca 2+ content or other molecules inside the ER has been a seemingly intractable question for decades, recent experimental studies have found further support for the hypothesis of luminal control. IP 3 Rs have been reported to be regulated by luminal [Ca 2+ ] ER and likely the widely-expressed luminal Ca 2+ buffer protein annexin A1 (ANXA1) which together inhibit IP 3 Rs at high [Ca 2+ ] ER [102].
New findings suggest that IP 3 Rs have two distinct modes of Ca 2+ release. A punctate liberation mode during the rise of the Ca 2+ transient which is then followed by a diffuse mode that sustains global Ca 2+ release. The punctate mode is terminated before reaching the peak, likely through an yet unknown mechanism regulated by [Ca 2+ ] ER . These two modes could also target different effector species, regulating different downstream elements of the IP 3 induced Ca 2+ signalling pathway. [60] 2.1.2 The dynamic regime of the local dynamics Intracellular Ca 2+ dynamics is a reaction diffusion system. The reactions comprise release of Ca 2+ from the ER, pumping by SERCAs, buffering and the binding/unbinding with other Ca 2+ binding sites. The reaction dynamics is local, diffusion provides the spatial coupling. The dynamic regime (excitable, bistable or oscillatory) of a reaction diffusion system is dominated by the the dynamic regime of the local dynamics. From a structural point of view, the local dynamics are the cluster dynamics.
[Ca 2+ ] profiles in the vicinity of single IP 3 Rs and within clusters cannot be measured directly, but can be simulated [96] or calculated analytically in good  While [Ca 2+ ] peaks at the cluster located at r=0 µm, [Ca 2+ ] distanced from the cluster will be one to two orders of magnitude smaller. Since IP3R dynamics are subject to [Ca 2+ ] in very close proximity to them, this makes meaningful cell wide spatial averaging difficult at best [6] approximation [6,96]. The [Ca 2+ ] at the cluster locations is about one or two orders of magnitude larger than spatially averaged concentration values, and decreases steeply with increasing distance from the channel, Fig. 2. This leads to the existence of microdomains of large [Ca 2+ ] at clusters with open channels, which are only weakly coupled to neighboring clusters by steep concentrations gradients. It is the local Ca 2+ dynamics that affects cluster dynamics the most.
The Ca 2+ concentration at closed clusters is the resting concentration in the range of ≤ 100 nM. Concentrations at open channels are >20 µM [6,96]. The dynamic range of the regulatory binding sites for both the positive and negative feedback of Ca 2+ to the open probability extends from a few hundred nM to micromolar values below 10 µM [43,49,95]. Oscillatory dynamics require concentration values in the dynamic range. However, with these large concentration changes, the system essentially never is in this dynamic range and the regime of the deterministic limit of the cluster dynamics is either excitable or bistable (except tiny parameter ranges) [97]. This conclusion is supported by an investigation into the time scales on cluster level. Typical interpuff intervals last a few seconds [25,26,53,99], interspike intervals are in the range from about 20 s to a few minutes. If the local dynamics were oscillatory and caused the sequence of spikes, the time scale of the ISI should be detectable as a temporal modulation of properties of the puff sequence at a given site. That has not been found [99]. A modulation of puff sequences on the ISI time scale could not be detected and no evidence of an oscillatory regime of the local dynamics has been observed [99]. The ISI time scale has only been observed on cell level.
Replacing local Ca 2+ concentrations with globally averaged [Ca 2+ ] values as the input for IP 3 Rs, even though their values differ by orders of magnitude, leads to misleading IP 3 R and Ca 2+ dynamics [97]. Averaged global concentrations during spikes are in the dynamic range of the IP 3 R regulatory binding sites thus allowing for cluster-cluster coupling. Using them in mathematical models as the Ca 2+ concentration experienced by the IP 3 R entailed oscillatory dynamics. However, that dynamic regime shrinks to negligible parameter ranges, high frequency and tiny global amplitudes with realistic local concentrations [97] and could not be verified by local measurements [99]. Thus IP 3 R Ca 2+ dissociation constants guarantee spatial coupling but do not allow oscillatory local dynamics.

Interspike intervals of global spikes are random
Once a cluster of IP 3 Rs opens to create a puff, the released Ca 2+ diffuses within the cell. If it reaches neighbouring clusters there is a probability of triggering follow up puffs. This can then become a self amplifying process, until a critical number of open clusters is reached, resulting in a cell-wide Ca 2+ release event, called a Ca 2+ spike [61]. These global spikes can be measured similar to measuring local puffs and can be described with the same quantities, like interspike interval (ISI), duration, or amplitude. Measuring a sequence of Ca 2+ spikes over a few minutes to hours yields a spike train from which we obtain the sequence of interspike intervals, Fig. 3. Just like blibs and puffs, spike times are inherently random, the ISI as a property of subsequent spike times is random as well. A global Ca 2+ spike has an inhibitory effect on subsequent puff events. The recovery form that inhibition takes tens of seconds, i.e., it is negative feedback on long time scales. It creates an absolute refractory period T min during which no puffs occur.
We can quantify how random spike timing of a given spike train is by the relation between the standard deviation σ of ISIs and the average ISI, T av . We see in Fig. 3 that they are linearly related like Such a linear relation has been found for all cases investigated (8 cell types and 10 conditions [17,28,31,81,100], see also [68]). The coefficient of variation of the stochastic part T av − T min of the ISIs is CV = σ/(T av − T min ) = α. The larger the CV, the more stochastic is the output of a given process. A CV value equal to 1 indicates a Poisson process, which is maximally random. A vanishing CV indicates a deterministic process. We determined the CV or α resp. as the slope of the linear approximation to population data as in Fig. 3, and from 2 experimental conditions with an individual cell. We found both values to agree [80,82] turning α into an observable not subject to cell variability (which is different from the results with puff sites in this respect [99]). Additionally, the value of α turned out to be robust against changes of buffering conditions [81], represents the data of one experiment, i.e., measured spike train of a cell. The wide spread indicates large cell-to-cell variability, but there is a functional σ-Tav moment relation visible as a linear fit. Plots from [81] stimulation strength and three pharmacological perturbations of the Ca 2+ signalling system. That surprising robustness turns Eq. (1) into one of the equations defining Ca 2+ signalling from the perspective of quantitative approaches. The value of α is set by the time scale of recovery from global negative feedback terminating the release spikes [98].

The relation between average interspike interval and stimulation
Cells are stimulated by extracellular agonists [A] binding to receptors in the cell membrane. The strength of stimulation controls the intracellular concentration of IP 3 . In general, we observe only puffs at low stimulation, spikes at intermediate agonist concentration and maintained high Ca 2+ concentration in some cell types and with some pathways at very strong stimulation. Within the spiking regime, cells respond to an increase of agonist concentration with a decrease of the average ISI, T av [36,100]. It was found for all pathways tested that the population averaged response could be well fit to a single exponential function which depends on the strength of the stimulus, given by the extracellular agonist concentration [A], Fig. 4, that is Here T min is the smallest ISI reached at strong stimulation, i.e., the absolute refractory period plus spike duration, and T ref the reference ISI at a reference agonist concentration [A ref ]. β is a constant for a given cell type and signalling pathway. We also found β to be the same for all individual cells. Hence, it is another observable defining Ca 2+ signalling from the perspec- That exponential dependency on stimulation in Eq. (2) follows from paired stimuli experiments, i.e., runs much deeper than a simple direct fit of an ansatz to experimental data. The change of the average stochastic part of the ISI ΔT av due to an agonist concentration step is proportional to the average stochastic part T av1 − T min at the lower agonist concentration T av1 [100]:

Long time scales from slow global processes and small spike probabilities
With some cell types, individual cells or experimental situations, T av is much longer than any time scale that is relevant for the state dynamics of clusters or even global cellular dynamics. From a dynamical systems point of view applying to deterministic models, this should not be possible, since each time scale requires a process setting it. However, long time scales might result simply from small probabilities and not from a slow process. Decay of a single radioactive atom for example happens at a random moment in time. If the atom is rather stable, decay is unlikely and it takes a long time on average to happen. But there is no process leading to the decay event. The state of the atom is stationary up to the time of the event. That may also apply to spike generation with the cell in the role of the atom and generation of a spike corresponding to the decay event. If the spike generation probability is small, we may observe long average ISIs and the state of the cell before the spikes is essentially stationary.
There is no process setting the long time scale in that case. Alternatively, there might be a slow process setting a long average ISI. The recovery from the negative feedback, which terminates spikes, is a prime candidate for such a slow process. The negative feedback might for example decrease [IP 3 ] [5], which then needs to recover before the next spike can occur. The inhibitory effect is a substantial decrease of the puff probability, which entails an absolute refractory period.
We can use the CV or α to assess the relative weight of small probability vs slow process in setting T av . If the CV is equal to 1, the ISIs follow an exponential distribution and are maximally random. There is no slow process setting the long time scale in that case, very similar to nuclear radioactive decay of an atom. A CV of 0 would indicate a purely deterministic and noisefree process with vanishing deviation. If the CV value is between 0 (deterministic) and 1 (pure randomness), a slow process changes the spike probability without rendering spike generation deterministic. Note, the average ISI is not simply the inverse of the recovery rate in that regime [39]. Measured CVs are between 0.2 (e.g., stimulated hepatocytes) and 0.98 (e.g., spontaneous spiking in microglia).

Open problems
We consider as open problems what is lacking for a theory able to derive the cellular signals from molecular properties. The large cell-to-cell variability defines here what is meaningful to be described by theory. An intuitive explanation for cell variability among many other possibilities might be the differences in the relative cluster positions. However, also this question has not been exhausted yet.
The puff property distributions for amplitude, duration and IPI have been simulated or described by ansatzes by a variety of groups [14][15][16], but have not been analytically derived yet. We cannot expect analytic expressions using realistic channel models (see below), but the distributions have not been written down even for strongly simplified models. Lock et al. recently demonstrated that all three IP 3 R isoforms generate similar puff property distributions sampled from many puff sites [59]. Hence, the distributions cannot depend on detailed molecular properties and a simplifying approach as common ground would make sense and would be a starting point providing conceptual understanding.
The situation with respect to global signals is similar. The interspike interval distribution for ISI sequences normalized by the average has been measured for HEK cells and spontaneously spiking astrocytes [84] and simulated [37], but it has not been derived yet. The robustness of the coefficient of variation CV has been very well confirmed experimentally [81,98,100] and has been simulated [72,82,83,98], but has neither been derived in some analytical work.
The concentration response curve of the average ISI shows an exponential dependency on the extracellular agonist concentration stimulating the cell [100]. The agonist sensitivity in the exponent is cell type and pathway specific [100]. The pre-factor of the exponent picks up all the cell variability. This detailed knowledge on the concentration response curve also awaits its theoretical explanation.
Open problems with respect to methods mainly concern the role of fluctuations. The large values of coefficients of variation on all structural levels demonstrates that fluctuation are not negligible. Their potential role becomes more tangible by considering intracellular Ca 2+ signalling as a deterministic reactiondiffusion system. The dynamic regime is then fixed by the local dynamics. We have no experimental evidence for an oscillatory local dynamics of intracellular Ca 2+ signalling [99], and the whole literature on puffs suggests the local dynamics to exhibit only time scales of a few seconds. The experimental results are compatible with an excitable regime of the local dynamics. Consequently, spikes are due to fluctuations. Concepts taking fluctuations along in systems of ordinary differential equations (ODEs) exist [109], but have not been applied to the system, yet. We will discuss them also below.

Modelling concepts from molecular properties to global dynamics including fluctuations and noise
The essence of the Ca 2+ signalling system is defined by its general properties, which are also the basic requirements models should meet: -The sequence of dynamic regimes with increasing stimulation: puffs, spikes, permanently elevated Ca 2+ . Pathway dependent also a bursting regime may follow or replace the spiking regime.  (2) and (3) with T min , α and γ being cell type and pathway specific but not subjected to cell variability. -ISIs depend sensitively on parameters of spatial coupling.
These general properties apply to all cells. Cells exhibit variability with respect to concentrations of the functional proteins, geometry of clusters and the cellwide cluster array, ER luminal Ca 2+ content etc. The general properties of Ca 2+ signalling cannot depend on the details of these highly variable cellular characteristics, which calls for models as simple as possible but meeting the above requirements.
Puff models should start from the molecular properties of the IP 3 R. Its random state changes are the source of noise. We will use one of the most recent models of the IP 3 R to describe concepts, the Siekmann model [79], which is a Markov model based on single-channel data. We will discuss in that context, how fluctuations might enter ODE-focused modelling approaches.
Puff property distributions form the basis for modelling of global dynamics. We will discuss concepts for calculating moments of ISI distributions. Most current models adapt molecular rate constants to global time scales to reproduce measured average ISI values. However, the origin of the time scales on global level are global processes. We will sketch how to introduce these global processes into the coupling between the puff dynamics and global dynamics which allows for using realistic molecular parameters.

IP 3 R clusters as ensembles of receptors described by the Siekmann model
Several channels (up to fifteen) form a cluster. The opening of one receptor channel within a cluster ('blib') increases the open probability of the other channels in the cluster due to strong channel coupling by Ca 2+ diffusion, which may cause a puff. We consider a cluster consisting of a stochastic ensemble of N channels.
We denote the number of channels in state i according to Fig. 5 as n i ≥ 0 with N = i n i , effectively removing one degree of freedom due to this require- Only the rates connecting these two modes are Ca 2+ dependent [40,79] ment. A puff occurs if some critical number of channels is in the open state, motivating to study the expectation value to be in a state i, n i .
The total change in probability for a set {n i } = {n 1 , n 2 , ..., n 6 } = {n 1 , n 2 , ..., n 5 , N − n 1 − · · · − n 5 } is given by the probability fluxes for each single channel transition from state j to state k, resulting in a change of {n i } like n j → n j − 1 ∩ n k → n k + 1.
We suggest to study whether higher moments may drive puff dynamics and where the hierarchy of moment equations can be cut off. This might lead to a set of ODEs as Ca 2+ signalling model with realistic parameter values, which establishes the ability to simulate time courses with all the computational comfort ODEs provide.

Spike generation as first passage process with time dependent transition probabilities
We would like to present a concept calculating the moments of the ISI distribution in this section as it naturally corresponds to the random spike timing. We also suggest a method to invoke global processes modulating the local dynamics.
The stochastic element of such a formulation of spike generation is a single cluster described by its IPI, puff duration and amplitude distributions. Such an approach dispenses with detailed intracluster concentration dynamics [98]. A model in the same spirit set up to simulate time courses has been developed by Calabrese et al. [14]. Clusters open sequentially. Once a critical number N cr of open clusters has been reached, the remaining ones will open with almost certainty due to coupling by Ca 2+ diffusion and the positive feedback by CICR. There are many (N paths ) paths from all clusters closed to this critical number (see Fig. 6). The ISI distribution is the distribution of first passage times from 0 to N cr open clusters with this approach.
The negative feedback terminating spikes entails a very small cluster open probability just after a global spike, from which all clusters slowly recover. Thus, slow time scales from global processes enter as a slow time dependence of the cluster IPI, puff duration and amplitude distributions.

Linear chain of states
We suggest to radically simplify the problem to reach a system which describes general properties not depending on assumptions restricting the validity of results too much and to reach possibly analytically tractable equations. We obtain a linear chain of states by averaging over all possible paths from 0 to N cr open clusters. That chain of states is indexed by the number of open clusters. The states are connected either by transition rate functions f i,i±1 (t, γ) or waiting time distributions Ψ i,i±1 (t, t − t , γ) that both result from puff properties.
The transition probabilities pick up slow time scales by their dependence on the time t since the last global spike. The probability to leave the initial state 0 and go further up at early times after a global spike is very small, such that no puffs occur early. One can almost only move to the left in the linear chain. Figure 8 shows exemplary waiting time distributions with recovery from global negative feedback for initial and later times.
Recovery from global negative feedback is described by a transient with rate γ. For the case of transition rates, we have After about t r = 5γ −1 the inhibitory effect vanishes and the system has recovered globally, i.e., The description with transition rates uses asymptotically markovian rates that are asymmetric in the sense of the recovery from global negative feedback only affecting the up rates. This is the case, because negative feedback influences the probability of clusters opening, not closing. For the case using waiting time distributions, they are the probability distributions from which a time value is drawn that determines when to jump to the next state, i.e., the time to the next opening or closing of a cluster. They depend on the time t since the last global spike, the relative time spent in a state Δt = t − t , where t is the time of entering the current state, and the current and target state i and i ± 1, respectively. The direction of the jump is drawn from the splitting probabilities, which are the relative weights, i.e., time integrals over Ψ i,i±1 at t, of possible outgoing transitions, and add up to one. This allows evaluating the system when using double exponentially distributed waiting times. The first time reaching the critical nucleus N cr is equivalent to generating a cell wide spike. We are, therefore, interested in the moments of the first passage time probability distribution to reach N cr .
Experiments show that puff times often do not strictly follow an exponential distribution, but rather a double exponential in some cases, Fig. 7. This requires using waiting time distributions instead of rate functions and to use the general master equation, which is formulated more generally in terms of probability fluxes.
Apart from choosing the state variables of the state scheme, the Ψ 's or f 's contain all the physics. This includes effects of stimulation, and positive and negative feedback from CICR on short time scales, but also recovery from global negative feedback on long time scales.
Positive feedback by CICR means the more clusters are open the larger is the open probability of the closed clusters. In mathematical terms, the λ i,i+1 are increasing functions of i. One possible choice is where stimulation strength is included via the IP 3 sensitive puff frequency λ 0 [26], N T is the total number of clusters, and k ∈ {1, 2, 3} some model parameter to quantify the strength of the positive feedback. Leftgoing rates in their most simple form account for the with a single cluster closing rate λ − . The first and second moments of the first passage time distribution from 0 to N cr open clusters can be calculated for very general f i,i±1 or Ψ i,i±1 with the method described in Falcke and Friedhoff [39]. The only requirement is that the f i,i±1 or Ψ i,i±1 can be Laplace transformed. This is possible for the Ψ i,i±1 despite their dependency on t and t − t if the t-dependency is exponential like [39]. Therefore, this method provides a basis for investigating a large variety of positive feedbacks by the choice of i-dependency of rates and waiting times, puff duration properties by the choice of left-going rates and t−t -dependency, pathway properties by the choice of [IP 3 ]-dependency, etc.

Calculating the CV
The state scheme presented in Fig. 6 and its transition rate functions (or waiting time distributions) define a (generalized) master equation, which can be solved using Laplace transforms to determine the moments of the first passage time distribution to reach state N cr [39]. The only requirement towards the waiting time distributions is that their Laplace transform exists. In case of transition rate functions, solving the master equation yields the Laplace transform (denoted by a tilde) of the probability vectorP i (s) for a process that started in state i at t = 0 as Solving the generalized master equation including the waiting time distributions gives as a solution the Laplace transforms of the probability flux vector, where A, B, E, and G are matrices that depend on the length of the chain of states N and the transition rates or waiting time distributions, f i,i±1 and Ψ i,i±1 , respectively, and r i andq contain the initial conditions, as explained in [39]. Application of this theory to a chain with state independent transitions, as it might result from a random walk, found that the CV has a minimum for a certain resonant lengthN . For a given set of parameters, CV(N ), and therefore, the value ofN can be controlled by varying the rate of recovery from negative feedback γ. The stochastic process to reach stateN for the first time is, therefore, more precise than reaching smaller or larger values of N for the first time. This is interesting in its own and for the general theory of stochastic physics, but does not resemble the robustness of the CV against changes in N cr found in Ca 2+ signalling. Here the CV is constant and independent of the various number of IP 3 R clusters per cell found in experiments, due to cell-to-cell variability. Hence, while the approach described in ref. [39] provides the tools it has not solved the problem, yet. Future modelling of Ca 2+ signalling, therefore, needs to properly define the transition rate functions f i,i±1 or the waiting time distributions Ψ i,i±1 to reproduce the measured properties of the CV, in particular its robustness against cell variability and variable conditions. CICR and spatial coupling of clusters have to be reflected by the transitions probabilities to model Ca 2+ spike generation. The probabilities for opening of more clusters, derived from Ψ i,i+1 (t, t − t , γ) or f i,i+1 (t), increases with the Ca 2+ concentration due to CICR, i.e., it increases with the number of open clusters. Due to spatial coupling by Ca 2+ diffusion, it also increases with the number of closed neighbors of open clusters and thus could pick up geometrical or spatial aspects. Ca 2+ binding molecules in the cytosol decreasing Ca 2+ diffusion would decrease the probability of opening more clusters. However, this still has to be worked out.

Similarities and differences to neural spiking
It is interesting and potentially useful to discuss in which respects Ca 2+ -spiking resembles or differs from the spiking activity of neurons, a biological problem that has been quantitatively explored by mathematical modeling to an impressive extent [48,52]. This concerns the single neuron's spontaneous activity and its characterization by interspike interval (ISI) histograms, ISI correlation coefficients, and spike train power spectra, the autonomous activity of many neurons in recurrent networks, and the encoding of time-dependent stimuli.
Obvious differences between the two forms of spiking are (i) the physical quantity that undergoes spiking (Ca 2+ concentration vs trans-membrane voltage), (ii) the time-scales and typical mean ISIs (several sec to minutes for Ca 2+ -spikes vs several to hundreds of ms for neurons), and (iii) the constancy of the spike form (the shape of Ca 2+ spikes is more variable than that of neural action potentials). A technical but important difference is the typical length of experimental recordings: neural spike trains may contain many thousands of spike pulses in a quasi-stationary setting, whereas Calcium spike trains are mostly limited to less than a hundred spikes. This poses a severe limitation for the determination of certain higher order statistics, such as interspike interval correlations. Related to this, for many sensory neurons, researchers can systematically explore the information transmission by presenting well-defined sensory (eg. acoustic, visual or electric) stimuli in the form of harmonic or broadband signals. This allows to study whether neurons preferentially encode information about slow, intermediate or fast stimulus components (see, for instance, [9,27,77]). In Ca 2+ experiments, one is mostly concerned with presenting a certain amount of signaling molecules in a step-like manner, which resembles the first experiments in neuroscience, see e.g. the famous work of Lord Adrian [1]). This, however, seems to be only a consequence of the current technical limitation and the question how the sequence of Ca 2+ spikes encode truly time-dependent signals may come into focus once more spikes can be recorded in experiment and stimuli can be better controlled.
Biophysically, it is interesting that both spiking phenomena rely on the opening and closing of ionic channels and the positive and negative feedback loops are mediated by the Ca 2+ or voltage-dependence of the opening and closing rates of these channels. The main players in the neural dynamics are the Na + and K + -selective voltage-dependent ion channels. This is described in the framework of the famous Hodgkin-Huxley model for the voltage across the neuron membrane V (t) (see standard textbooks on the topic, e.g., [24,52]) where C is the membrane capacitance and I Ext is an external current that can serve as an stimulus. The variables I K , I Na , I L describe ionic potassium, sodium and leak currents, respectively. The parameters g K , g Na , g L and E K , E Na , E L are the corresponding maximal conductances and reversal potentials. The remaining variables m(t), n(t) and h(t) are the gating variables that are of particular importance for the generation of an action potential. They are described by where x can be substituted by n, m or h. Much like in early modeling of CIRC by Ca 2+ channels the variables m and n describe two fast binding processes that activate certain channels, while h describes a slow process that inactivates the sodium-selective channels after a depolarization of the membrane. In a Ca 2+ channel model this would correspond to the fast activation due to the binding of activating Ca 2+ and IP 3 and the slow inactivation due the binding of inhibitory Ca 2+ . The positive feedback loop that is essential to understand the upstroke of the action potential is the sodium dynamics: sodium is in excess outside the cell, the opening probability of the Na + selective channels increases upon depolarization. A small depolarization will thus lead to the opening of some channels, which causes Na + ions to rush into the cell, which depolarizes the membrane further, leading to more channel openings and so forth. This positive feedback loop can be compared to the puff generation via Ca 2+ -induced Ca 2+ release (CIRC) but also to the accelerated puff generation via the global Ca 2+ concentration prior to a cell-wide Ca 2+spike.
Inactivation of Na + channels and activation of K +selective channels (with potassium being in excess inside the cell) leads to the termination of the action potential. Put in mathematical terms, negative feedback loops on a somewhat slower timescale explain the second half of the neural spike-this again is very similar on a mathematical level to the mechanism at work in Ca 2+ spiking.
There are features in the neural membrane dynamics that are sensitive to time-dependent input currents and there are features which are not. Among the latter is the exact shape of the action potential-once the voltage is sufficiently depolarized, a largely stereotypical action potential is generated. To simplify the description, one may cut out this stereotypical part of the dynamic response as it cannot contribute to the signal transmission property of a neuron and focus on what is really the signal-dependent part. This is what is done in an Integrate-and-Fire (IF) model: where the more involved dynamics of the different ion channels and corresponding currents are subsumed in a simplified function f (V ) that describes the currents up to some threshold V T . Interestingly, the particular shape of f (V ) can be obtained experimentally [3,4]. Brette [11] argues that the positive Na + feedback that sets in after a particular voltage is crossed is so abrupt that a simple linear model f (V ) = μ − V with constant parameter μ, i.e., the famous leaky Integrate-and-Fire (LIF) model, describes the sub-threshold dynamics of a real neuron best. The function s(t) could be a time-dependent signal or a stochastic processes accounting for intrinsic and/or external noise. Indeed, especially the generation of the action potential is a stochastic process due to the presence of multiple sources of noise. This includes channel noise, quasi random input from other neurons (network noise) and the unreliability of synapses [101]. Many of these noise sources can be approximated by a Gaussian stochastic process, and often the simplifying assumption of strictly uncorrelated (white) Gaussian noise is made and so will we do in the remainder of this paper. We would like to mention the limitations of this assumption. Some sorts of channel noise [42,75] and very often for network noise [30,103], fluctuations display significant correlations. Furthermore, for strong synaptic connections, the shot-noise character of neural network noise invalidates the Gaussian approximation in some cases [70].
Hence, when we want to mimic spontaneous stochastic spiking, a simple choice for the driving current is to set s(t) = √ 2Dξ(t), i.e., to use a white Gaussian noise of intensity D with ξ(t) = 0 and ξ(t)ξ(t + τ ) = δ(τ ). For concreteness, we state again the standard stochastic model, the leaky integrate-and-fire model with white Gaussian noise: Note that the spike is not explicitly modelled, instead if V (t) reaches the threshold a spike is said to be emitted at time t i = t and the voltage variable is reset to the reset value V R . The abstract spikes are described by delta-functions δ(t − t i ) and form the spike-train, i.e., the sum of all spikes: δ(t − t i ) The spike train is the essential output of an IF model and its different statistics under the influence of noisy stimulation currents has been the subject of many studies (see [13,44,51,101] for reviews of stochastic IF models). We note that the reset after a spike may occur instantaneously or with some refractory period τ ref that accounts for the temporal extent of the action potential in a conductance-based model. Several statistics of neural spike trains are also routinely studied for Ca 2+ spikes. The stationary spike rate is given by an ensemble average, r 0 = x(t) but can be also determined via a time average, is the number of spikes in the time interval T ). Statistics of the interspike interval I i = t i −t i−1 (the time between to consecutive spikes) have been already discussed for Ca 2+ spikes: the mean interval I = 1/r 0 , the coefficient of variation CV = (I − I ) 2 / I , and, of course the most complete description of the single interval, the full ISI probability density function (PDF) ρ(I). There are, however, also a number of statistics that are not as common in the study of Ca 2+ spikes but well established in the computational neuroscience community. These include (i) count statistics, especially the Fano factor F (T ) = (N (T ) − N (T ) ) 2 / N (T ) that compares the growth of the spike count's variance to its mean (see, e.g., [21,105] for studies that highlight the importance of the Fano factor and [84] for a study that investigates the Fano factor in the context of Ca 2+ spiking), (ii) the spike-train correlation function C(τ ) = x(t)x(t + τ ) − x(t) 2 that describes the probability to find a spike at time t i + τ if a reference spike occurred at time t i . This statistics bears information of the spike generation process, for instance experimentally and theoretically obtained spike-train correlation functions usually have a decreased firing probability right after a spike has occurred reflecting refractory processes similar to what is observed in Ca 2+ puffs. Often, oscillatory activity is better characterized in the Fourier domain by the spike-train power spectrum: According to the first defining equation, the power spectrum is given by the variance of the Fourier coefficientsx(f ) of the spike train in a time window T . However, according to the second equation and the Wiener-Khinchine theorem [46], it is also given by the Fourier transform of the autocorrelation function. Turning back to the interspike intervals, we mention finally the serial correlation coefficient (SCC): The high frequency limit reflects the mean firing rate r0, while the low frequency limit bears information about the variability of the spike train. In the considered case of strong mean input μ vT the interspike interval PDF can be approximated by an inverse Gaussian distribution that is fully characterized by r0 and CV . Parameters: μ = 10, D = 1, vR = 0 and vT = 1 that puts the covariance between two ISIs that are lagged by an integer k in relation to the variance of the single interval providing a number between −1 and 1. Correlations among ISIs may reflect slower processes that are at work in the driving input or in the intrinsic dynamics of the neuron. For instance, a negative SCC of adjacent intervals indicates that an ISI longer than the mean is on average followed by an interval shorter than the mean and/or the other way around. Such correlations have been found in many neurons (see [2,41] for reviews) and may lead to an improved information transmission [18,19]. Many of these statistics are related, as it can be easily demonstrated by means of the power spectrum.
We have already pointed out the relation between spike-train correlation function and spike-train power spectrum via the Wiener-Khinchine theorem. The power spectrum, however, also contains information on the interval statistics (see [22]). In the high-frequency limit of a stationary stochastic spike train, the spectrum saturates at the firing rate (the inverse of the mean ISI), lim f →∞ S(f ) = r 0 = 1/ I . If intervals are independent, i.e., if we deal with a renewal spike train, the spectrum attains also a simple form in the opposite limit of vanishing frequency: which means that by comparing the high-and lowfrequency limits we can read off how regular the renewal spike train is. More generally, the full spike-train power spectrum of a renewal point process can be obtained from the knowledge of the interspike interval probability density, more specifically, its one-sided Fourier transform,ρ(f ), via the expression [91] The spectrum can thus be calculated for the leaky IF model driven by white Gaussian noise [57] (using much earlier results for the Laplace transform of the firstpassage-time density of an Ornstein-Uhlenbeck process [23]), because in this model, the reset of the voltage erases any memory about previous ISIs and the driving noise is uncorrelated and thus does not carry memory either. Since the exact result for the power spectrum uses higher mathematical functions (the parabolic cylinder functions), it is instructive to look for a further simplification, which can be achieved if the system is in the strongly mean-driven regime μ V T − V R . In this case, the statistics of the LIF model is close to that of a perfect IF model with f (V ) = μ (omitting the leak term on the right hand side of Eq. (17)). For this model the ISI density is an inverse Gaussian probability density [47] the Fourier transform of which is a simple exponential function: In Fig. 9, we display a simulated spike-train power spectrum for an LIF model with strong mean input (μ = 10 V T − V R = 1) and highlight the limit cases, from which both the firing rate r 0 and coefficient of variation CV can be readily obtained. The simulation is also compared to Eq. (22), using as an approximate descriptionρ(f ) from Eq. (24); the approximation agrees very well for this example, because the constant drift μ is dominating the subthreshold dynamics so strongly that the LIF dynamics is close to that of a perfect IF model. Comparable power spectra have indeed been reported in Ca 2+ spiking [82].
Many extensions of the simple one-dimensional IF model have been studied analytically, such as models with time-dependent threshold [56] or models with colored [12,64,74] or non-Gaussian noise [29,65,70]. In higher-dimensional stochastic IF models we can also reproduce non-renewal behavior observed in many neurons in the form of a non-vanishing serial correlation coefficient, ρ k = 0 for k > 0. As in Ca 2+ spiking, also in neural spiking there are often slower processes at work that steer the pulse-generating process-either as a simple external control or as a feedback of the spike train onto the spike generator. This can be easily incorporated into the Integrate-and-fire framework by adding a slow variable. Consider the following example, where the membrane potential is affected by an addi- val correlations and leads to a renewal spike train). Accordingly the low frequency limit for the power spectrum of the renewal spike train has a higher power compared to the nonrenewal case. Generally, a decreased power at low frequency may improve the signal to noise ratio for potential low frequency signals, hence, improving the information transmission properties of the neuron tional negative adaptation current a(t): In the last line, we complemented the usual reset rule for the voltage by an incrementation rule for the adaptation variable a(t): it is increased by a value Δ when a spike occurs. In between spikes, according to Eq. (26) the adaptation variable will decay exponentially with the time constant τ a that is typically larger than the membrane time constant or the mean interspike interval and ranges between 50ms and several seconds. A sequence of spikes occurring in rapid succession (as we for instance observe when the neuron is subject to a depolarizing current step) leads to a large value of the adaptation variable which has an inhibiting effect on the voltage dynamics in Eq. (25)-the response to the current step will be initially a rapid increase that is followed by an adaptation to a much lower value. The IF model endowed with an adaptation current is conveniently termed adaptive IF models and can be thought of as a simplification of a conductance-based model with a Ca 2+ gated K + current. Since the adaptation current is not reset to a fixed value but increased upon spiking it carries information of past ISIs and may lead to ISI correlations. Indeed the model and some related ones as for instance the IF model with dynamical threshold have been shown to generate nonrenewal spike trains and, specifically, negative interspike interval correlations [20,58]; analytical methods to calculate these correlations have been worked out over the last years [76,78].
As a consequence of the nonrenewal character of the spike train, the power spectrum is not described anymore by Eq. (22) but by a more complex expression involving higher order interval distributions (see, e.g., [51]). While the high-frequency limit of the spectrum is still given by the firing rate (black dashed line in Abb. 10), the zero-frequency limit now involves also the serial correlation coefficients: which in the renewal case reduces to Eq. (21). The effect of the negative ISI correlations is thus to reduce power at low frequencies, which can be clearly seen by comparing the original spectrum to the power spectrum of the shuffled spike train. This effect is especially important for the transmission of slow stimuli (so far not included in the model Eqs. (25), (26)): the power spectrum of the spontaneous state is a good approximation for the noise background in the case of a time-dependent signal (e.g., a cosine signal with low frequency) being present. If noise power is reduced at low frequencies, this can increase the signal-to-noise ratio [18,19]. It is conceivable that some of the concepts reviewed here for neural spike trains may become relevant and applicable to Ca spike trains once longer recordings and better temporal control of stimuli become possible. Since in particular slower processes are at work in the intracellular Ca 2+ dynamics, models like the adaptive integrate-and-fire model that we discussed may serve as an inspiration to capture the cumulative refractoriness of Ca 2+ spikes.

Conclusion
Modelling of Ca 2+ signalling has taken place in the tradeoff between models accounting for the randomness of puffs and spikes, cell variability and measured parameter dependencies on one side and rate equation models convenient to simulate time courses on the other side in recent years. Rate equation models need further development to reproduce measured parameter dependencies. We suggest to include higher moment's dynamics derived from the Master equation to account for spike generation by fluctuations. Alternatively, approaches like spike generation as first passage of a random walk on a linear chain of states as presented in Sects. 5 and 5.1 might be used.
Stochastic theory of neuronal spiking might serve as a role model for what can be achieved with stochastic theory of Ca 2+ spiking. The main challenges ahead are to go beyond simple renewal approaches for spike generation towards ISI correlations, cumulative refractoriness and other phenomena comprising several ISI, explanation of the concentration response relation of the ISI and the robustness properties of the moment relation.
The task of mechanistic mathematical modelling in cell biology is to identify mechanisms on the basis of formulating them as hypothesis in mathematical models. Here, the agreement with experimental data serves as part of the hypothesis verification. Rate equation models fail here with respect to the agonist concentration response relation of the average interspike interval, the sensitive dependence of the average interspike interval on parameters of spatial coupling (diffusion, buffers, geometry) and of course the moment relation as a defining property of Ca 2+ spiking. Stochastic models still have to be developed to address these problems.
Such a model development may lead to answers to obvious questions in the field. Frequency encoding is one of the generally accepted and experimentally supported concepts providing meaning to Ca 2+ signals [66]. However, spike timing is random. The spectrum of a spike train with exponentially distributed ISI is flat. The absolute refractory period introduces frequencies with moderately larger power in the spectrum than the average power [82], but essentially there is no typical frequency in many IP 3 induced Ca 2+ spike sequences. Taking the large cell-to-cell variability at a given agonist concentration into account, there is no defined relation between agonist concentration and Ca 2+ spike frequency applying to all cells of a given type, but each cell has its own relation. How can frequency encoding work with these properties of spiking? What are the reasons for the large cell variability and what does it mean? Addressing these questions requires models that faithfully reproduce the properties of spike sequences including their fluctuations but also have predictive power, e.g., by telling us how the spike statistics will vary if biophysical parameters are changed.
Funding Open Access funding enabled and organized by Projekt DEAL.

Author contribution statement
All authors wrote the introduction (Sect. 1) and conclusion (Sect. 7). Sections 2-5 were written by VNF and MF, and Sect. 6, including simulations and figures, by LR and BL. All authors edited the manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.