Hierarchic Stochastic Modelling Applied to Intracellular Ca2+ Signals

Important biological processes like cell signalling and gene expression have noisy components and are very complex at the same time. Mathematical analysis of such systems has often been limited to the study of isolated subsystems, or approximations are used that are difficult to justify. Here we extend a recently published method (Thurley and Falcke, PNAS 2011) which is formulated in observable system configurations instead of molecular transitions. This reduces the number of system states by several orders of magnitude and avoids fitting of kinetic parameters. The method is applied to signalling. is a ubiquitous second messenger transmitting information by stochastic sequences of concentration spikes, which arise by coupling of subcellular release events (puffs). We derive analytical expressions for a mechanistic model, based on recent data from live cell imaging, and calculate spike statistics in dependence on cellular parameters like stimulus strength or number of channels. The new approach substantiates a generic model, which is a very convenient way to simulate spike sequences with correct spiking statistics.


Semi-Markov processes
We define an appropriate stochastic process X n (t) we can use as a model for the intracellular calcium dynamics. We start with a finite sample space Ω = {0, 1, ..., S} consisting of all possibly macroscopic states the process can visit. By T n ∈ R + we denote the subsequent transition epochs between these states with T 0 < T 1 < ... < T n . We can now formulate transition probabilities according to The process is temporally homogeneous by noting the independence of Q i,j (t) from n, and we also assume Q i,j (0) = 0.
We will briefly summarise the reasoning behind that terminology. By defining one gets the transition probabilities for the embedded Markov chain with normalisation condition j p i,j = 1. Next we define the following probability distributions: The successive visits of the process X n form a Markov chain with transition probabilities p i,j , whereas the length of the sojourn time intervals [T n , T n+1 ) are given by the distribution functions G i,j (t). If these distributions can be written as G i,j (t) = 1 − exp(− j q i,j t) ≡ G i (t), than the process is a pure Markov process with rates q i,j , the sojourn times are exponentially distributed and independent of the next state. In that case the process is memoryless, which means that for s, t > 0 only holds for exponential distributions.
In conclusion, a semi-Markov process (X n , T n ) still fulfils the Markov property with respect to the subsequent state transitions, but allows arbitrary (with respect to j lim t→∞ Q i,j (t) = 1) sojourn time distributions and hence looses its memorylessness with respect to the transition times T n+1 − T n . It is therefore the ideal frame work to apply the desired non-exponential waiting times often found experimentally, now exactly defined by: Therefore, as introduced in the main text, the term conditioned waiting time seems appropriate. By noting that the transition probabilities of the embedded Markov chain are given by p i,j = ∞ 0 Ψ i,j (t)dt, the semi-Markov process is completely defined by a set of conditioned waiting time densities. Also note that by definition of the transition probabilities in Eq. 1, the time variable t does not correspond to a system time, but describes the wait after the last transition T n .

Detailed derivation of the first passage time density
We want to outline the computational strategy to compute mean first passage times (FPT) for arbitrary discrete state semi-Markovian systems. The particular simple case of a linear chain is used for the stochastic Ca 2+ model, and we will show its solution as a special case at the end of this section.
In the main text we introduced the Laplace transformed non-Markovian master equation: containing the Laplace transformed probability fluxesĨ j il . The FPT is given by (see also main text)F with the moments: This equation yields for n=1 the mean FPT.
Next we calculate the Laplace transformed fluxes. We start with the formula for a single flux [4]: which is a convolution of the conditioned waiting time to go to state l and all incoming fluxes to the state i. The f j il are the initial functions, dealing with the initial state j of the whole system (see also main text). By writing all possible fluxes as a vector I j (t) and the conditioned waiting time densities in an appropriate matrix Ψ(t), the system of integral equations can be written as with f j (t) as vector of the initial functions.
The transition network of the system is completely determined by the fluxes (see also [5]). A standard technique for solving such integral equations with a convolution kernel is by using the Laplace Transform Equation 10 reads in Laplace space: Basic algebra yieldsf This is an inhomogeneous linear system of equations for the Laplace transformed fluxesĨ il which can be solved exactly by standard algebraic methods.
The solution yields all Laplace transformed fluxes as functions of the Laplace transformed waiting times. As can be seen in Eq. 7 for the FPTF i,j (s), there are two sets of fluxes needed, reflecting the initial states i or j respectively.
Inserting these fluxes into Eq. 7 gives the desired Laplace transformed FPT densityF i,j (s).
In case of a linear chain with N = K + 1 states, we write for the flux vector I j = (I j 01 , I j 10 , I j 12 , ..., I j K−1,K , I j K,K−1 ) and the matrix elements Ψ il are determined by the set of equations given in Eq. 10 (see also section below for an example). The equation 7 for the FPT simplifies to: For our minimal stochastic Ca 2+ model we have N = 5 and the FPT is given Note that we did not have to specify the explicit form of the conditioned waiting times for the entire calculation up to here. In order to calculate the mean FPT or higher moments of the FPT density Eq. 8 can be applied.

