Calcium Signals Driven by Single Channel Noise

Usually, the occurrence of random cell behavior is appointed to small copy numbers of molecules involved in the stochastic process. Recently, we demonstrated for a variety of cell types that intracellular Ca2+ oscillations are sequences of random spikes despite the involvement of many molecules in spike generation. This randomness arises from the stochastic state transitions of individual Ca2+ release channels and does not average out due to the existence of steep concentration gradients. The system is hierarchical due to the structural levels channel - channel cluster - cell and a corresponding strength of coupling. Concentration gradients introduce microdomains which couple channels of a cluster strongly. But they couple clusters only weakly; too weak to establish deterministic behavior on cell level. Here, we present a multi-scale modelling concept for stochastic hierarchical systems. It simulates active molecules individually as Markov chains and their coupling by deterministic diffusion. Thus, we are able to follow the consequences of random single molecule state changes up to the signal on cell level. To demonstrate the potential of the method, we simulate a variety of experiments. Comparisons of simulated and experimental data of spontaneous oscillations in astrocytes emphasize the role of spatial concentration gradients in Ca2+ signalling. Analysis of extensive simulations indicates that frequency encoding described by the relation between average and standard deviation of interspike intervals is surprisingly robust. This robustness is a property of the random spiking mechanism and not a result of control.


