On the phase space structure of IP3 induced Ca2+ signalling and concepts for predictive modeling

The correspondence between mathematical structures and experimental systems is the basis of the generalizability of results found with specific systems and is the basis of the predictive power of theoretical physics. While physicists have confidence in this correspondence, it is less recognized in cellular biophysics. On the one hand, the complex organization of cellular dynamics involving a plethora of interacting molecules and the basic observation of cell variability seem to question its possibility. The practical difficulties of deriving the equations describing cellular behaviour from first principles support these doubts. On the other hand, ignoring such a correspondence would severely limit the possibility of predictive quantitative theory in biophysics. Additionally, the existence of functional modules (like pathways) across cell types suggests also the existence of mathematical structures with comparable universality. Only a few cellular systems have been sufficiently investigated in a variety of cell types to follow up these basic questions. IP3 induced Ca2+signalling is one of them, and the mathematical structure corresponding to it is subject of ongoing discussion. We review the system's general properties observed in a variety of cell types. They are captured by a reaction diffusion system. We discuss the phase space structure of its local dynamics. The spiking regime corresponds to noisy excitability. Models focussing on different aspects can be derived starting from this phase space structure. We discuss how the initial assumptions on the set of stochastic variables and phase space structure shape the predictions of parameter dependencies of the mathematical models resulting from the derivation.

(Received 31 December 2017; accepted 12 March 2018; published online 13 April 2018)   The correspondence between mathematical structures and experimental systems is the basis of the generalizability of results found with specific systems and is the basis of the predictive power of theoretical physics.While physicists have confidence in this correspondence, it is less recognized in cellular biophysics.On the one hand, the complex organization of cellular dynamics involving a plethora of interacting molecules and the basic observation of cell variability seem to question its possibility.The practical difficulties of deriving the equations describing cellular behaviour from first principles support these doubts.On the other hand, ignoring such a correspondence would severely limit the possibility of predictive quantitative theory in biophysics.Additionally, the existence of functional modules (like pathways) across cell types suggests also the existence of mathematical structures with comparable universality.Only a few cellular systems have been sufficiently investigated in a variety of cell types to follow up these basic questions.IP 3 induced Ca 2þ signalling is one of them, and the mathematical structure corresponding to it is subject of ongoing discussion.We review the system's general properties observed in a variety of cell types.They are captured by a reaction diffusion system.We discuss the phase space structure of its local dynamics.The spiking regime corresponds to noisy excitability.Models focussing on different aspects can be derived starting from this phase space structure.We discuss how the initial assumptions on the set of stochastic variables and phase space structure shape the predictions of parameter dependencies of the mathematical models resulting from the derivation.V C 2018 Author(s).All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/).https://doi.org/10.1063/1.5021073IP 3 induced Ca 2þ signalling is one of the most versatile and universal cellular signalling systems and a popular model system in non-linear dynamics for pattern formation in noisy systems.We discuss the experimental evidence allowing for identification of the mathematical structure to which it corresponds, and a variety of concepts for deriving simplified models from it.

