Han-Yu Jiang^1,2, Jun He^,11Department of Physics and Institute of Theoretical Physics, Nanjing Normal University, Nanjing, 210097, China
²Sino-U.S. Center for Grazingland Ecosystem Sustainability/Pratacultural Engineering Laboratory of Gansu Province/ Key Laboratory of Grassland Ecosystem, Ministry of Education/College of Pratacultural Science, Gansu Agricultural University, Lanzhou, 730070, China

Received:2020-09-10Revised:2020-11-2Accepted:2020-11-4Online:2021-01-07


Abstract
Ca²⁺ plays an important role in cell signal transduction. Its intracellular propagation is the most basic process of Ca²⁺ signaling, such as calcium wave and double messenger system. In this work, with both numerical simulation and mean field ansatz, the three-dimensional probability distribution of Ca²⁺, which is read out by phosphorylation, is studied in two scenarios with boundaries. The coverage of distribution of Ca²⁺ is found at an order of magnitude of μm, which is consistent with experimental observed calcium spike and wave. Our results suggest that the double messenger system may occur in the ER-PM junction to acquire great efficiency. The buffer effect of kinase is also discussed by calculating the average position of phosphorylations and free Ca²⁺. The results are helpful to understand the mechanism of Ca²⁺ signaling.
Keywords： Ca ²⁺ signaling;diffusion;Gillespie algorithm;mean field ansatz