List of Figures
1 IP 3 R model (Supporting Figure 1 and Table 1) We use a modified version of the DeYoung-Keizer model (DKM), but the following derived reaction diffusion system and its solution can also be used with other channel models. The DKM assumes each IP 3 R to consist of four identical subunits with 3 binding sites each.
One binding site for IP 3 , one for Ca 2+ activating the subunit and another one for Ca 2+ , which inhibits the subunit. The two binding sites for Ca 2+ with a higher affinity for the activating site and a lower affinity for the dominant inhibiting site is a minimal choice to generate nonlinearities which are essential for Ca 2+ induced Ca 2+ release.
Since each of the 3 binding sites can be free or occupied, a single subunit has 2 3 different states X ijk and 12 possible transitions, which can be visualized on a cube as shown in Supporting Figure 1A.  Table 1 for values of rates b i and rate constants a i . B: Stationary open probability P o defined by Equation (2) for the rate constants a i and rates b i given in Supporting Table 1.
The first index of X ijk specifies IP 3 binding and is 1 if IP 3 is bound and 0 otherwise.
Analogously, the second index indicates Ca 2+ binding to the activating site, and the last one corresponds to Ca 2+ binding to the dominant inhibiting site. A subunit is active in the a 1 20 (µMs) −1 IP 3 binding with no inhibiting Ca 2+ bound b 1 20 s −1 IP 3 dissociation with no inhibiting Ca 2+ bound a 2 0.001 (µMs) −1 Ca 2+ binding to inhibiting site with IP 3 bound b 2 0.03 s −1 Ca 2+ dissociation from inhib. site with IP 3 bound a 3 2.6 (µMs) −1 IP 3 binding with inhibiting Ca 2+ bound b 3 20 s −1 IP 3 dissociation with inhibiting Ca 2+ bound a 4 0.025 (µMs) −1 Ca 2+ binding to inhib. site with no IP 3 bound b 4 0.1 s −1 Ca 2+ dissociation from inhib. site without IP 3 a 5 10 (µMs) −1 Ca 2+ binding to activating site b 5 1.225 s −1 Ca 2+ dissociation from activating site Supporting where γ −1 With these relations we can determine the stationary open probability P o in dependence on Ca 2+ and IP 3 . Since a channel opens if three or four subunits are in the state X 110 , the open probability takes the form which is shown in Supporting Figure 1B in dependence on Ca 2+ and IP 3 for the rates given in Supporting 2 Derivation of the linearized system (Supporting Figure 2 and Table 2) To reflect the most important cytosolic properties, we take free [Ca 2+ ] one mobile [B] and one immobile buffer [B i ] in the cytosol into account leading to the following reaction diffusion system where we applied the buffer conservation and used a linear pump and leak flux with the flux constants P p and P l . J j (t) is the stochastically in time varying channel cluster current of the jth cluster, k + and k − denote the capture and dissociation rates of the buffers, and We introduce dimensionless concentrations by where K B denotes the dissociation constant of the mobile buffer. For one channel cluster, point leads to a quantitatively well approximation over a large concentration range. For regions, where the concentration is even higher than 2 K D , i.e. close to open channels, the dynamics is mostly determined by the diffusion and channel terms.
this leads to To obtain a dimensionless system we choose defining the diffusion length L and reaction time T . They are used to rescale time t → t/T and space r → r/L. Similar the remaining quantities in equations (5) dissociation constants ratio of the mobile and immobile buffer Supporting Table 2: Definition of non-dimensional parameters as given in Table 1 in the main text.
in dimensionless parameters given in Supporting Table 2 and in Table 2 of the main text.
Thus equations (5) take the form In order to solve these equations analytically we linearize them around the initial state which is assumed to be the stationary state when all channels are closed. The scaled initial conditions read denotes the ratio of the dissociation constants of the mobile and immobile buffer respectively. Replacing c = c 0 + δc,ē =ē 0 + δē, b = b 0 + δb and b i = b i,0 + δb i in equations (7) we get We now neglect all nonlinear terms in (9) in δc, δb and δb i and end up with the linearized dimensionless system of the form For more convenient reading we drop the δs of the concentrations and find the system given in the main text with the definitions of σ m = (1 + c 0 ), σ im = R (1 + κc 0 ) and 3 Deriving the solution (Supporting Figure 3) The model uses Green's functions [3,4] to determine the concentrations. The solution of a partial differential equation calculated by Green's function is where V denotes the cell volume and S the cell surface. The solution depends on the initial concentration distribution F 0 (r ), the time dependent boundary condition Φ(r , τ ) and the volume production term F (r , τ ).
Before we solve the system (11) by coupled Green's functions for a spherical cell with radius R, we first determine the Green's function for a single component system, i.e.
neglect buffer reactions and coupling with the ER.
The general equation in spherical coordinates reads where g(r, θ, t) is a source density depending on r and θ and t and f (r, θ) denotes the initial condition. The boundary conditions given by j(θ, t) for no-flux conditions (13c) specifies influx through the cell membrane or is given in case of Dirichlet boundary conditions (13d) by the concentration at the surface.
The Green's function is the response of a system at point P (r) at time t due to a δ source in time (t ) and space at point P (r ). Hence, for two points we can use the symmetry and neglect first the φ dependence which can be incorporated later by trigonometric properties.
The corresponding equation of Equation (13) for the Green's function G = G(r, θ, t|r , θ , t ) takes the form where G has to fulfil the corresponding boundary condition in (13). This problem can be solved by Laplace transformation and separation ansatz. After Laplace transform with respect to t the governing equation of the transformed Green's functionG(r, θ, s|r , θ , t ) We first solve the homogeneous problem being the Helmholtz equation where the λs are determined by the particular boundary condition.
It reads in spherical coordinates which can be solved by a standard separation ansatz. We expand the space dependent part D∇ 2G in eigenfunctions of the Laplace operator ∇ 2 . The radial part leads to Bessel's differential equation and the angle dependent part obeys a Legendre differential equations.
Due to convergence restrictions, the solution of the Helmholtz equation (16) takes the form ψ lp (r, θ) = J l+1/2 (λ lp r) where J l+1/2 (x) denotes the Bessel function of the first kind, P l (cos θ) is the Legendre polynomial and λ lp is determined for no-flux boundary condition by ∂ ∂r J l+1/2 (λ lp r) Thus we can solve Equation (15) by inserting the ansatz leading to By applying the integral-operators we get due to orthogonality where the norms N are given by Equation (24) determines the unknown coefficients β lp . The solution in Laplace space is thus given bỹ It can be transformed back easily to time by the residual theorem, since we have first order poles, s + Dλ 2 , along the negative real axis only. The Green's function of the inhomogeneous diffusion problem (13) without the φ dependence finally is For simulation of a cell the spherical symmetry is not valid and we have to introduce explicitly the φ dependence in the Green's function (27). This only depends on the cosines of the angles of the two points P (r) and P (r ), and hence we can rotate the coordinate system such as one of the angles is zero leading to P l (cos θ) = 1. The angle Θ between the points is given by as shown Supporting Figure 3.
The Green's function takes the form where the φ dependence gives another normalization factor of 1/(2π).  Figure 3: Sketch of the angles of two points P and P in spherical coordinates. For more than two asymmetrically distributed points we have to incorporate the φ dependence. Since the solution (27) does only depend on the cosine of the angle between the two points P and P , we can use Equation (28) to rotate the system with respect to the coordinates of the two points.