I. INTRODUCTION
In spring 1995, I (MF) joined John (Jack) L. Hudson's lab in Charlottesville, Virginia, to work with him on dynamic clustering of globally coupled non-linear oscillators or a topic from pattern formation far from thermodynamic equilibrium.James D. Lechleiter and Patricia Camacho were in Charlottesville at this time, too.James had just published his results on the effect of energizing mitochondria on Ca 2þ waves in Xenopus oocytes, 1 which had several aspects very interesting for the theory of pattern formation.According to that theory, free ends of waves in excitable systems should either form a spiral or recede.The free ends of Ca 2þ waves with energized mitochondria neither formed spirals nor receded but showed different dynamics.Jack suggested working on these patterns.This was my first biophysical project and it redirected my career.Jack worked experimentally and developed also the mathematical models explaining his experiments.His high standards and expectations towards theory close to experiments substantially influenced all of my later scientific and educational work.
The first years of this biophysical research led to results on spiral instabilities, spiral pattern regimes, and generation and annihilation dynamics, 2 but could not explain Lechleiter's experiments.The underlying mathematical structure of the model did not correspond to the experimental system.When we replaced the model with a direct transition from excitability to an oscillatory regime by a model with a direct transition from excitability to bistability, 3 it explained not only the mitochondria experiments 4 but also experiments which were not taken into account when it was developed. 5his exemplifies how the basic mathematical structure of a non-linear dynamical system defined by its set of bifurcations and their relation, often called the bifurcation diagram and phase space structure, is essential for the predictive power of a theoretical description.
In physics, the fundamental equations like Newton's first law, the variational principles of classical mechanics, or the Schr€ odinger equation of quantum mechanics have been developed with simple examples.Nonetheless, they are stunningly predictive far beyond the systems used in their formulation.This predictive power originates from a correspondence between the experimental objects and mathematical structures.The mechanics of macroscopic objects corresponds to variational principles and differential equations, and the behaviour of microscopic objects corresponds to operator theory in Hilbert spaces.The identification of the correct mathematical structure corresponding to an observation provides predictive power to a mathematical theory in science.Mathematical models formulated within mathematical structures not corresponding to the observations still may reproduce the measurements used for their development but rarely are predictive beyond them as illustrated by the history of atom models.
In general, the biophysics of cells has to obey the basic laws of physics-the first principles.But cells consist of many components and interactions and therefore specifying the fundamental equations of physics to a living cell is close to impracticable.The approach of theoretical biophysics is consequently to consider the components and interactions assumed to be most relevant for a specific process of interest and to verify the assumptions retrospectively by contrasting model predictions with experimental results.But does the lack of models derived from first principles for cellular behavior also mean that the correspondence of mathematical structures to observations has no meaning in cellular biophysics?The predictive power growing out of it makes it worth to follow up on this only seemingly philosophical question.
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.][8] 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þ binding, such that Ca 2þ released from the ER through one channel increases the open probability of neighboring channels.This positive feedback of Ca 2þ on 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, which is in the range of 3-4 orders of magnitude.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.0][11][12][13] 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 lm 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. 146][17] Openings of single IP 3 Rs (blips) may trigger collective openings of IP 3 R within a cluster (puffs), while Ca 2þ diffusing from a puff site can then activate neighboring clusters, eventually leading to a global, i.e., cell wide, Ca 2þ spike. 13,16,18,19Repetitive sequences of these Ca 2þ spikes encode information that is used to regulate many processes in various cell types. 6,20,216][27] The molecular mechanism of this feedback is pathway and cell type specific and not always known although a negative feedback on the IP 3 concentration might be involved. 28,29Hence, the negative feedback that determines the time scale of interspike intervals is different from the feedback contributing to interpuff intervals and requires global (whole cell) release events.
At very strong stimulation, cells exhibit a raised Ca 2þ concentration of much longer duration than spikes which may oscillate, 30,31 burst, 32,33 or is rather constant. 1,34,35ypically, the amplitude of these oscillations is smaller than the spike amplitude.In the following, we review our current understanding of experimental results on Ca 2þ signaling and how it illustrates the relation between mathematical structures and observations in biophysics.