Minimal Markovian Example
The analytic approach presented in the last section for solving the first passage time problem shall be illustrated in a simple example. The system we want to consider is a three state continuous time Markov chain, with transitions We recall from the theory of continuous time Markov chains [6] that the sojourn time distributions are the following: By considering the transition probabilities of the embedded (discrete time) Markov chain, namely p 01 = 1, p 10 = γ α+γ , p 12 = α α+γ and p 21 = 1, we can define our conditioned waiting times according to the last section about semi-Markov processes: After Laplace transform we obtain: We now solve for the vector of the fluxesĨ j = (Ĩ j 01 ,Ĩ j 10 ,Ĩ j 12 ,Ĩ j 21 ) the matrix equation 13. With the initial function vectorsf 0 = (Ψ 01 , 0, 0, 0) andf 0 = (0, 0, 0,Ψ 21 ) and the matrix of conditioned waiting times Next we plug these fluxes into the equation for the Laplace transformed FPT density Eq. 14 calculated in the last section: .
The mean first passage time is given by Because the structure of the solution in terms of the Laplace transformed FPT density is quite simple for that minimal system, it is possible to use the inverse Laplace transform to obtain the FPT density directly: where The result is a bi-exponential density. If we impose a time-scale separation, meaning α λ, the second exponential goes to zero for growing λ. This leads to a convergence to an exponential distribution of the FPT, as observed in the parameter studies for the stochastic calcium model.

Laplace transformation of the generalised exponential distribution
The generalised exponential (GE) distribution is a generalised form of the human mortality distribution discovered by Gompertz and Verhulst in the first half of the 19th century [7,8]. For the analytic solution of the first passage time problem formulated in the main text, the Laplace transformation of the GE density function and its survival function are needed.
The GE density itself reads We want to solve its Laplace transform whis is formally given by the following integral: We start by substituting y = e −t , which gives Next we use the substitution x = y λ : Now we use the definition of the beta function, written as an integral, it reads: There is also a form involving the Euler Γ-function: By using l = s λ + 1 and substituting k = α, we obtain from Eq. 25: .
The Laplace transformation for the survival function can be done in analogy. This time we want to solve where we already used the trivial Laplace transformation of a constant. We substitute y = e −t followed by the second substitution x = y λ again to end up Now we recall the beta function (Eq. 27) and finally substitute l = s λ and k = α + 1 to obtain the solution: Note that the powers ofΨ o needed for the construction of the conditioned waiting times can be written as These are just sums of the original terms with a rescaled shape parameter kα, and therefore their Laplace transforms are sums of the accordingly rescaled results given above in Eq. 32.
5 Details for the minimal mechanistic Ca 2+ model 5.1 Closing probability ψ c (t).
The probability density of the waiting time until an IP 3 R closes again, ψ c (t), can be derived from the experimental finding that the individual channels in a cluster close independently with closing rate γ, which does not depend on the Ca 2+ concentration. It is a molecular property of IP 3 R [5,9]. Based on this, we can write the waiting time density for closing of an IP 3 R cluster with, on average, N ch channels involved in a puff: The channel closing rate γ has recently been determined by total internal reflection fluorescence (TIRF) microscopy and is 17 s −1 in SH-SY5Y cells [9].
Equation 34 fits data from TIRF microscopy [9] and leads to Ca 2+ spike statistics close to experimental data if incorporated into the hierarchic stochastic model [5].
The form of ψ c (t) exactly constitutes a GE distribution with parameters N ch and γ substituting for α and λ from the formulation of the previous section. Using its result we readily write for the Laplace transform of the closing probability: However, in contrast to the general GE distribution, the parameter N ch takes only postive integers. This allows for a simpler expression, i.e. by using the binomial theorem we can write: where we dropped 1 N ch −1−k = 1 inside the sum. After multiplication we get: With this we just have to Laplace transform a sum of exponentials, and we finally obtain:Ψ  To determine the other parameters, we started by the random opening of only one single cluster (a puff), which is always the first step in the stochastic process generating a global Ca 2+ spike. These openings can be modelled by an exponential distribution with puff rate λ 0 [5,13], so that λ 0 is determined by are computed from an appropriate Ca 2+ diffusion profile as in [5]. The Ca 2+ concentration depends on the number of open clusters N o : where c 1 = 0.426 µM is the Ca 2+ concentration resulting from one nearby open cluster. The simple form of Eq. 41 results from the tetrahedral geometry and from the assumption of free boundaries [5]. It can be replaced by a diffusion problem taking more cellular details into account without difficulty. Based on Eq. 41, we obtain all fitting parameters needed to reproduce the minimal stochastic Ca 2+ model analysed in the main text (Tab. S2).