Ca 2+ dynamics
We now introduce the coupling with buffer reactions. Therefore we write the RDS (11) in matrix form a 11 a 12 a 13 a 21 a 22 0 In order to solve this system of coupled PDEs by coupled Green's functions or a Green's dyadic G [5], we have to solve similar to Equation (14) the following problem, g 11 g 12 g 13 g 21 g 22 g 23 Analogously to (15) the time derivative can be replaced by s due to a Laplace transform leading toLG With the same boundary conditions for calcium and the buffers, the system (30) can be solved by the spectral ansatzG where ψ lp (r, θ) is the solution for the Helmholtz equation (18), which respects the appropriate boundary conditions. Thus, the Green's matrix is determined by the amplitude matrix α lp . By inserting (33) into equation (32) we get ∞ l=0,p=0 By applying the integral operators (23) on both sides, the amplitude matrix is given by with the coupling matrix with the previously introduced shortcuts σ m = (1 + c 0 ), σ im = R (1 + κc 0 ) and σ c = (34) can be transformed back into real space by the property of the matrix inversion the concentrations at point r due to release at point r c and uptake are given by with the relation (28) for P l (cos Θ) (See Supporting Figure 3). The time and boundary dependent response functions read

Boundary conditions
Intracellular Ca 2+ dynamics can be determined by Neumann (no-flux) boundary condition or by Dirichlet boundary condition in dependence on the cell type and kind of experiment.
The boundary conditions at the plasma membrane r = b c (the scaled cell radius R) are reflected by the λ lp s and the two last terms in Equation (40). No-flux boundary conditions require ∂ ∂r J l+1/2 (λ lp r)

Average concentrations
The channel currents are set by the difference of the average Ca 2+ concentrations in the cytosol and in the ER (see Equation 1 in main paper). These concentrations can be determined by spatial integration over the whole cell.
If we do not assume Ca 2+ entry through the cell membrane, the total amount of Ca 2+ will stay constant With the assumption that at time t = 0 the cell has no open channels and is in equilibrium we have also the relation where γ is the volume ratio between the cytosol and the ER. Hence, to calculate the average Ca 2+ concentration within the ERē(t), we rely on the average concentrationsc of all three components. Since V cytc = V cell dV c we have to integrate the solution (40) over the entire cell:c The φ integration simply gives a factor of 2π whereas the other two integrations lead to the equations R = bc 0 J l+1/2 (λ lp r) where 1 F 2 [x, y, z] denotes the hyper-geometric function [6]. Hence, only modes with l = 0 contribute to the global concentration and Equation (48a) can be simplified to With this solution we can calculate the total cell response for one cluster by what can be used to determineē(t) by relations (45) and (46). In the case of no Ca 2+ conservation,ē(t) is given bȳ this means by the difference of the initial ER concentrationē 0 and the difference of the released Ca 2+ and Ca 2+ pumped back into the ER. This method requires the calculation of the average cytosolic Ca 2+ concentration as well.
The real concentrations are obtained by rescalingc(t) andē(t) according to Equation (7).