PDF (1319KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite
Cite this article
Han-Yu Jiang, Jun He. Three-dimensional cytoplasmic calcium propagation with boundaries. Communications in Theoretical Physics, 2021, 73(1): 015601- doi:10.1088/1572-9494/abc7af

1. Introduction

Cytoplasmic Ca²⁺ is the most basic second messenger of cell signal transduction [1, 2]. It widely involves in the regulation of cell life activities, including respondence to external stimuli, cell membrane permeability, cell secretion, metabolism and differentiation [3]. Within the cellular signaling network, the accurate decoding of diverse Ca²⁺ signal is a fundamental molecular event [4, 5]. In the double messenger system, the activity phospholipase C (PLCβ) induces the hydrolysis of the phospholipid phosphatidylinositol-4, 5-bisphosphate (PIP₂) to generate inositol 1,4,5-trisphosphate (IP₃) and diacylglycerol (DAG) in the plasma membrane (PM). IP₃ mobilizes the endogenous Ca²⁺, transfers the Ca²⁺ stored in endoplasmic reticulum (ER) to the cytoplasm, and increases the intracellular Ca²⁺ concentration. Ca²⁺ combines with protein kinase C (PKC), which is a kind of calcium dependent kinase, and participates in many physiological processes, short-term physiological effects such as cell secretion and muscle contraction, or long-term physiological effects such as cell proliferation and differentiation. PKC is translocated to the inner surface of the PM, which is activated by DAG, another second messenger of phosphatidylinositol signal bound to PM. Then, the serine and threonine residues of different substrate proteins in different cell types are phosphorylated. In the current work, we define the response strength of double messenger system to stimuli as the efficiency of such signal pathway, which is described by the phosphorylation events induced near the position of the stimuli in PM. Ca²⁺ also combines with calmodulin (CaM) to form a Ca²⁺-CaM complex, which then activates the target enzymes [6]. CaM kinases are highly conserved in all eukaryotes, and are important target enzymes. It mediates many functional activities in animal and plant cells, plays an important role in the synthesis of cAMP and cGMP and degradation of glycogen, and smooths muscle contraction and neurotransmitter secretion and synthesis. Furthermore, since CaM may diffuse slowly with a speed smaller than Ca²⁺, it plays a role of buffer of the Ca²⁺ propagation to enlarge its influence coverage [7, 8].

The intracellular propagation and distribution of Ca²⁺ ion affects the activation and operation of downstream components of signaling, such as CaM and PKC [9]. According to the amplitude of the Ca²⁺ spikes, the frequency of the Ca²⁺ waves, and the microdomains of the Ca²⁺ flickers, Ca²⁺ signaling displays different temporal and spatial patterns [10–13]. The distribution of intracellular Ca²⁺ is strictly regionalized, while maintaining the homeostasis of cytoplasmic Ca²⁺ is the prerequisite for normal cell growth [14]. Therefore, it is of great significance to analyze the spatial distribution and dynamic changes of intracellular Ca²⁺ for the study of cell physiological characteristics and signal transduction mechanism. In the past few decades, a variety of experimental methods were developed to measure spatial distribution of Ca²⁺ concentration in cytoplasm, ER, and other organelles. In [15], using digital imaging microscopy and the dye fura-2, distribution of intracellular cytoplasm-free Ca²⁺ were studied. The laser scanning confocal microscope technology was also developed as a research tool for observing the spatial and temporal changes of intracellular calcium distribution [16, 17]. Fluorescent indicators for Ca²⁺ are important tool to detect the Ca²⁺ [18, 19]. By using a laser scanning confocal microscope and the fluorescent calcium indicator fluo-3, Cheng and his collaborators detected the calcium sparks in heart cells [20, 21]. A plasmonic-based electrochemical impedance microscope was also developed to provide valuable information without a fluorescent labeling [22]. On the other hand, the establishment of spatial network model of calcium ion, by biochemical reaction modeling or parameter fitting, is also important to understand the spatial and temporal characteristics of calcium signal and the mechanism of signal transduction. In 2010, Rudiger and Shuai used deterministic and stochastic simulation methods to simulate a network model in which several IP₃ receptors channels release local calcium signals in a single cluster [23]. Qi established a model determined by a simple calcium kinetic equation and a Markov process [24]. The full width at half-maximum of calcium spikes was also analyzed and well reproduced in the simulation [25].

Although a lot of work has been done to study calcium spatial distribution, most of them are about calcium concentration, which is important to understand the calcium spikes and calcium waves, which usually propagate along the membrane. However, for the double messenger system, we need to know the propagation and diffusion of the calcium ion from EM to PM instead of concentration. Hence, the current understanding of the decoding mechanism of CaM and the signal network of calcium PKC still is fragmentary. Besides, little is known about the spatial location of above kinase phosphorylation events for calcium signal transduction study. In [26, 27], the authors developed a frame to estimate the accuracy of the position of external stimuli. In their work, the Ca²⁺ induced by stimuli enters the cytoplasm through a channel in membrane, and is read out by the spatial distribution of phosphorylation events in cytoplasm or membrane. It is an important reaction mechanism of the external stimuli. However, the double messenger system happens between ER and PM, which is quite different from the scenario studied in [26, 27]. Theoretically, the double messenger system should be a diffusion problem with two boundaries. Moreover, with the increase of the distance between ER and PM, the probability distribution of the Ca²⁺ becomes small, which makes the signal propagation in such pathway more difficult. The efficiency of the double messenger system should be dependent on the distance between ER and PM. The study about such issue is scarce in the literature. Hence, it is interesting to preform a study about the propagation of the Ca²⁺ perpendicular to the membrane through the phosphorylation.

By studying the random binding of calcium ions with calmodulin and relevant kinases distributed in different cytoplasmic regions, we can estimate the position of phosphorylation events in the cytoplasm or in the PM, as well as the position of binding of calcium ion to CaM/PKC and the influence domain of calcium ion through phosphorylation. Due to the short-term characteristics of calcium puff and the complexity of intracellular environment, it is difficult and cumbersome to locate the distribution area and binding probability of calcium ions by biological experiments. It is interesting to introduce theoretical model to simulate such process. Based on such new idea, by adopting mathematical simulation and mean-filed ansatz to solve non-linear diffusion equation, we start from the level of single calcium ion to study: (1) the spatial probability of Ca²⁺ read out by phosphorylation from ER to PM, which provides an overall picture of the Ca²⁺ distribution; (2) average positions of phosphorylation and the free Ca²⁺ in a direction perpendicular to the membrane, which is related to buffer effect of CaM; (3) estimation error of the readout position of phosphorylation event, which assesses the effective region of Ca²⁺. Large efficiency of double messenger system benefits from large number of phosphorylation events and small estimation error in PM. With these results, we can analyze the dependence of the efficiency of the double messenger system on the distance between ER and PM. This study provides a theoretical basis for further exploring the mechanism of Ca²⁺ mobilization induced calcium signaling pathway in cell biology, and the correlation between spatial heterogeneity of calcium ions and calcium signaling pathway.

2. Theoretical frame

In the current work, we consider CaM and PKC scenarios (see figure 1). Ca²⁺ enters cytoplasm through calcium pump on ER. In the CaM scenario, the phosphorylation happens in cytoplasm while the PKC attached by Ca²⁺ should be bound onto PM first in PKC scenario. In the followings, we will describe these two scenarios explicitly when presenting master equation and simulation method.

Figure 1.

New window|Download| PPT slide
Figure 1.The CaM and PKC scenarios. The explicit description is presented in text.

2.1. Master equation

First, we present analytical description of Ca²⁺ by the master equation, with which the spatial probability distribution of phosphorylation events is introduced.

2.1.1. Diffusion of free Ca²⁺ until first attachment

In the literature, the permeation of Ca²⁺ ion in water and cytosolic was discussed [28]. Moving ion will be damped because its energy is lost to surrounding medium through viscous resistance to that movement. It can be described by using the Langevin equation [29], $m\tfrac{{\rm{d}}{\boldsymbol{v}}}{{\rm{d}}t}=-m\gamma {\boldsymbol{v}}+{\boldsymbol{R}}$ where the m, v, and γ are the mass, velocity and the friction coefficient of the ion, and R is random force. Average velocity deceases exponentially as ⟨v⟩=⟨v⟩₀e^−γt. Based on the Einstein relation, the friction coefficient of the calcium ion can be evaluated as γ⁻¹≈1 fs [30]. It suggests that initial velocity of calcium ion will be lost in an extremely short time after entering cytoplasm. Hence, permeation of Ca²⁺ behaves as a diffusion.

We define a probability distribution function as ${{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]$ for the state with Ca²⁺ ion unattached to a kinase, and we will call the Ca²⁺ in such state as ‘free Ca²⁺’ here and hereafter. It satisfies initial condition as ${{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t=0]=\delta ({\boldsymbol{r}})$ which suggests that Ca²⁺ ion enters system at r=0 on ER in both scenarios, that is, we define calcium pump as the origin of coordinate. If not losing into cytosolic, the moving Ca²⁺ ion attaches and activates a kinase, CaM or PKC, and it obeys diffusion equation as(1)$\begin{eqnarray}{\partial }_{t}{{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]={D}_{C}{{\rm{\nabla }}}^{2}{{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]-({\nu }_{a}+{\nu }_{l}){{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t],\end{eqnarray}$where the D_C is the diffusion constant for Ca²⁺ in cytosolic, and ν_a and ν_l are the rates of attachment and losing, respectively.

In [26, 27], Ca²⁺ ion enters and kinase binds to the same PM, which is important in the determination of position of external signal. In the current work, we consider three-dimensional PKC scenario (see figure 1), which happens in the important double messenger system. It leads to obvious difference. Ca²⁺ ion enters from ER idealized as an infinite plane and diffuses in three-dimensional space. Furthermore, PKC needs to propagate to PM as another plane at a distance l. Hence, the Ca²⁺ ion will be reflected by boundary, PM and/or ER, which leads to boundary condition as(2)$\begin{eqnarray}{\partial }_{z}{{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]{| }_{z=0}=0,\quad \mathrm{for}\ \mathrm{CaM}\ \mathrm{and}\ \mathrm{PKC},\end{eqnarray}$(3)$\begin{eqnarray}{\partial }_{z}{{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]{| }_{z=l}=0,\quad \mathrm{for}\ \mathrm{PKC}\ \mathrm{only}.\end{eqnarray}$

2.1.2. Attachment and detachment

After the first attachment, a function ${{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]$ is introduced to reflect two probabilities: having kinase activated by Ca²⁺ ion at position r and time t, and having a distribution of phosphorylation event n(ξ) at time t. ${{ \mathcal P }}_{i}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]$ is defined analogously but with Ca²⁺ detached from the kinase. The master equation is given by(4)$\begin{eqnarray}\begin{array}{l}{\partial }_{t}{{ \mathcal P }}_{i}[(n({\boldsymbol{\xi }});{\boldsymbol{r}},t]={D}_{C}{{\rm{\nabla }}}^{2}{{ \mathcal P }}_{i}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]\\ \quad +\,{\nu }_{d}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]-({\nu }_{a}+{\nu }_{l}){{ \mathcal P }}_{i}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t],\end{array}\end{eqnarray}$where ν_d is the rate of detachment.

Because Ca²⁺ ion and active kinase will be also reflected by boundary, we have boundary condition as(5)$\begin{eqnarray}{\partial }_{z}{{ \mathcal P }}_{i,a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]{| }_{z=0}=0,\quad \mathrm{for}\ \mathrm{CaM}\ \mathrm{and}\ \mathrm{PKC},\end{eqnarray}$(6)$\begin{eqnarray}{\partial }_{z}{{ \mathcal P }}_{i}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]{| }_{z=l}=0,\quad \mathrm{for}\ \mathrm{PKC}\ \mathrm{only}.\end{eqnarray}$

2.1.3. Phosphorylation in CaM scenario

Difference appears after activation of kinase in two scenarios. In CaM scenario (see figure 1), the active kinase phosphorylates directly in cytosolic at a rate of ν_p, which can be described as(7)$\begin{eqnarray}\begin{array}{l}{\partial }_{t}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]={D}_{K}{{\rm{\nabla }}}^{2}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]-{\nu }_{d}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]\\ \quad +\,{\nu }_{a}\left\{{{ \mathcal P }}_{i}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]+{{ \mathcal P }}_{i0}[{\boldsymbol{r}},t]\right\}\\ \quad +\,{\nu }_{p}\left\{{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }})-\delta ({\boldsymbol{\xi }}-{\boldsymbol{r}});{\boldsymbol{r}},t]-{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]\right\}.\end{array}\end{eqnarray}$The active kinase diffuses by constant D_K. The Dirac delta function δ(ξ−r) means that a phosphorylation event happens at r, then, the distribution jumps from ${{ \mathcal P }}_{a}[n({\boldsymbol{\xi }})-\delta ({\boldsymbol{\xi }}-{\boldsymbol{r}});{\boldsymbol{r}},t]$ to ${{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]$ at a rate of ν_p.

2.1.4. Binding and phosphorylation in PKC scenario

In PKC scenario, phosphorylation does not happen directly after attachment. The PKC attached by a Ca²⁺ diffuses and obeys a master equation(8)$\begin{eqnarray}\begin{array}{l}{\partial }_{t}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]={D}_{K}{{\rm{\nabla }}}^{2}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]-{\nu }_{d}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]\\ +\,{\nu }_{a}\{{{ \mathcal P }}_{i}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]+{{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]\}.\end{array}\end{eqnarray}$