II. EXPERIMENTAL RESULTS ON THE PHASE SPACE STRUCTURE AND DYNAMICAL PROPERTIES OF IP 3 INDUCED Ca 21 RELEASE
The pathway exhibits local Ca 2þ release through individual channel clusters at low [IP 3 ], spiking at intermediate [IP 3 ] and an elevated cytosolic [Ca 2þ ] i at high [IP 3 ].A basic observation in all experiments is that cell-to-cell variability with respect to Ca 2þ spiking behavior is large but not completely arbitrary.It obeys some preserved characteristics, which have been confirmed for all cell types in which they have been investigated.We will focus on these general characteristics since they obviously reflect essential system properties.
It is convenient for the presentation of experimental results to introduce also a few mathematical concepts.In mathematical terms, intracellular Ca 2þ dynamics are described by reaction-diffusion equations like @X @t ¼ DDX þ FðX;r; t; pÞ; where X is a vector of concentrations, t is time, r is the space coordinate, D is a diagonal diffusion matrix, D the Laplace operator, FðÁÞ is a non-linear function describing the local dynamics, and p is a vector of parameters.X comprises free cytosolic Ca 2þ , Ca 2þ bound to Ca 2þ -binding molecules, IP 3 , and free and bound Ca 2þ in the lumen of the ER and mitochondria in a rather general formulation of the dynamics.
In general, non-linear dynamics reaches asymptotically attractors in phase space which may be stationary states or manifolds of higher dimension.Attractors with higher dimension like limit cycles, tori, or even chaotic attractors potentially describe the Ca 2þ spiking behaviour.They may be caused by the dynamics of spatial modes [eigenfunctions of the linearized rhs of Eq. ( 1)] or by the local dynamics, 37 i.e., may occur with DX 0 also.Spatial modes have been observed with the Ca 2þ dynamics of excitation contraction coupling in cardiac myocytes, 38,39 which is a driven system in terms of dynamical systems theory.However, there is no experimental evidence for attractors of the autonomous and/ or IP 3 induced intracellular Ca 2þ dynamics caused by spatial modes, and hence we can focus on properties of the local dynamics.
The local dynamics of Eq. ( 1) are the behaviour of the IP 3 R clusters.The majority of the modelling literature assumes oscillatory local dynamics in the spiking regime since measured spikes are repetitive.Indeed, spike sequences even with a CV of 0.3-0.4 of the ISI appear surprisingly regular in visual inspection.However, a closer look could not confirm this assumption. 24,40lusters are dynamically coupled by Ca 2þ diffusion, which needs to be reduced for investigations focussing on the local dynamics.Such an uncoupling can be achieved by high intracellular concentrations of the Ca 2þ buffer EGTA.The elemental event of the local dynamics is the stochastic opening of channels in a cluster.The first open channel entails with some probability opening of more channels in the cluster causing a puff.Puffs last typically a few tens of ms but with large scatter. 13,41The probability of triggering calcium puffs is linearly related to the number of IP 3 R in a cluster. 42Puff sequences at a given cluster exhibit some correlation between amplitude and subsequent interpuff intervals, a weak correlation between interpuff intervals and subsequent amplitude, but no detectable correlation between consecutive amplitudes. 41Both puff amplitude and frequency increase and saturate with increasing stimulation of cells. 42ypical interpuff intervals last a few seconds, 13,24,41,42 and 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. 24A 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. 24The ISI time scale has only been observed on cell level.Consequently, spikes are a collective phenomenon requiring coupling of clusters.Another set of experiments demonstrated that indeed the average ISI depends sensitively on the intracellular buffer concentration modulating the strength of spatial coupling. 40This confirms the results of the analysis of the local dynamics.
These experimental results are supported by theoretical investigations.The Ca 2þ concentration at closed clusters is the resting concentration in the range of 100 nM.Detailed simulations of the concentration dynamics in the immediate vicinity of channels 14 showed that concentrations at open channels are high (>20 lM).4][45] 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). 17f channels are sufficiently sensitized for Ca 2þ binding, puffs may cooperate to set off a global release spike spreading from the initiating site into the cell in a wave like manner.Waves occur if a critical number of releasing clusters is reached. 16,46,47The randomness of puffs causes randomness of spike timing with a linear relation between the standard deviation r of interspike intervals and the average T av as shown in Fig. 2 and further for 8 cell types and 10 conditions 27,40,[48][49][50] (see also Ref. 51).The slope a of this relation between SD and average is the same for all cells of the same type stimulated with the same agonist 27,40,48,49,52 and robust against changes in stimulation strength, 27 pharmacological perturbations, 27 changes in buffering conditions, 40 and the large cell variability.It has been verified even in cells not exhibiting clustering of channels and puffs. 49Values of a are, for example, about 0.2 for hepatocytes stimulated with vasopressin, 0.25 for human embryonic kidney (HEK) cells stimulated with carbamyl choline (CCh), 0.37 for hepatocytes stimulated with phenylephrine, 27 0.7 for processed lipoaspirate (PLA) cells, 52 and close to 1 for spontaneously spiking astrocytes. 40Consequently, the standard deviation is of the same order of magnitude as the average ISI.
The standard deviation of ISI of oscillatory systems moving on a limit cycle in phase space and perturbed by noise is typically smaller than the values measured for Ca 2þ spiking, 53 and/or the cumulant relation may exhibit a negative slope. 53Varying parameter values across the range covered by cell variability and the perturbations applied in two studies 27,40 cause loss of a unique relation between r and T av 53 with these oscillatory systems since the period and the noise causing the standard deviation are determined by different processes.Thus the robustness of a against cell variability and perturbations can hardly be reconciled with an oscillatory dynamics since all these parameter variations against which a is robust would need to affect the processes setting the average and the processes setting the SD in exactly the way conserving the CV.But since spike generation is stochastic, the parameters control only the spike generation probability, and the type of stochastic process-e.g., inhomogeneous Poisson-fixes the relation between T av and r. 54 The second parameter of Eq. ( 2), the absolute refractory period T min , was also found to be the same for all individual cells of the same type stimulated with the same agonist. 27,40hen T min has passed, the puff probability recovers from 0 gradually to its asymptotic value.This slow recovery delays initiation of the next spike.That spike may occur during recovery, if the asymptotic spike generation probability is large compared to the recovery rate, or after recovery in the opposite case.The contribution of this stochastic part of the ISI to the total average ISI has been thoroughly investigated and is well known.9][50] The slower the recovery, the smaller is the ratio of SD to average ISI (coefficient of variation CV). 54he wave-nucleation like generation of global release spikes as well as the ISI statistics strongly suggests excitability as the dynamic regime of IP 3 induced Ca 2þ spiking in agreement with the analysis of the local dynamics.Excitable systems exhibit a stationary state which is stable against small perturbations.Perturbations above the excitation threshold are amplified to a transition to the excited state.The stochastic behavior of channel clusters causes incidental local transitions to the excited state, which then spreads with some probability into the whole cell.The resulting large fraction of open clusters-i.e., a release spike-causing a high Ca 2þ concentration and high open probability are the excited state of Ca 2þ dynamics.This state is terminated by negative feedback acting on a slower time scale than the excitation.The probability for generating this supercritical local excitation fixes the average stochastic part T av À T min and the standard deviation r.
The complete distribution of ISI cannot be easily determined from experimental data since measured spike trains are not longer than about 60 ISI.Fusion of ISI sequences normalised by T av have been used as a surrogate data set and led to skewed distributions with an absolute refractory period. 55More sophisticated methods based on the time rescaling theorem and Kolmogorov-Smirnov tests for comparison of measured and hypothetical distributions identified an inhomogeneous Gamma distribution as the most likely experimental ISI distribution with time dependent stimulus. 56Distributions of ISI obtained with these methods and constant stimulation are shown in Fig. 3.
The response of the average ISI to stimulation with extracellular agonists has features applying to all of the four plasma membrane receptors for which it has been investigated. 27On that basis, we assume them also to be general features of the system.T min is not affected by stimulation, as we have learned from the robustness properties of Eq. ( 2), already.Stimulation controls the average stochastic part T av À T min of the ISI.The concentration response has been  established by applying steps in the concentration a of the stimulating agonist. 27The change of the average stochastic part of the ISI due to this concentration step is proportional to the average stochastic part at the lower agonist concentration T av1 (Ref.27) Analysis of measurements revealed that b does not depend on the agonist concentration, 27 which entails an exponential dependency on a T ref st is the average stochastic part measured at the reference concentration a ref .This prefactor of the exponential is cell specific and picks up all the cell variability.The constant c in the exponent is the same for all cells of a given cell type stimulated with the same agonist.Equation ( 4) does not bear directly information on the dynamic regime of IP 3 induced Ca 2þ spiking, but it defines clear constraints to its theory.

