1. Introduction
Cytoplasmic Ca2+ 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 Ca2+ 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 (PIP2) to generate inositol 1,4,5-trisphosphate (IP3) and diacylglycerol (DAG) in the plasma membrane (PM). IP3 mobilizes the endogenous Ca2+, transfers the Ca2+ stored in endoplasmic reticulum (ER) to the cytoplasm, and increases the intracellular Ca2+ concentration. Ca2+ 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. Ca2+ also combines with calmodulin (CaM) to form a Ca2+-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 Ca2+, it plays a role of buffer of the Ca2+ propagation to enlarge its influence coverage [7, 8].
The intracellular propagation and distribution of Ca2+ ion affects the activation and operation of downstream components of signaling, such as CaM and PKC [9]. According to the amplitude of the Ca2+ spikes, the frequency of the Ca2+ waves, and the microdomains of the Ca2+ flickers, Ca2+ signaling displays different temporal and spatial patterns [10–13]. The distribution of intracellular Ca2+ is strictly regionalized, while maintaining the homeostasis of cytoplasmic Ca2+ is the prerequisite for normal cell growth [14]. Therefore, it is of great significance to analyze the spatial distribution and dynamic changes of intracellular Ca2+ 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 Ca2+ concentration in cytoplasm, ER, and other organelles. In [15], using digital imaging microscopy and the dye fura-2, distribution of intracellular cytoplasm-free Ca2+ 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 Ca2+ are important tool to detect the Ca2+ [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 IP3 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ read out by phosphorylation from ER to PM, which provides an overall picture of the Ca2+ distribution; (2) average positions of phosphorylation and the free Ca2+ 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 Ca2+. 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 Ca2+ 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). Ca2+ enters cytoplasm through calcium pump on ER. In the CaM scenario, the phosphorylation happens in cytoplasm while the PKC attached by Ca2+ 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. The CaM and PKC scenarios. The explicit description is presented in text. |
2.1. Master equation
First, we present analytical description of Ca2+ by the master equation, with which the spatial probability distribution of phosphorylation events is introduced.
2.1.1. Diffusion of free Ca2+ until first attachment
In the literature, the permeation of Ca2+ 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⟩0e−γt. Based on the Einstein relation, the friction coefficient of the calcium ion can be evaluated as γ−1 ≈ 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 Ca2+ behaves as a diffusion.
We define a probability distribution function as ${{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t]$ for the state with Ca2+ ion unattached to a kinase, and we will call the Ca2+ in such state as ‘free Ca2+’ here and hereafter. It satisfies initial condition as ${{ \mathcal P }}_{{i}_{0}}[{\boldsymbol{r}},t=0]=\delta ({\boldsymbol{r}})$ which suggests that Ca2+ 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 Ca2+ ion attaches and activates a kinase, CaM or PKC, and it obeys diffusion equation as
$\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 DC is the diffusion constant for Ca2+ in cytosolic, and νa and νl are the rates of attachment and losing, respectively.In [26, 27], Ca2+ 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. Ca2+ 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 Ca2+ ion will be reflected by boundary, PM and/or ER, which leads to boundary condition as
$\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}$
$\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 Ca2+ 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 Ca2+ detached from the kinase. The master equation is given by
$\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 Ca2+ ion and active kinase will be also reflected by boundary, we have boundary condition as
$\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}$
$\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
$\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 DK. 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 Ca2+ diffuses and obeys a master equation
$\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 Ca2+ 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],
$\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
$\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 Ca2+ ion is lost, which obeys
$\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 DK 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. Ca2+ 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 Ca2+ ion detaches from the active kinase, simulation continues to a new loop. When simulation ends, we have a trajectory of n0 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 zold, which certainly locates between two boundaries, goes to a new position znew after a random movement at time Δ t. If the znew 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 znew 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 znew 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
$\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 $\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 Ca2+.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,
$\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 $\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/DK and b = νlνd/DK.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
$\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
Figure 2. Spatial probability distribution of phosphorylation with νa = 10, νd = 10, νl = 20, and DK = 0.02 [35–38] from 107 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. Spatial probability distribution of phosphorylation with νa = 10, νd = 10, νl = 20, νb = νu = 1, and DK = 0.02 [35–38] from 2 × 105 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
$\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. Average positions of phosphorylation and free Ca2+ in CaM scenario. The results are obtained by varying one of the parameters, νa = 10, νd = 10, νl = 20, and DK = 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 DK < 0.1. The result is more sensitive to the rate of loosing νl. With larger νa and smaller νd, the Ca2+ will attached into the CaM more quickly and enduringly, which reflects by larger ⟨z⟩. The large DK of the kinase also leads to larger ⟨z⟩. Such results suggest that the CaM is buffer of Ca2+ diffusion as suggested in [7, 8].
To give the buffer effect of CaM more clearly, we also present the average position of free Ca2+ 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
$\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 Ca2+ 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 Ca2+ is only about 5 μm. Larger νa and DK, 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 Ca2+ and compare them with CaM results with MF ansatz as shown in figure 5. Since Ca2+ will be reflected by PM and ER, the distance between two membranes will effect the average position for free Ca2+. 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 Ca2+ 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, Ca2+ will be attached or lost quickly, so the buffer effect becomes more important.
Figure 5. Average position of free Ca2+ in PKC scenario. The results are obtained by varying one of the parameters, νa = 10, νd = 10, νl = 20, and DK = 0.02. |
3.3. Estimation error of readout position of Ca2+
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
$\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 nz. 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 DK. Estimation error ℓ2 is about 0.2 at the z = 0.2 with parameters we adopted, which corresponding 5 μm.
Figure 6. The variation of estimation error against distance with different parameters from 108 stochastic simulations in CaM scenario. The results are obtained by varying one of the parameters, νa = 10, νd = 10, νl = 20, and DK = 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. The variation of estimation error against distance with different parameters from 108 stochastic simulations in the PKC scenario. Except the varied parameter, νa = 10, νd = 10, νl = 20, νb = 1, νu = 1, and DK = 0.02. |
4. Summary and discussion
Ca2+ 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 Ca2+ signal, such as calcium spike, calcium wave, and second messenger system, the propagation of Ca2+ ion in cytoplasm and its readout are in the most basic level. The study of spatial probability distribution of a Ca2+ ion is helpful to understand intracellular Ca2+ signal transduction. In the current work, we study three-dimensional Ca2+ propagation in two typical CaM and PKC scenarios, and the phosphorylation is adopted to read out Ca2+.
In CaM scenario, after released from the calcium channel in ER, Ca2+ diffuses in cytoplasm, which corresponds to realistic case, such as a Ca2+ event in calcium spike and wave [12]. In such case, Ca2+ 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 Ca2+ 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 DK, 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 ℓ2 is calculated to show the coverage of influence of a Ca2+ 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 Ca2+ is read out by phosphorylation, coverage of influence of a Ca2+ ion is at an order of μm. The rates and diffusion constant of the CaM will effect coverage of influence of Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ 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.