The kinase should be bound onto PM and activated by DAG to make phosphorylation possible. A new state function ${{ \mathcal P }}_{b}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t]$ is defined as the probability distribution, when Ca²⁺ attaches to a kinase and the complex binds into PM. Because PM is at a fixed distance from ER, we adopt $\bar{{\boldsymbol{r}}}$ to denote the r with z equaling to a distance l. The binding and unbinding of PKC to PM can be written as a boundary condition is [31],(9)$\begin{eqnarray}-{D}_{K}{\partial }_{z}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t]={\nu }_{u}{{ \mathcal P }}_{b}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t]-{\nu }_{b}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t],\end{eqnarray}$where ν_b and ν_u are the rates of binding and unbinding from PM, respectively.

The phosphorylation at a rate of ν_p on PM follows master equation as(10)$\begin{eqnarray}\begin{array}{l}{\partial }_{t}{{ \mathcal P }}_{b}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t]={\nu }_{b}{{ \mathcal P }}_{a}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t]-{\nu }_{u}{{ \mathcal P }}_{b}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t]\\ +\,{\nu }_{p}\{{{ \mathcal P }}_{b}[n(\bar{{\boldsymbol{\xi }}})-\delta (\bar{{\boldsymbol{\xi }}}-\bar{{\boldsymbol{r}}});\bar{{\boldsymbol{r}}},t]-{{ \mathcal P }}_{b}[n({\boldsymbol{\xi }});\bar{{\boldsymbol{r}}},t]\}.\end{array}\end{eqnarray}$