III. BASIC REQUIREMENTS AND CONCEPTS FOR MODELLING OF IP 3 INDUCED INTRACELLULAR Ca 21 DYNAMICS
A comprehensive monograph reviewing modelling of intracellular Ca 2þ dynamics has recently been published. 57ere, we would like to fill a void in the literature by a critical reflection on the framework of model derivation and the approximations coming with modelling concepts used in the biophysical literature.
The essence of the 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.• The dynamics of individual clusters are not oscillatory on the time scale of ISI.
• Cell-to-cell variability of average ISI is large.
• The spiking regime obeys Eqs. ( 2), (3), and (4) with T min , a, and c being cell type and pathway specific but not subjected to cell variability. 58 ISIs depend sensitively on parameters of spatial coupling.
The high stimulation regime is not in this list since the behavior is cell type dependent-it might be stationary or oscillatory.
The hierarchical organization of CICR carries the randomness of individual channels onto the level of cell-wide spikes via the stochastic puff dynamics of clusters (Fig. 1).The random channel state changes are the source of noise.Consequently, the master equation for the probability of the microscopic states of the system is the starting point for an exact theory.In the most general case, that would comprise the position and number of Ca 2þ ions and Ca 2þ binding molecules and the state of the channels and pumps.While the master equation as the starting point for a theory defines concepts and methods to be used, solving it is not practical in the end.Hence, probabilistic theories usually start from formulations of the state dynamics eligible for simulating trajectories in phase space.