Green's cell algorithm implementation
The analytical solution for the concentration dynamics (40)

Gillespie algorithm
The Gillespie algorithm allows for simulation of stochastic processes [9]. Given the actual time t, the probability that the next stochastic event occurs in the infinitesimal time interval [t + τ, t + τ + dt] and is an event Ξ i , is given by where α 0 = α j is the sum of all propensities. The event probability P (τ, i) can be realized by drawing two random numbers r 1 and r 2 from a uniform distribution in the interval [0, 1]. Then τ and i are determined by The original Gillespie method assumes that the propensities during two events stay constant. This is not valid for our problem, since opening and closing of channels change the Ca 2+ concentration respectively the propensities by up to three orders of magnitudes.
To resolve this problem we use the method described in [8], which adopt the hybrid version of the Gillespie algorithm [7] to Ca 2+ dynamics. The time of the next stochastic event is determined by solving  is that we only have to calculate the concentration at points where IP 3 R clusters are localized. Another advantage is the linearity, which enables us to calculate the concentrations of clusters independently and to superpose them in the following.
These properties can be used for a sufficiently parallelized algorithm depicted in Supporting Figure 4. Therefore the problem is split into two parts leading to two different kinds of processes. The master process is the "coordinator" of the algorithm, which determines the channel transitions, reaction times and collects global properties. For this tasks it requires the concentration at the IP 3 R clusters. These are determined by worker processes, which calculate the concentration at a cluster according to Equation (40).
The Ca 2+ concentration of each cluster is calculated by a single worker process for which the channel state history of all clusters is required. Therefore each worker process has a copy of the source term σ j (t) of each cluster. These source • If g new < ln r 1 (no stochastic event occurs), the master sets t old = t new and determines the next time step τ , which is broadcasted to the worker processes.   Table 1 of the main text.
The prototype of a cell was simulated with different values of the cytosolic Ca 2+ base level and IP 3 concentration leading to distinct channel behaviors as shown in Supporting If both concentrations are high as in panel A, the channels do not exhibit a cooperative signal, since due to high Ca 2+ 0 most channels are inhibited, and as soon as they are excitable, they will open again and return into the inhibited state leading to the shown noisy signal. This mechanism also holds for very low [IP 3 ] (D), but now the total amount of sensitized, inhibited and open channels is decreased.
If we switch Ca 2+ 0 to physiological concentrations at high [IP 3 ], i.e. going from A to B, we observe very regular oscillations. Due to the high IP 3 concentration, most channels are in the excitable state, and as soon as a single IP 3 R opens, a global wave travels through the system, synchronizing the inhibition of the channels and terminating the Ca 2+ release.
As soon as the inhibiting Ca 2+ dissociates from the corresponding binding site, one channel will open again since Ca 2+ 0 respectively the open probability is high. For a further decrease of Ca 2+ 0 , means going to B and F, the oscillations become slower and more irregular, as the probability of an initial event decreases (note the different time ranges for each column in Supporting Figure 5). In both cases (B and F) only rare single events, as blips or puffs, are observed, caused by full inhibition and additionally by low Ca 2+ 0 in panel F. Lowering the IP 3 concentration for fixed Ca 2+ 0 , i.e. going from B to C, H and E or from F to G, K and J, causes an increase of T av and shrinks the amplitudes, as the channels are less sensitized and the nucleation probability decreases. In these less sensitized  Analogously, we determined also individual slopes for cells which do not vary in the buffer concentration but in the IP 3 concentration (between 60 nM and 0.3 µM), spatial arrangements (minimal intercluster distances between 1 µm and 2.5 µm) and pump strength (between 22 s −1 and 250 s −1 ). It turns out that the observed dependence of the slope on the current strength is stable under these circumstances as well. This points to a functional robustness as explained in the main text.