Implementation of the tetrahedron model.
where the prefactors account for the multiplicity of the transitions. Knowing this, we separated the whole stochastic process generating the global Ca 2+ spikes into one inhomogeneous Poisson process describing the puffs, and a Bernoulli trial with success probability C 14 . For the analytically tractable case, the puff process was a homogeneous Poisson process with a constant puff rate λ 0 . We now include a negative feedback to capture phenomena like Ca 2+ store depletion. So we recast the puff process in the following form with λ(t) as time dependent puff rate tending to the asymptotic rate λ 0 and the intensity Λ(t) = t 0 λ(τ )dτ One can now ask for the probability that at a time point t a puff occurs and triggers a global Ca 2+ spike. Conditioned on the ground state S 0 we are looking for the time for the first Ca 2+ spike to occur, the ISI. Formally this leads to: This expression contains the probabilities for all possible puffs occurring before t not leading to a spike, the failed Bernoulli trials with probability (1 − C 14 ).
The very first term is the probability that at t a puff occurs and becomes a global spike and no other puff occurs in [0, t), the second term is the probability for one failed puff in [0, t), the third term handles two failed puffs and so on.

Following permutation and symmetry arguments outlined in the book by Van
Kampen [14] we may write: The resulting stochastic process for the Ca 2+ spiking is also a inhomogeneous Poisson process with the spike rate κ(t) = C 14 λ(t). The new rate is just the puff rate refined with the splitting probability C 14 , this property is known as Poisson splitting. It may be generalised to heterogeneous puff sites with individual splitting probabilities C 1N,i and puff rates λ 0,i , see also main text.

Error analysis for the generic model
The generic model can only be an approximation for the real stochastic process governing the Ca 2+ spike generation. The applied Poisson splitting completely neglects the dynamics of the system between the individual states, instead it implicitly sets all transition times except for the S 0 → S 1 transition to zero.
So we expect to underestimate the true average ISI, T gm av < T tr av . The time the system spends in processing the failed puffs is a major source for the under-estimation. The number of failed puffs is on average given by 1/C 14 , which is for standard parameters ∼ 100. To exactly quantify the approximation error of the generic model we compared its results for ξ → ∞ (no feedback) with the analytic results given by the hierarchic stochastic model. As expected the error grows linear with 1/C 14 (Fig. S1A). We also computed the relative error (T tr av − T gm av )/T tr av as function of IP 3 , and state that it approaches zero for IP 3 → 0 and saturates for IP 3 → ∞ (Fig. S1B). Our main argument for the generic model is the time scale separation between puff and spike events. In accordance with that the approximation error gets bigger with weaker time scale separation, realised by a smaller channel closing rate γ (Fig. S1A-B).