A. Simulations
The diffusion coefficients of Ca 2þ and Ca 2þ binding molecules are sufficiently large to establish the deterministic concentration profile on the time scale of typical channel state changes due to the frequent sampling of space by thermal motion.The number of SERCA molecules is orders of magnitude larger than the number of Ca 2þ channels.Hence, we can describe diffusion, the reactions involving cytosolic Ca 2þ binding molecules, and the SERCA flux by reactiondiffusion equations like Eq. ( 1).The random opening and closing of channels causes time dependent source terms in the partial differential equation for the Ca 2þ concentration.We illustrate that with a simple model comprising cytosolic Ca 2þ c, one Ca 2þ buffer b (Ca 2þ bound form), and the ER Ca 2þ concentration e Here, b t denotes the total buffer concentration, k þ and k -the binding and dissociation rate, V p is the maximum SERCA pump flux, and the ratio of cytosol to ER volume.We have approximated the shape of a channel mouth by a spatial d-function and the time course of a single opening by a temporal d-function.N is the number of channels, ri is the location of the ith channel, and t i;j f g the sequence of its openings before time t.The sequence of time points of openings is determined by Markov chain Monte Carlo simulations for the state of each individual channel.The simulations are based on state schemes; an example is shown in Fig. 4.[64] B. Distributions and their moments Probability distributions for stochastic variables are the natural way to characterize stochastic systems.They are the solutions of the master equation.However, we need to simplify the system to obtain equations we can solve.These simplified systems can be informed by the general properties listed above.We know about the ISI distribution that it should exhibit an absolute refractory period and a linear relation between standard deviation and average.
The formulation of the problem in terms of Eq. ( 5) and Markov chains can also serve as starting point for analytical calculations or derivation of simplified models.The robustness of spike generation with respect to cell variability and perturbations demonstrates that it cannot depend on very specific parameter values or other details.Hence, simplifications should not destroy the basic characteristics of the system.At the same time, the large cell variability entails requirements on the theory.With each experiment comprising a population of cells, we sample a phase space volume large enough for accommodating this cell variability.Hence, the qualitative properties of IP 3 induced Ca 2þ dynamics listed above must not depend sensitively on the value of parameters distinguishing individual cells.These parameters comprise protein concentrations, 65 the number of clusters, their spatial arrangement, diffusion properties, and more. 13 suggestion for calculating the ISI distribution has been made in this spirit. 54It starts from the wave nucleation character of spike generation.All clusters are closed at the end of a spike.Each opening cluster entails a sphere of increased Ca 2þ concentration around it.We indicate that by the orange spheres in the red round cells above scheme (8)  8) can be directly calculated from interpuff interval and puff duration distributions. 54Such an approach is able to explain the cumulant relation Eq. (2). 54 lot remains to be done even with such a simple approach.The dependency of the transition probability on the numbers of open clusters and the parameters of spatial coupling has not been worked out analytically, yet.Also, the effect of the recovery from the negative feedback terminating spikes has not yet been described analytically in this approach but with phenomenological Ans€ atze or stochastic simulations only. 40,53,55,66A new approach to this problem has been suggested recently but has not been specified to Ca 2þ spiking, yet. 67Derivation of the concentration response relation Eq. ( 4) with this approach has neither been attempted, yet.