In both PKC and CaM scenarios, we define a probability distribution ${{ \mathcal P }}_{0}[(n({\boldsymbol{\xi }});{\boldsymbol{r}},t]$ when the Ca²⁺ ion is lost, which obeys(11)$\begin{eqnarray}{\partial }_{t}{{ \mathcal P }}_{0}[n({\boldsymbol{\xi }});t]={\nu }_{l}\int {\rm{d}}{\boldsymbol{r}}\left\{{{ \mathcal P }}_{i}[n({\boldsymbol{\xi }});{\boldsymbol{r}},t]+{{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]\right\}.\end{eqnarray}$

In summary, the CaM scenario shown in upper panel of figure 1 is described by diffusion equations in equations (1), (4), (7), (11) and boundary conditions in equations (2), (5), and the PKC scenario shown in lower panel of figure 1 is described by diffusion equations in equations (1), (4), (8), (11) and boundary conditions in equations (2), (3), (5), (6), (9), (10). Hereafter, we will make a dimensionless treatment. The time and space are scaled by ν_p and $\sqrt{{\nu }_{p}/{D}_{C}}$, respectively. The scaled ν_a, ν_d, ν_l, and D_K can be obtained from scaled time and space.

2.2. Stochastic simulation of master equation

With the master equation above, theoretically, we can obtain the distribution of phosphorylation event, n(ξ). However, the master equations are difficult to be solved analytically. In the current work, we adopt the Gillespie algorithm to do the stochastic simulation [26, 27, 32] in three-dimensional space with the boundary conditions.

The simulation follows scenarios shown in figure 1. Ca²⁺ ion enters system at r=0. An event occurs in a step of time Δt=−lnr/(ν_a+ν_l) with r being a random number in a range from 0 to 1. This event is attachment or leaving system, which is determined by drawing a random number from 0 to ν_l+ν_a. If the random number is smaller than ν_l, simulation stops. If not, the position of active kinase is changed in three directions by drawing three independent random numbers Δr from the Gaussian distribution with zero mean and variance 2Δt. If obtained new position is out of the boundaries, it will be reflected by boundaries to a position in allowed region, z>0 for CaM scenario and l>z>0 for PKC scenario.

In CaM scenario, kinase may phosphorylate target protein or detach in a new step of time. The type of event is still determined by drawing a random number. At the same time, the kinase moves randomly. If phosphorylation happens, the position is recorded as ξ_i. If the Ca²⁺ ion detaches from the active kinase, simulation continues to a new loop. When simulation ends, we have a trajectory of n₀ phosphorylation events. With N simulations, we obtain n trajectories with phosphorylation events. We would like to remind that if a simulation does not raise any phosphorylation event, we do not count it into the number of trajectories n.

In PKC scenario, after the kinase is activated, the situation is more complex because it has possibility to be bound onto PM. Let us consider an active kinase in z_old, which certainly locates between two boundaries, goes to a new position z_new after a random movement at time Δ t. If the z_new is still between the two boundaries, we draw a random number r between 0 and 1. If $r\lt {{\rm{e}}}^{-[(l-{z}_{{\rm{new}}})(l-{z}_{{\rm{old}}})]/({D}_{K}{\rm{\Delta }}t)}{\nu }_{b}\sqrt{\pi {\rm{\Delta }}t}/(2\sqrt{{D}_{K}})$, the kinase is bound to the membrane [33, 34]. If the z_new is larger than l, the kinase is bound if $r\gt 1-{\nu }_{b}\sqrt{\pi {\rm{\Delta }}t}/(2\sqrt{{D}_{K}})$. Otherwise, it will be reflected. If the z_new is smaller than 0, kinase is reflected only. If after reflecting, the kinase is still out of boundaries, the above step should be repeated until it goes to a position between the boundaries. After binding to the membrane, a random number is drawn to determine whether phosphorylation happens or kinase unbinds from the membrane. The position of phosphorylation will be recorded as in CaM scenario.

2.3. Mean-field (MF) ansatz

Now we consider MF ansatz where phosphorylation rate at position r is assumed to be proportional to the probability of finding a kinase at this position. In CaM scenario, the expected number $\hat{n}({\boldsymbol{r}})$ of phosphorylation events in the limit t→∞ is $\hat{n}({\boldsymbol{r}})={\bar{p}}_{a}({\boldsymbol{r}})={\int }_{0}^{\infty }{\rm{d}}{t}^{\prime} {p}_{a}({\boldsymbol{r}},t^{\prime} )$, where the barred quantities indicate time-integrated quantities [26, 27]. For equations (1), (4), (7), to obtain the distribution of expected phosphorylation events, we are then left with solving(12)$\begin{eqnarray}\begin{array}{rcl}-\delta ({\boldsymbol{r}}) & = & {{\rm{\nabla }}}^{2}{\bar{p}}_{{i}_{0}}({\boldsymbol{r}})-({\nu }_{a}+{\nu }_{l}){\bar{p}}_{{i}_{0}}({\boldsymbol{r}}),\\ 0 & = & {{\rm{\nabla }}}^{2}{\bar{p}}_{i}({\boldsymbol{r}})+{\nu }_{d}{\bar{p}}_{a}({\boldsymbol{r}})-({\nu }_{a}+{\nu }_{l}){\bar{p}}_{i}({\boldsymbol{r}}),\\ 0 & = & {D}_{K}{{\rm{\nabla }}}^{2}{\bar{p}}_{a}({\boldsymbol{r}})-{\nu }_{d}{\bar{p}}_{i}({\boldsymbol{r}})+{\nu }_{a}[{\bar{P}}_{i}({\boldsymbol{r}})+{\bar{P}}_{i0}({\boldsymbol{r}})].\end{array}\end{eqnarray}$The first equation can be easily solved as(13)$\begin{eqnarray}{\bar{p}}_{{i}_{0}}(r)={{\rm{e}}}^{-\omega r}/2\pi r\equiv {\hat{n}}_{0}({\boldsymbol{r}}),\end{eqnarray}$which is also the probability distribution of last position of free Ca²⁺.

Now, we need solve the two equations left which are coupled to each other. The z is larger than zero for these two equations. However, due to the boundary at z=0 is reflective, it is reasonable to make an even extension. Hence, we can preform the Fourier transformation for three direction as,(14)$\begin{eqnarray}\begin{array}{rcl}0 & = & -{D}_{K}{q}^{2}{\bar{p}}_{a}({\boldsymbol{q}})-{\nu }_{d}{\bar{p}}_{a}({\boldsymbol{q}})+{\nu }_{a}[{\bar{p}}_{i}({\boldsymbol{q}})+{\bar{p}}_{{i}_{0}}({\boldsymbol{q}})],\\ 0 & = & -{q}^{2}{\bar{p}}_{i}({\boldsymbol{q}})+{\nu }_{d}{\bar{p}}_{a}({\boldsymbol{q}})-{\omega }^{2}{\bar{p}}_{i}({\boldsymbol{q}}).\end{array}\end{eqnarray}$After inverse transformation, we have(15)$\begin{eqnarray}\hat{n}({\boldsymbol{r}})={\bar{p}}_{a}(r)=\displaystyle \frac{{\nu }_{a}}{{D}_{K}}\displaystyle \frac{1}{{\lambda }_{+}^{2}-{\lambda }_{-}^{2}}\displaystyle \frac{{{\rm{e}}}^{-{\lambda }_{-}r}-{{\rm{e}}}^{-{\lambda }_{+}r}}{2\pi r},\end{eqnarray}$where ${\lambda }_{\pm }=\sqrt{(a\pm \sqrt{{a}^{2}-4b})/2}$ with a=(ν_a+ν_l)+ ν_d/D_K and b=ν_lν_d/D_K.

As in the CaM scenario, the distribution of phosphorylation events in PKC scenario at time t→∞ is given by $\hat{n}={\bar{p}}_{b}(\bar{{\boldsymbol{r}}})$, which is related to kinase distribution as, ${\partial }_{t}{p}_{b}(\bar{{\boldsymbol{r}}})={\nu }_{b}{p}_{a}(\bar{{\boldsymbol{r}}})-{\nu }_{u}{p}_{b}(\bar{{\boldsymbol{r}}})$ from equation (9). After integrating, we have(16)$\begin{eqnarray}{\hat{n}}_{{\rm{PKC}}}={\bar{p}}_{b}(\bar{{\boldsymbol{r}}})=\displaystyle \frac{{\nu }_{b}}{{\nu }_{u}}{\bar{p}}_{a}(\bar{{\boldsymbol{r}}}).\end{eqnarray}$

One can finds that MF master equation in PKC scenario are the same as equation (12) in CaM scenario. Difference is that PKC scenario has additional boundary condition at z=l. It makes the Fourier transformation adopted in CaM scenario failure to solve these coupled differential equations because it only works in a region from 0 to ∞. Hence, in this work, we do not try to give analytical results, but present the simulation results in PKC scenario. However, for large l, effect of boundary at z=l on master equations becomes small. If neglecting the boundary effect, the ${p}_{a}(\bar{{\boldsymbol{r}}})$ can be solved as in the CaM scenario. Combined with equation (16), one can expect a similar result except a factor for PKC scenario at large l. It can be used to check our simulation results.

3. Results

With above preparation, spatial probability distribution of phosphorylation can be obtained, which will be presented first. Based on such distribution, we will analyze average position and estimation error also.

3.1. Spatial probability distribution of phosphorylation

Because spatial probability distribution of phosphorylation is in three-dimensional space, we will present the results with slice map. In figures 2 and 3, we present the results for CaM scenario. Here, the values of parameters are cited from [35–38], which were measured in experiment.

Figure 2.

New window|Download| PPT slide
Figure 2.Spatial probability distribution of phosphorylation with ν_a=10, ν_d=10, ν_l=20, and D_K=0.02 [35–38] from 10⁷ stochastic simulations in CaM scenario. The panels in the first row are for results of phosphorylation events (Pi) in the mean field ansatz (MF). The panels in the second and fourth rows are for the simulation (Simu) results of phosphorylation events (Pi) and trajectories (Tr), respectively. And comparisons between simulation and mean field ansatz are given in third and fifth rows. The color of the scattering points is determined by Gaussian kernel density estimation. The first three columns are for slices at z=0.0, 0.3, and 0.6, with thickness δz=0.01, and the last column is for x=0 with δx=0.01.


In the first row of figure 2, the distribution with mean field (MF) ansatz in equation (15) is presented. At small z, the phosphorylation concentrates at entrance. When deviating from ER, phosphorylation events distribute in a larger region. In the second and fourth rows, we present simulation results for both phosphorylation events and trajectories. In order to determine event distribution at z, we select events (Pi) at position ξ or trajectories (Tr) at mean position $\hat{{\boldsymbol{\xi }}}$ which are in a range form z−δz/2 to z+δz/2. The simulation results exhibit the same picture as the MF ansatz. To give a more obvious comparison, we present the results accumulated at x direction in the third and fifth rows. The MF results are given by $\hat{n}({\boldsymbol{r}})\delta z\delta {xN}$, which is for events not trajectories. The results suggest that the MF results for events fit the simulation as expected. The phosphorylation events concentrates at the entrance, and diminishes rapidly with the increase of z and x. It is interesting to observe that the results for trajectories are almost the same as these for events. It suggests that with the current parameters adopted, the number of the events in a trajectory is very small, which ensures the effectiveness of MF ansatz.

Different from CaM scenario, all phosphorylations happen on PM in PKC scenario. In figure 3, we present the results with different distances, l=0.05, 0.3 and 0.6 between PM and ER. The simulations can not performed when PM and ER overlaps. Here we adopt a small distance l=0.05 instead of zero. At small distance, phosphorylations concentrate at center as in CaM scenario. With the increase of distance, the numbers of events and trajectories diminish rapidly. As shown in third and fifth rows, the phosphorylation events and trajectories are relatively more than in CaM scenario. As shown in the last column, the number is about two times large than the line for CaM, that is, $\hat{n}({\boldsymbol{r}})\delta {xN}$. At larger distance, the MF lines in CaM scenario fit the PKC simulation very well as expected. Different from the CaM scenario, the number of trajectories is about four times larger than the number of events. In the last rows, we give the scaled results for MF, it is interesting to see that the lines also fit the simulation results very well at large distances. It suggests that the MF ansatz is still established, though a trajectory contains about 4 events averagely.

Figure 3.

New window|Download| PPT slide
Figure 3.Spatial probability distribution of phosphorylation with ν_a=10, ν_d=10, ν_l=20, ν_b=ν_u=1, and D_K=0.02 [35–38] from 2×10⁵ stochastic simulations in PKC scenario. The first three columns are for distances at l=0.05, 0.3, and 0.6, with thickness δz=0.01, and the last column is for x=0 with δx=0.01. Other notations are analogous to those in figure 2.

3.2. Average position of phosphorylation and buffer effect of kinase

As shown in the above, phosphorylation consternates near ER in CaM scenario, and the number of phosphorylation events in PKC scenario is also much larger at small distance than at large distance. To give a numerical description, here we present average position of phosphorylation in CaM scenario in figure 4. Because average positions in x and y directions are zero, we only give ⟨z⟩ here. In mean field ansatz, it is defined as(17)$\begin{eqnarray}\langle z\rangle =\displaystyle \frac{{\int }_{-\infty }^{\infty }{\rm{d}}{x}{\rm{d}}{y}{\int }_{0}^{\infty }{\rm{d}}{z}\,z\,\hat{n}({\boldsymbol{r}})}{{\int }_{-\infty }^{\infty }{\rm{d}}{x}{\rm{d}}{y}{\int }_{0}^{\infty }{\rm{d}}{z}\,\hat{n}({\boldsymbol{r}})}=\displaystyle \frac{{\lambda }_{-}^{2}+{\lambda }_{-}{\lambda }_{+}+{\lambda }_{+}^{2}}{{\lambda }_{-}^{2}{\lambda }_{+}+{\lambda }_{+}{\lambda }_{+}^{2}}.\end{eqnarray}$In simulation, we will collect all positions of events or mean positions of trajectories ξ_zk and the average position $\langle z\rangle ={\sum }_{k=1}^{m}{\xi }_{{zk}}/m$ with m is the total number of events or trajectories.

Figure 4.

New window|Download| PPT slide
Figure 4.Average positions of phosphorylation and free Ca²⁺ in CaM scenario. The results are obtained by varying one of the parameters, ν_a=10, ν_d=10, ν_l=20, and D_K=0.02.


With the parameters adopted in previous subsection, the average ⟨z⟩ is about 0.2, which corresponds to about 3 μm before dimensionless treatment. Such value is consistent with the order of the magnitude of a calcium spark [12]. We also discuss the dependence of ⟨z⟩ on the parameters. The values are very stable with a rate of attachment ν_a<100, a rate of detachment ν_d>1, and a diffusion constant D_K<0.1. The result is more sensitive to the rate of loosing ν_l. With larger ν_a and smaller ν_d, the Ca²⁺ will attached into the CaM more quickly and enduringly, which reflects by larger ⟨z⟩. The large D_K of the kinase also leads to larger ⟨z⟩. Such results suggest that the CaM is buffer of Ca²⁺ diffusion as suggested in [7, 8].

To give the buffer effect of CaM more clearly, we also present the average position of free Ca²⁺ when it attaches to kinase or leaving the system before an attachment as described in equation (13). The average position with MF ansatz can be obtained as(18)$\begin{eqnarray}\langle z\rangle =\displaystyle \frac{{\int }_{-\infty }^{\infty }{\rm{d}}{x}{\rm{d}}{y}{\int }_{0}^{\infty }{\rm{d}}{z}\,z\,{\hat{n}}_{0}({\boldsymbol{r}})}{{\int }_{-\infty }^{\infty }{\rm{d}}{x}{\rm{d}}{y}{\int }_{0}^{\infty }{\rm{d}}{z}\,{\hat{n}}_{0}({\boldsymbol{r}})}=\displaystyle \frac{1}{\omega }.\end{eqnarray}$The simulation results are also presented and compared. We would like to note that the ⟨z⟩ for free Ca²⁺ is only dependent on ν_a and ν_l. In the stable region, the buffer effect is relatively small. When the ⟨z⟩ increases, the buffer effect of the CaM becomes important, and the phosphorylations will happen at a position far away from the ER, for example, with ν_l=1, the average position of phosphorylations is at about 30 μm while the ⟨z⟩ for free Ca²⁺ is only about 5 μm. Larger ν_a and D_K, small ν_d and ν_l are beneficial to the buffer function of CaM.

In PKC scenario, the average position of phosphorylation is fixed at the PM. Here, we only give the results for free Ca²⁺ and compare them with CaM results with MF ansatz as shown in figure 5. Since Ca²⁺ will be reflected by PM and ER, the distance between two membranes will effect the average position for free Ca²⁺. For a large distance l=1.2, as expected the effect of PM becomes small, the results are almost the same as the results in CaM scenario with MF ansatz. In such cases, the buffer effect of the PKC is obvious. Phosphorylations happen at PM which distance is about six time larger than the ⟨z⟩ for free Ca²⁺ with ν_a<1. With decrease of the distance, the ⟨z⟩ becomes small due to the suppression of PM. Meantime, the buffer effect of the kinase becomes relatively small too. For example, at a distance l=0.3, the ⟨z⟩ is about 0.15 with ν_a<1. With increase of both ν_a and ν_l, Ca²⁺ will be attached or lost quickly, so the buffer effect becomes more important.

Figure 5.

New window|Download| PPT slide
Figure 5.Average position of free Ca²⁺ in PKC scenario. The results are obtained by varying one of the parameters, ν_a=10, ν_d=10, ν_l=20, and D_K=0.02.

3.3. Estimation error of readout position of Ca²⁺

In the above, we present the average position in z direction. In x and z directions, average position is zero. Here we will give the estimation error to describe the distribution of phosphorylation parallel to ER. With the mean field ansatz it is written as(19)$\begin{eqnarray}\begin{array}{rcl}{{\ell }}^{2} & = & \displaystyle \frac{{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{x}{\rm{d}}{y}\,({x}^{2}+{y}^{2})\,\hat{n}({\boldsymbol{r}})}{{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{x}{\rm{d}}{y}\,\hat{n}({\boldsymbol{r}})}\\ & = & \displaystyle \frac{2{{\rm{e}}}^{{\lambda }_{-}z}{\lambda }_{-}^{3}(1+{\lambda }_{+}z)-2{{\rm{e}}}^{{\lambda }_{+}z}{\lambda }_{+}^{3}(1+{\lambda }_{-}z)}{{\lambda }_{-}^{2}{\lambda }_{+}^{2}({{\rm{e}}}^{{\lambda }_{-}z}{\lambda }_{-}-{{\rm{e}}}^{{\lambda }_{+}z}{\lambda }_{+})}.\end{array}\end{eqnarray}$In simulation, in order to determine estimation error at z, we will collect the mean position of a simulation $\hat{{\boldsymbol{\xi }}}$ which satifies $z-\delta z/2\lt {\hat{\xi }}_{z}\lt z+\delta z/2$, with a number n_z. The number distribution is defined as $\hat{n}(z)={n}_{z}/(n\delta z)$ and the estimation error as ${{\ell }}^{2}(z)={\sum }_{k=1}^{{n}_{z}}({\hat{\xi }}_{x,k}^{2}+{\hat{\xi }}_{y,k}^{2})/{n}_{z}$. Here, the number distribution $\hat{n}(z)$ is for the phosphorylation trajectories, not for the phosphorylation events.

We present CaM results for the estimation error in figure 6. Dependence of parameters and comparison between MF and simulation results are also presented. For all cases, estimation error increases with the increasing of z. The results are not sensitive to variation of ν_a in a large region from 3 to 300. For ν_a larger than 100, the simulation results begin to deviate from MF results. Dependence on the parameter ν_d is also small, and appears at small value of ν_d. The effect of ν_l is quite large from 3 to 300 while the MF and simulation results fit each other when ν_l larger than 10. The estimation error becomes larger with the increasing of diffusion parameter D_K. Estimation error ℓ² is about 0.2 at the z=0.2 with parameters we adopted, which corresponding 5 μm.

Figure 6.

New window|Download| PPT slide
Figure 6.The variation of estimation error against distance with different parameters from 10⁸ stochastic simulations in CaM scenario. The results are obtained by varying one of the parameters, ν_a=10, ν_d=10, ν_l=20, and D_K=0.02.


In figure 7, We present PKC results for the estimation error, and compare it with the analytical results with equation (19) in the CaM scenario. The results suggest the PKC simulation fit MF results in CaM scenario at distance l>0.2 as expected. Generally speaking, the PKC results exhibit the same picture as the CaM scenario. Besides, estimation error is independent on ν_b and ν_u as expected.

Figure 7.

New window|Download| PPT slide
Figure 7.The variation of estimation error against distance with different parameters from 10⁸ stochastic simulations in the PKC scenario. Except the varied parameter, ν_a=10, ν_d=10, ν_l=20, ν_b=1, ν_u=1, and D_K=0.02.

4. Summary and discussion

Ca²⁺ signal is very important in cell signal transduction, which attracts much attentions from scientific communities in biology, medicine, and physics. For all the research topics of the Ca²⁺ signal, such as calcium spike, calcium wave, and second messenger system, the propagation of Ca²⁺ ion in cytoplasm and its readout are in the most basic level. The study of spatial probability distribution of a Ca²⁺ ion is helpful to understand intracellular Ca²⁺ signal transduction. In the current work, we study three-dimensional Ca²⁺ propagation in two typical CaM and PKC scenarios, and the phosphorylation is adopted to read out Ca²⁺.

In CaM scenario, after released from the calcium channel in ER, Ca²⁺ diffuses in cytoplasm, which corresponds to realistic case, such as a Ca²⁺ event in calcium spike and wave [12]. In such case, Ca²⁺ takes effect near membrane, and effect of other membranes is very small. The result suggests that the probability distribution is only dependent of the radius from entrance. With the parameters chosen [35–38], the phosphorylation events concentrate at the entrance and diminish exponentially with the radius. At the direction perpendicular to the membrane, the average position is only about several μm, which is at the same order of magnitude of the width of the calcium spike. The CaM is considered as buffer of the Ca²⁺ propagation in cytoplasm [7, 8]. With the parameters chosen, the buffer effect is not very large. However, there exist many types of CaM which has different properties. The results suggests that the average position in z direction will increase with the increase of the rate of attachment ν_a and diffusion constant D_K, and decrease of the rate of detachment ν_d and the rate of losing ν_l. At the direction parallel to membrane, the average position is zero, error estimation ℓ² is calculated to show the coverage of influence of a Ca²⁺ ion. At small z, the error estimation is also at an order of several μm, and it will increase with the increase of z almost linearly. Generally speaking, if Ca²⁺ is read out by phosphorylation, coverage of influence of a Ca²⁺ ion is at an order of μm. The rates and diffusion constant of the CaM will effect coverage of influence of Ca²⁺ as a buffer. We would like to note that the results in CaM scenario can be applied to all type of membranes, not only the ER considered here.

The PKC scenario is quite different from the CaM scenario. In such scenario, the PKC attached by Ca²⁺ must bind into the PM, which is the important part of the double messenger system. All phosphorylation events happen in PM, and concentrate at zero point. When the distance between ER and PM becomes smaller, the phosphorylation events will decrease very rapidly, especially at small distance. At the same time, the estimation error becomes smaller also. If we recall that the PKC should bind into the PM and activates by the DAG, a small distance between PM and ER, smaller than 1 μm suggested by current result, will benefit the double message system. In fact, the importance of the ER-PM junction were studied in many works, such as Ca²⁺ influx [39, 40]. The current work suggest that in the ER-PM junction, the efficiency of the double messenger system may be improved greatly.

In summary, the three-dimensional probability distribution of Ca²⁺ which is read out by phosphorylation in cytoplasm is studied in two scenarios. The results in PKC scenario with two boundaries are found consistent with the results in CaM scenario especially at large distance as suggested by MF ansatz. The coverage of distribution of Ca²⁺ is found at an order of μm, which is consistent with experimental observed calcium spike and wave. Our results also suggest that the double messenger system may occurs in the ER-PM junction to acquire great efficiency.

Acknowledgments

This project is partially supported by the National Natural Science Foundation of China under Grant No. 11 675 228 and China postdoctoral Science Foundation under Grant No. 2015M572662XB.

Reference View Option By original order
By published year
By cited within times
By Impact factor

[1]Clapham D E 2017 Calcium signaling
Cell 131 1047
DOI:10.1016/j.cell.2007.11.028 [Cited within: 1]

[2]Berridge M J Bootman M D Lipp P 1998 Calcium-a life and death signal
Nature 395 645
DOI:10.1038/27094 [Cited within: 1]

[3]Berridge M 2016 The inositol trisphosphate/calcium signaling pathway in health and disease
Physiol. Rev. 96 1261
DOI:10.1152/physrev.00006.2016 [Cited within: 1]

[4]Berridge M J 1983 Rapid accumulation of inositol trisphosphate reveals that agonists hydrolyse polyphosphoinositides instead of phosphatidylinositol
Biochem. J. 212 849
DOI:10.1042/bj2120849 [Cited within: 1]

[5]Streb H Irvine R F Berridge M J Schulz I 1983 Release of Cd+from a nonmitochondrial store in pancreatic cells by inositol 1,4, 5-trisphosphate
Nature 306 67
DOI:10.1038/306067a0 [Cited within: 1]

[6]Cheung W 1980 Calmodulin plays a pivotal role in cellular regulation
Science 207 19
DOI:10.1126/science.6243188 [Cited within: 1]

[7]Palecek J Lips M B Keller B U 1999 Calcium dynamics and buffering in motoneurones of the mouse spinal cord
J. Physiol. 520 485
DOI:10.1111/j.1469-7793.1999.00485.x [Cited within: 3]

[8]John L M Mosquera-Caro M Camacho P Lechleiter J D 2001 Control of IP3-mediated Ca²⁺ puffs in Xenopus laevis oocytes by the Ca²⁺-binding protein parvalbumin
J. Physiol. 535 3
DOI:10.1111/j.1469-7793.2001.t01-2-00003.x [Cited within: 3]

[9]Yang D Kim J 2020 Emerging role of transient receptor potential (TRP) channels in cancer progression
BMB Rep. 53 125
DOI:10.5483/BMBRep.2020.53.3.016 [Cited within: 1]

[10]Jaffe L F 1993 Classes and mechanisms of calcium waves
Cell Calcium 14 736
DOI:10.1016/0143-4160(93)90099-R [Cited within: 1]

[11]Wei C Wang X Chen M Ouyang K Song L S Cheng H 2009 Calcium flickers steer cell migration
Nature 457 901
DOI:10.1038/nature07577

[12]Cheng H Lederer W J 2008 Calcium sparks
Physiol. Rev. 88 1491
DOI:10.1152/physrev.00030.2007 [Cited within: 2]

[13]Dani J W Chernjavsky A Smith S J 1992 Neuronal activity triggers calcium waves in hippocampal astrocyte networks
Neuron 8 429
DOI:10.1016/0896-6273(92)90271-E [Cited within: 1]

[14]Bush D S 1995 Calcium regulation in plant cells and its role in signaling
Annu. Rev. Plant. Physiol. Plant. Mol. Biol. 46 95
DOI:10.1146/annurev.pp.46.060195.000523 [Cited within: 1]

[15]Yamashita T Amano H Harada N Su Z L Kumazawa T Tsunoda Y Tashiro Y 1990 Calcium distribution and mobilization during depolarization in single cochlear hair cells: imaging Microscopy and Fura-2
Acta Otolaryngol. 109 256
DOI:10.3109/00016489009107441 [Cited within: 1]

[16]MO’Malley D Wang S 1995 Visualization of calcium activity in nerve cells
IEEE Comput. Graph. Appl. 15 55
DOI:10.1109/38.391491 [Cited within: 1]

[17]Oshima T Ikeda K Furukawa M Ueda N Suzuki H Takasaka T 1996 Distribution of Ca²⁺ channels on cochlear outer hair cells revealed by fluorescent dihydropyridines
Am J. Physiol. 271 C944
DOI:10.1152/ajpcell.1996.271.3.C944 [Cited within: 1]

[18]Miyawaki A et al. 1997 Fluorescent indicators for Ca²⁺ based on green fluorescent proteins and calmodulin
Nature 388 882
DOI:10.1038/42264 [Cited within: 1]

[19]Schäferling M 2012 The art of fluorescence imaging with chemical sensors
Angew. Chem., Int. Ed. 51 3532
DOI:10.1002/anie.201105459 [Cited within: 1]

[20]Cheng H Lederer W J Cannell M B 1993 Calcium sparks: elementary events underlying excitation-contraction coupling in heart muscle
Science 262 740
DOI:10.1126/science.8235594 [Cited within: 1]

[21]Cheng H Lederer M R Lederer W J Cannell M B 1996 Calcium sparks and [Ca²⁺]_i waves in cardiac myocytes
Am. J. Physiol. 270 C148
DOI:10.1152/ajpcell.1996.270.1.C148 [Cited within: 1]

[22]Lu J Li J 2015 Label-free imaging of dynamic and transient calcium signaling in single cells
Angew. Chem., Int. Ed. 54 13576
DOI:10.1002/anie.201505991 [Cited within: 1]

[23]Rüdiger S Shuai J Sokolov I 2010 Law of mass action, detailed balance, and the modeling of calcium puffs
Phys. Rev. Lett. 105 048103
DOI:10.1103/PhysRevLett.105.048103 [Cited within: 1]

[24]Qi H Huang Y D Rüdiger S Shuai J W 2014 Frequency and relative prevalence of calcium blips and puffs in a model of small IP3R clusters
Biophys. J. 106 2353
DOI:10.1016/j.bpj.2014.04.027 [Cited within: 1]

[25]Tan W C Fu C Q Fu C J Xie W J Cheng H P 2007 An anomalous subdiffusion model for calcium spark in cardiac myocytes
Appl. Phys. Lett. 91 183901
DOI:10.1063/1.2805208 [Cited within: 1]

[26]Wasnik V H Lipp P Kruse K 2019 Positional information readout in Ca²⁺ signaling
Phys. Rev. Lett. 123 058102
DOI:10.1103/PhysRevLett.123.058102 [Cited within: 5]

[27]Wasnik V H Lipp P Kruse K 2019 Accuracy of position determination in Ca²⁺ signaling
Phys. Rev. E 100 022401
DOI:10.1103/PhysRevE.100.022401 [Cited within: 5]

[28]Kuyucak S Andersen O S Chung S-H 2001 Models of permeation in ion channels
Rep. Prog. Phys. 64 1427
DOI:10.1088/0034-4885/64/11/202 [Cited within: 1]

[29]Chandrasekhar S 1943 Stochastic problems in physics and astronomy
Rev. Mod. Phys. 43 1
DOI:10.1103/RevModPhys.15.1 [Cited within: 1]

[30]Hille B 1992 Ionic Channels of Excitable Membranes2nd edn Sunderland, MA Sinauer Associates
[Cited within: 1]

[31]Crank J 1975 The Mathematics of Diffusion Oxford Oxford University Press
[Cited within: 1]

[32]Gillespie D T 1977 Exact stochastic simulation of coupled chemical reactions
J. Phys. Chem. 81 25
DOI:10.1021/j100540a008 [Cited within: 1]

[33]Erban R Chapman S J 2007 Stochastic modelling of reaction-diffusion processes: algorithms for bimolecular reactions
Phys. Biol. 4 16
DOI:10.1088/1478-3975/4/1/003 [Cited within: 1]

[34]Andrews S S Bray D 2004 Stochastic simulation of chemical reactions with spatial resolution and single molecule detail
Phys. Biol. 1 137
DOI:10.1088/1478-3967/1/3/001 [Cited within: 1]

[35]Nalefski E A Newton A C 2001 Membrane binding kinetics of protein kinase C betaII mediated by the C2 domain
Biochemistry 40 13216
DOI:10.1021/bi010761u [Cited within: 4]

[36]Smith G D Keizer J E Stern M D Lederer W J Cheng H 1998 Detection in cardiac myocytes
Biophys. J. 75 15
DOI:10.1016/S0006-3495(98)77491-0

[37]Donahue B S Abercrombie R F 1987 Free diffusion coefficient of ionic calcium in cytoplasm
Cell Calcium 8 437
DOI:10.1016/0143-4160(87)90027-3

[38]Schaefer M Albrecht N Hofmann T Gudermann T Schultz G 2001 Diffusion-limited translocation mechanism of protein kinase C isotypes
FASEB J. 15 1634
DOI:10.1096/fj.00-0824fje [Cited within: 4]

[39]Jing J et al. 2015 Proteomic mapping of ER-PM junctions identifies STIMATE as a regulator of Ca²⁺ influx
Nat. Cell Biol. 17 1339
DOI:10.1038/ncb3234 [Cited within: 1]

[40]Prakriya M Lewis R S 2015 Store-operated calcium channels
Physiol. Rev. 95 1383
DOI:10.1152/physrev.00020.2014 [Cited within: 1]