C. Rate equations
][72][73][74] The derivation of rate equations implies averaging over the state distribution dynamics defined by the master equation.The spatial character of spike initiation renders the averaging difficult.Another (related) conceptual problem arises from the fact that the dynamics on cell level is still noisy.In contrast, the more frequent situation in the derivation of cellular dynamics encounters noise on the molecular level only.The population average carried out in the master equation of such systems during the derivation of rate equations is an average over the molecules in a single cell.The large number limit guaranteeing the validity of deterministic rate equations applies to the cell level.
With IP 3 induced Ca 2þ spiking, this limit does not apply to the cell level since cell behavior is noisy.The average needs to be carried out across an ensemble of identical cells.Consequences of these considerations can be illustrated by a comparison to existing rate equation models.We prepare this comparison by reconsidering the rate equation derivation of the most simple stochastic process-radioactive decay of atoms.The stochastic variable is the number N a of atoms.We denote the probability per unit time for decay of a single atom with k.The atom number N a obeys for large initial numbers N i the exponential function N a ¼ N i e Àkt .Each decaying atom is in a stationary state till it decays, and there is no process setting the time point of its decay.However, if we ask for the time t r required till a specific number of atoms N r remains in the deterministic limit, it is set by N i , N r and k (t r ¼ k À1 lnðN i =N r Þ).The process setting the time scale is the continuous decrease in N a down to N r .
Rate equation models derived by averaging on the molecular level and assuming deterministic behavior on cell level usually require specific processes to set the time scale of ISIs.That might be a rising fraction of channels recovered from inhibition, an approach to a critical Ca 2þ concentration or the rise of receptor sensitization. 57,74However, the noisy behavior of Ca 2þ spiking entails different determinants of the average ISI. Figure 5 illustrates some differences between the rate equations obtained by assuming deterministic cell behavior and noisy behavior on cell level.The time courses were obtained from simulations of a purely deterministic model 75 (black) and a noisy excitable version of it 53 (red).Both systems respond with a spike to the perturbation.The noisy system generates subsequent spikes some time after the stimulated one [red line in Fig. 5(a)].The deterministic rate equation model stays in a stationary state after the initial perturbation without generating a second spike [black line in Fig. 5(a)].During the time T div , the dynamics is completely noise dominated.This illustrates that completely analogous to radioactive decay, there is no deterministic process on the level of the individual cell setting its ISI after recovery from the previous spike.
Spiking is lost in the rate equations since T div diverges due to averaging on the molecular level.Thus also the dependency of the ISI on the parameters characterizing the noise and spatial coupling is lost.Most rate equation models tune parameters to an oscillatory regime to establish spiking [Fig.5(b)].The interspike interval is then dominated by the time required to reach the threshold for CICR.This entails parameter dependencies of the ISI different from the ones of noise driven dynamics.
The sketch of T div for the excitable model in Fig. 5(a) applies when the asymptotic spike generation probability reached after recovery is smaller than the recovery rate from negative feedback.The medium and long ISI data in Xenopus oocytes 18 and spontaneously spiking astrocytes and microglia cells 40,66 are experimental realizations.Their recovery phase from negative feedback is substantially shorter (a % 1) than the average ISI. 40,66The effect of noise on time scales and parameter dependencies is also substantial if the recovery phase and the average ISI are of comparable length. 27,40,66n summary, averaging on the single cell level across molecules and clusters eliminates the noise generating the spike.The rate equations for this average do not reflect the spike generating mechanism since usually an oscillatory regime is then used to "rescue" spiking.However, averaging over a stochastic ensemble of cells defined by a cellular spike generation probability distribution allows for including the average of the noise generated time scale and its parameter dependencies, and can thus reflect the spike generation mechanism.
Deriving rate equations in a way reflecting the spike generation mechanism is an open problem and has not been attempted, yet.Suitable concepts might be inspired by the integrate-and-fire models of neuronal dynamics starting from an expression for the spike generation probability on cell level.Investigations on globally coupled noisy excitable systems might be specified to Ca 2þ dynamics. 76Another very promising approach includes higher moments in the derivation. 77arameter dependencies and the mathematical structure of models can also be restricted by Eqs. ( 3) and (4).Stochastic simulations of the excitable regime of the frequently used DeYoung-Keizer-model reproduced Eq. ( 3) but not Eq.(4).FIG. 5. Time scales set by noise are not captured by current deterministic rate equations.(a) Caricature of a Ca 2þ time course as produced by deterministic rate equations (black) and by a corresponding noisy system (red) after an initial perturbation (arrow) based on model simulations. 53,75The noisy system generates subsequent spikes some time after the previous one.During the time T div , the deterministic rate equations are in a stationary state without generating a second spike.(b) The interspike interval is dominated by the time required to reach the threshold of CICR (blue) in the oscillatory regime of deterministic rate equation models.The dependency of the ISI on the parameters characterizing the noise is lost.

IV. CONCLUSION
While detailed multiscale simulations can mimic experimental observations in a rather flexible manner, 19,53,[59][60][61][62][63][64] neither the current state of the stochastic theory nor the rate equation models live up to the request for predicting experimental outcome beyond the examples used for model derivation.This indicates that we have not yet understood how to derive the appropriate models.Based on the accordance of experimental and multiscale simulation results, we come to the conclusion that a reaction diffusion system with a local dynamics in a noisy excitable regime must be the starting point of the derivation of predictive models since it is the mathematical structure corresponding to the observations.9][80] On the basis of early interpretations of experimental results, it became one of the prototypical cellular limit cycle oscillators.The recent experimental results reviewed in this study revealed that the repetition of spikes is caused by noise instead of a limit cycle or torus in phase space.Derivation of predictive and simple models starting from this noisy spatially extended excitable system is a task reaching beyond the specific biological system.Hence, this classic still poses theoretical problems interesting and challenging for the whole field of nonlinear dynamics.

2 .
Variability in Ca 2þ signals.(a) The transient cytosolic Ca 2þ concentration of an astrocyte stimulated with 10 lM adenosintriphosphat (ATP) (upper panel) exhibits some variability as indicated by the variable individual ISIs (lower panel).(b) An astrocyte of the same experiments shows slower and more irregular spiking illustrating cell-to-cell variability.(c) The systematic analysis of the standard deviation r of ISI versus the average ISI T av for HeLa cells stimulated with 100 lM histamine reveals a linear dependence in accordance with the moment relation (2) where each data point corresponds to the characteristic of an individual cell.(d) The r-T av relation of astrocytes stimulated with 10 lM ATP exhibits also a linear dependence with a different slope than HeLa cells.T av À T min is the average stochastic part of the ISI.
The local rise in Ca 2þ increases the open probability of the open cluster's neighbours.A spike occurs when a critical number N cr of open clusters is reached via one of many possible paths of cluster openings.Hence, the ISI calculation can be formulated as a first passage problem from 0 to N cr open clusters.The first passage time distribution corresponds to the ISI distribution for stationary spike trains.This approach radically simplifies the system into a state space defined by the number of open clusters only. 54That implies averaging over all N path possible paths from 0 to N cr open clusters. 54(8) The transition probabilities from k to k þ 1 open clusters are determined by the probability that k open clusters open another one, and from k to k À 1 that a cluster closes.The transition probabilities W i;k in state scheme (

FIG. 4 .
FIG. 4. This state scheme of the IP 3 R originally published by Siekmann et al. 68 is comprised of two modes.One is the drive mode containing three closed states C 1 , C 2 , C 3 and one open state O 6 .The other is the park mode which includes one closed state C 4 and one open state O 5 .The rates of statetransitions within each mode are constants.a and b are the rates connecting the two modes and depend on Ca 2þ in a highly dynamic manner.Reprinted with permission from P. Cao, M. Falcke, and J. Sneyd, Biophys.J. 112, 2138-2146 (2017).Copyright 2018 Elsevier.69 36 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).Reprinted with permission from A. Skupin, H. Kettenmann, and M. Falcke, PLoS Comput.Biol.6,e1000870 (2010).Copyright 2010 ICSB.36 FIG. 3. ISI distributions P(ISI) for two spike trains measured with HEK cells stimulated with 100 lM CCh.The differences between the distributions illustrate cell variability.The experimental data are from the experiments published in Ref. 27, and the fitting method is explained in Ref. 56.