A Dynamic Model of Cytosolic Calcium Concentration Oscillations in Mast Cells

Abstract

In this paper, a dynamic model of cytosolic calcium concentration (${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$) oscillations is established for mast cells (MCs). This model includes the cytoplasm (Cyt), endoplasmic reticulum (ER), mitochondria (Mt), and functional region (μd), formed by the ER and Mt, also with ${\mathrm{Ca}}^{2+}$ channels in these cellular compartments. By this model, we calculate ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations that are driven by distinct mechanisms at varying $k_{\mathrm{deg}}$ (degradation coefficient of inositol 1,4,5-trisphosphate, ${\mathrm{IP}}_3$ and production coefficient of ${\mathrm{IP}}_3$), as well as at different distances between the ER and Mt (ER–Mt distance). The model predicts that (i) Mt and μd compartments can reduce the amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations, and cause the ER to release less ${\mathrm{Ca}}^{2+}$ during oscillations; (ii) with increasing cytosolic ${\mathrm{IP}}_3$ concentration (${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$), the amplitude of oscillations increases (from 0.1 μM to several μM), but the frequency decreases; (iii) the frequency of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations decreases as the ER–Mt distance increases. What is more, when the ER–Mt distance is greater than 65 nm, the μd compartment has less effect on ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations. These results suggest that Mt, μd, and ${\mathrm{IP}}_3$ can all affect the amplitude and frequency of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations, but the mechanism is different. The model provides a comprehensive mechanism for predicting cytosolic ${\mathrm{Ca}}^{2+}$ concentration oscillations in mast cells, and a theoretical basis for calcium oscillations observed in mast cells, so as to better understand the regulation mechanism of calcium signaling in mast cells.

1. Introduction

In recent years, more and more studies have found that mast cells (MCs) play a major role in the mechanism of acupuncture effect, and substances such as histamine and leukotriene, secreted by mast cells in the process of acupuncture, may be the key factors affecting acupuncture [1]. As an important second messenger, calcium signaling widely exists in various cell physiological processes, participating in the regulation of neurotransmitters released by neurons and astrocytes, metabolic processes, cell maturation, differentiation, and death [2,3,4,5]. The increase in cytosolic ${\mathrm{Ca}}^{2+}$ concentration (${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$) can be divided into the following two pathways: (i) the release of ${\mathrm{Ca}}^{2+}$ from the intracellular ${\mathrm{Ca}}^{2+}$ stores, mainly the endoplasmic reticulum (ER, the largest ${\mathrm{Ca}}^{2+}$ store), or (ii) extracellular ${\mathrm{Ca}}^{2+}$ influx to cytosol (Cyt), through the opening of plasma membrane ${\mathrm{Ca}}^{2+}$ channels. It has been widely accepted that ${\mathrm{Ca}}^{2+}$ release-activated ${\mathrm{Ca}}^{2+}$ (CARC) channels are the main mode of ${\mathrm{Ca}}^{2+}$ influx in electrically non-excitable cells, including MCs [6]. Meanwhile, it is also known that ER calcium depletion activates CRAC channels on the plasma membrane, leading to extracellular ${\mathrm{Ca}}^{2+}$ influx and endoplasmic reticulum ${\mathrm{Ca}}^{2+}$ supplementation. ${\mathrm{IP}}_3$ interacts with Ca²⁺ channels in the ER, causing the release of stored ${\mathrm{Ca}}^{2+}$, and the depletion of ${\mathrm{Ca}}^{2+}$ in the ER triggers ${\mathrm{Ca}}^{2+}$ entry through CRAC channels. CRAC channels in MCs are non-voltage-gated and show a characteristic inward rectification [7]. Therefore, the ${\mathrm{Ca}}^{2+}$ flow of CRAC channels is related to the ${\mathrm{Ca}}^{2+}$ concentration of the ER and Cyt. In contrast to the CRAC channels, the plasma membrane ${\mathrm{Ca}}^{2+}$-ATPase (PMCA) channels extrude ${\mathrm{Ca}}^{2+}$ to the extracellular space, to maintain calcium concentration balance. They can transfer ${\mathrm{Ca}}^{2+}$ against the concentration gradient in the presence of ATP [8]. Previously, we believed that the mitochondria (Mt, the second largest calcium store) only play a role in regulating cytosolic calcium concentration under the pathological condition of very high intracellular calcium ion concentration [9]. Until the 1990s, there were studies demonstrating that the non-pathological increase of cytosolic calcium concentration was accompanied by ${\mathrm{Ca}}^{2+}$ concentration increasing in the mitochondrial matrix (${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}$) [10]. Moreover, recent studies have shown that the Mt regulate ${\mathrm{Ca}}^{2+}$ oscillations by firstly uptaking ${\mathrm{Ca}}^{2+}$, and subsequently releasing it [11,12,13,14]. The evidence was supported by the discovery of the functional region formed by the ER and Mt (μd) [15,16,17]. The functional region is composed of the Mt membrane, ER membrane, and Cyt between them. However, the assumption distances of μd vary hugely, from less than 10 nm to more than 200 nm [18]. When mast cells are mechanically stimulated, the intracellular ${\mathrm{Ca}}^{2+}$ concentration will increase, then leukotriene ${\mathrm{C}}_4$ (${\mathrm{LTC}}_4$) will be produced, which can activate phospholipase C ($\mathrm{PLC}$), to promote ${\mathrm{PIP}}_2$ decomposition to ${\mathrm{IP}}_3$ [19]. Subsequently, ${\mathrm{IP}}_3$ binds to the inositol 1,4,5-trisphosphate receptor (${\mathrm{IP}}_3\mathrm{R}$) on the ER membrane, which induces ${\mathrm{Ca}}^{2+}$ release from the ER. ${\mathrm{IP}}_3\mathrm{R}$ is regulated by ${\mathrm{Ca}}^{2+}$ in a biphasic manner (stimulatory at low levels/inhibitory at high levels) [19,20,21,22]. The mitochondrial ${\mathrm{Ca}}^{2+}$ uniporter (MCU) in mitochondria uptakes ${\mathrm{Ca}}^{2+}$ quickly when exposed to high ${\mathrm{Ca}}^{2+}$ concentration environments around the opened ${\mathrm{IP}}_3\mathrm{R}$ channel pore. MCs are activated by mechanical stimulation, and the ${\mathrm{Ca}}^{2+}$ concentration in μd can reach more than 10 times that in the cytoplasm, which was enough to activate the mitochondrial MCU channel and allow ${\mathrm{Ca}}^{2+}$ uptake [23,24,25]. This explains why high ${\mathrm{Ca}}^{2+}$ is observed when global ${\mathrm{Ca}}^{2+}$ is low.

In order to better explore the physiological mechanism of calcium oscillations, researchers have carried out a large number of experiments. For example, Joseph Di Capite et al. [26] recorded calcium waves in mast cells, and studied the effects of CARC channels on them. Osipchuk et al. [27] studied the effect of ATP on calcium signaling spread between mast cells, and recorded the calcium signaling. At the same time, mathematical modeling can establish the internal relationship between experimental data and parameters, and predict the possible phenomena, so as to save time and cost. In the early years, Goldbeter et al. [28], Hofer [29], and Li and Rinzel [30] described calcium oscillations that only consider the function of the endoplasmic reticulum, while little consider the influences of the mitochondria. Based on the ${\mathrm{Ca}}^{2+}$ dynamic model proposed by Othmer-Tang et al., Falcke et al. [31] added the mitochondrial ${\mathrm{Ca}}^{2+}$ cycle equation into the model. Shi [32] and Qi et al. [33] established a theoretical model, considering the influences of the mitochondria. They explored the effect of the interaction between the mitochondria and the endoplasmic reticulum on calcium oscillations. Arash Moshkforoush et al. [34] developed a compartmental closed-cell mathematical model of ${\mathrm{Ca}}^{2+}$ dynamics that includes a functional region between the ER and Mt. However, they do not consider the effect of plasma membrane calcium channels and the extracellular ${\mathrm{Ca}}^{2+}$ concentration on intracellular calcium oscillations. Although there are many mathematical models that describe calcium oscillations, most of them only consider the effects of the endoplasmic reticulum. Therefore, in order to explain the calcium signaling observed in mast cells, and explore the influence of each compartment on oscillations more accurately, a comprehensive dynamic model of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations is established in this paper. This model takes the following cellular compartments into account: plasma membrane (Mem), cytoplasm (Cyt), endoplasmic reticulum (ER), and mitochondria (Mt). The major ${\mathrm{Ca}}^{2+}$ channels and ${\mathrm{Ca}}^{2+}$ buffering in these compartments are considered. The functional region formed by the ER and Mt (μd) is explicitly assumed as a ${\mathrm{Ca}}^{2+}$ pool. The degradation and production of ${\mathrm{IP}}_3$ is added to this model, to investigate the effect of ${\mathrm{IP}}_3$ dynamics on ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations.

2. Mathematical Model

The full model includes plasma membrane channels, and degradation and production of ${\mathrm{IP}}_3$, Cyt, ER, Mt, and μd. MCs will release ${\mathrm{IP}}_3$ after they are activated by mechanical stimuli. Then, ${\mathrm{IP}}_3$ bines to ${\mathrm{IP}}_3\mathrm{R}$ to trigger the intracellular ${\mathrm{Ca}}^{2+}$ signal. The whole progress is shown in Figure 1. Calcium dynamics in each compartment are governed by a balance of ${\mathrm{Ca}}^{2+}$ fluxes, leaks, and buffering processes.

2.1. Cross-Membrane ${\mathit{Ca}}^{\mathit{2}\mathit{+}}$ Current

According to previous researches, we accept that CRAC channels and plasma membrane ${\mathrm{Ca}}^{2+}$-ATPase (PMCA) channels are the main ${\mathrm{Ca}}^{2+}$ channels in MCs [33]. CRAC channels are ${\mathrm{Ca}}^{2+}$ influx channels, and PMCA channels are ${\mathrm{Ca}}^{2+}$ outflux channels. The CRAC current is given by the Hodgkin–Huxley (HH) model [35], as follows:

$$I_{\mathrm{CRAC}}=g_{\mathrm{CRAC}}\cdot P_{\mathrm{CRAC}}\cdot (E_{\mathrm{m}}-E_{\mathrm{Ca}})$$

(1)

where $g_{\mathrm{CRAC}}$ is the conductance, $E_{\mathrm{m}}$ is the membrane potential, and $E_{\mathrm{Ca}}=\varphi \cdot \log \frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{e}}}{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}}$ is the Nernst potential for ${\mathrm{Ca}}^{2+}$, $\varphi =\frac{RT}{zF}$, where $R$ is the universal gas constant, $T$ is the absolute temperature, $z$ = 2 is the valence of ${\mathrm{Ca}}^{2+}$, and $F$ is the Faraday constant. $P_{\mathrm{CRAC}}$ is the proportion of CRAC channels in open state, and it is assumed as follows [35]:

$$P_{\mathrm{CRAC}}=\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{act}\frac{1}{2}}}{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{act}\frac{1}{2}}+{{[\mathrm{Ca}}^{2+}]}_{\mathrm{ER}}}$$

(2)

The PMCA current is given by the following [36]:

$$I_{\mathrm{PMCA}}=I_{\mathrm{PMCA},\mathrm{M}}\cdot \frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}}{K_{\mathrm{PMCA}}+{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}}$$

(3)

where $I_{\mathrm{PMCA},\mathrm{M}}$ is the maximum PMCA current, and $K_{\mathrm{PMCA}}$ is the ${\mathrm{Ca}}^{2+}$ concentration for the half activation of PMCA channels.

2.2. ${\mathit{Ca}}^{\mathit{2}\mathit{+}}$ Outflows from ER

The calcium outflow from the ER to Cyt or $\mathsf{\mu }\mathrm{d}$, through ${\mathrm{IP}}_3\mathrm{R}$ channels, is defined as follows:

$$J_{\mathrm{IP}3\mathrm{R}}=(1-C_{\mathrm{IP}3\mathrm{R}})\cdot (V_{\mathrm{IP}3\mathrm{R}}{P}_{\mathrm{oIP}3\mathrm{R}})\cdot {([\mathrm{Ca}}^{2+}{]}_{\mathrm{ER}}-{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}})$$

(4)

$$J_{{\mathrm{IP}3\mathrm{R}}_{\mathsf{\mu }\mathrm{d}}}=C_{\mathrm{IP}3\mathrm{R}}\cdot (V_{\mathrm{IP}3\mathrm{R}}{P}_{{\mathrm{oIP}3\mathrm{R}}_{\mathsf{\mu }\mathrm{d}}})\cdot {([\mathrm{Ca}}^{2+}{]}_{\mathrm{ER}}-{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}})$$

(5)

where $V_{\mathrm{IP}3\mathrm{R}}$ is the maximum total flux through ${\mathrm{IP}}_3\mathrm{R}$ channels [33], and $P_{\mathrm{oIP}3\mathrm{R}}$ and $P_{{\mathrm{oIP}3\mathrm{R}}_{\mathsf{\mu }\mathrm{d}}}$ are the open probabilities of ${\mathrm{IP}}_3\mathrm{R}$ channels facing the Cyt and μd, respectively, and they are defined as follows [33]:

$$P_{\mathrm{oIP}3\mathrm{R}}=S_{\mathrm{act}}^4+4S_{\mathrm{act}}^3\cdot (1-S_{\mathrm{act}})$$

(6)

$$S_{\mathrm{act}}=(\frac{{{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}}{{{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}+d_1})\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}}{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}+d_5})\cdot h$$

(7)

$$P_{{\mathrm{oIP}3\mathrm{R}}_{\mathsf{\mu }\mathrm{d}}}=S_{{\mathrm{act}}_{\mathsf{\mu }\mathrm{d}}}^4+4S_{{\mathrm{act}}_{\mathsf{\mu }\mathrm{d}}}^3\cdot (1-S_{{\mathrm{act}}_{\mathsf{\mu }\mathrm{d}}})$$

(8)

$$S_{\mathrm{a}\mathrm{c}{\mathrm{t}}_{\mathsf{\mu }\mathrm{d}}}=(\frac{{{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}}{{{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}+d_1})\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}}{{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}+d_5})\cdot h_{\mathsf{\mu }\mathrm{d}}$$

(9)

where $S_{\mathrm{act}}$ and $S_{\mathrm{a}\mathrm{c}{\mathrm{t}}_{\mathsf{\mu }\mathrm{d}}}$ express the probability of the activated subunit, respectively, and are defined by sigmoidal functions of ${[\mathrm{IP}}_3]$ and ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$, and $h$ is the slow inactivation gating variable, defined as follows:

$$\frac{dh}{dt}={\alpha }_{\mathrm{h}}(1-h)-{\beta }_{\mathrm{h}}h$$

(10)

$${\alpha }_{\mathrm{h}}=a_2{d}_2\cdot (\frac{{{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}+d_1}{{{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}+d_3})$$

(11)

$${\beta }_{\mathrm{h}}=a_2\cdot {{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$$

(12)

for $\mathsf{\mu }\mathrm{d}$, it is the following:

$$\frac{d{h}_{\mathsf{\mu }\mathrm{d}}}{dt}={\alpha }_{\mathrm{h}}(1-h_{\mathsf{\mu }\mathrm{d}})-{\beta }_{{\mathrm{h}}_{\mathsf{\mu }\mathrm{d}}}{h}_{\mathsf{\mu }\mathrm{d}}$$

(13)

$${\beta }_{{\mathrm{h}}_{\mathsf{\mu }\mathrm{d}}}=a_2\cdot {{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$$

(14)

where $a_2$, $d_1$, $d_2$, $d_3$, and $d_5$ are parameters.

For ${\mathrm{IP}}_3$ dynamics, the production speed of ${\mathrm{IP}}_3$ is related to PLC, and the production of phospholipase C isoforms depends on ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$, so the production speed of ${\mathrm{IP}}_3$ is defined as follows:

$$J_{{\mathrm{IP}}_3{\mathrm{pro}}_{\mathrm{Cyt}}}=V_{\mathrm{PLC}}\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^2}{K_{\mathrm{PLC}}^2+{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^2}$$

(15)

for $\mathsf{\mu }\mathrm{d}$, it is as follows:

$$J_{{\mathrm{IP}}_3{\mathrm{pro}}_{\mathsf{\mu }\mathrm{d}}}=V_{\mathrm{PLC}}\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}^2}{K_{\mathrm{PLC}}^2+{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}^2}$$

(16)

where $V_{\mathrm{PLC}}$ is the maximal production rate of PLC isoforms, and $K_{\mathrm{PLC}}$ is the sensitivity of PLC to ${\mathrm{Ca}}^{2+}$.

${\mathrm{IP}}_3$ is degraded through phosphorylation by ${\mathrm{IP}}_3$ kinases. The kinetic equation can be written as follows:

$$J_{{\mathrm{IP}}_3{\mathrm{deg}}_{\mathrm{Cyt}}}=k_{\mathrm{deg}}\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^2}{K_{\mathrm{deg}}^2{ + [\mathrm{Ca}}^{2+}{]}_{\mathrm{Cyt}}^2}{{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$$

(17)

for $\mathsf{\mu }\mathrm{d}$, it as follows:

$$J_{{\mathrm{IP}}_3{\mathrm{deg}}_{\mathsf{\mu }\mathrm{d}}}=k_{\mathrm{deg}}\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}^2}{K_{\mathrm{deg}}^2{ + [\mathrm{Ca}}^{2+}{]}_{\mathsf{\mu }\mathrm{d}}^2}{{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}$$

(18)

where $k_{\mathrm{deg}}$ represents the phosphorylation rate constants, and $K_{\mathrm{deg}}$ is the half-saturation constant of ${\mathrm{IP}}_3$ kinases.

${\mathrm{IP}}_3$ leaks from μd to Cyt can be defined as follows:

$$J_{{\mathrm{IP}}_3\mathrm{leak}}=k_{{\mathrm{IP}}_3\mathrm{leak}}{([\mathrm{IP}}_3{]}_{\mathsf{\mu }\mathrm{d}}-{{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}).$$

(19)

Therefore, the change in ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ and ${{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}$ can be written as follows:

$$\frac{d{{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}}{dt}=J_{{\mathrm{IP}}_3{\mathrm{pro}}_{\mathrm{Cyt}}}-J_{{\mathrm{IP}}_3{\mathrm{deg}}_{\mathrm{Cyt}}}+J_{{\mathrm{IP}}_3\mathrm{leak}}$$

(20)

$$\frac{d{{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}}{dt}=J_{{\mathrm{IP}}_3{\mathrm{pro}}_{\mathsf{\mu }\mathrm{d}}}-J_{{\mathrm{IP}}_3{\mathrm{deg}}_{\mathsf{\mu }\mathrm{d}}}-J_{{\mathrm{IP}}_3\mathrm{leak}}$$

(21)

SERCA pumps transport ${\mathrm{Ca}}^{2+}$ into the ER, and the flux from Cyt to the ER is defined as follows:

$$J_{\mathrm{SERCA}}=(1-C_{\mathrm{SERCA}})\cdot V_{\mathrm{SERCA}}\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^2}{k_{\mathrm{SERCA}}^2+{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^2})$$

(22)

the flux from μd to the ER is defined as follows:

$$J_{{\mathrm{SERCA}}_{\mathsf{\mu }\mathrm{d}}}=C_{\mathrm{SERCA}}\cdot V_{\mathrm{SERCA}}\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}^2}{k_{\mathrm{SERCA}}^2+{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}^2})$$

(23)

where $V_{\mathrm{SERCA}}$ is the maximal flux through SERCA, and $k_{\mathrm{SERCA}}$ is the ${\mathrm{Ca}}^{2+}$ activation constant for SERCA.

${\mathrm{Ca}}^{2+}$ leaks from the the ER are driven by the concentration gradients between the ER and either the Cyt or μd. The leak from the ER into the Cyt is defined as follows:

$$J_{{\mathrm{leak}}_{\mathrm{Cyt}}^{\mathrm{ER}}}=k_{\mathrm{Cyt}}^{\mathrm{ER}}\cdot {([\mathrm{Ca}}^{2+}{]}_{\mathrm{ER}}-{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}})$$

(24)

and the leak from the ER into the μd is defined as follows:

$$J_{{\mathrm{leak}}_{\mathsf{\mu }\mathrm{d}}^{\mathrm{ER}}}=k_{\mathsf{\mu }\mathrm{d}}^{\mathrm{ER}}\cdot {([\mathrm{Ca}}^{2+}{]}_{\mathrm{ER}}-{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}})$$

(25)

Although the μd is not a membrane-bound compartment, we similarly define the leak from the μd into the Cyt as follows:

$$J_{{\mathrm{leak}}_{\mathrm{Cyt}}^{\mathsf{\mu }\mathrm{d}}}=k_{\mathrm{Cyt}}^{\mathsf{\mu }\mathrm{d}}\cdot {([\mathrm{Ca}}^{2+}{]}_{\mathsf{\mu }\mathrm{d}}-{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}})$$

(26)

2.3. ${\mathit{Ca}}^{\mathit{2}\mathit{+}}$ Outflows from Mt

In the Mt, the mNCX channels exchange one ${\mathrm{Ca}}^{2+}$ for three ${\mathrm{Na}}^{+}$. Flux through mNCX channels facing the Cyt is defined as follows:

$$J_{\mathrm{mNCX}}=(1-C_{\mathrm{mNCX}})\cdot V_{\mathrm{mNCX}}\cdot (\frac{{{[\mathrm{Na}}^{+}]}_{\mathrm{Cyt}}^3}{k_{\mathrm{Na}}^3+{{[\mathrm{Na}}^{+}]}_{\mathrm{Cyt}}^3})\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}}{k_{\mathrm{mNCX}}+{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}})$$

(27)

and for mNCX channels facing the μd, it is as follows:

$$J_{{\mathrm{mNCX}}_{\mathsf{\mu }\mathrm{d}}}=C_{\mathrm{mNCX}}\cdot V_{\mathrm{mNCX}}\cdot (\frac{{{[\mathrm{Na}}^{+}]}_{\mathsf{\mu }\mathrm{d}}^3}{k_{\mathrm{Na}}^3+{{[\mathrm{Na}}^{+}]}_{\mathsf{\mu }\mathrm{d}}^3})\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}}{k_{\mathrm{mNCX}}+{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}})$$

(28)

where ${{[\mathrm{Na}}^{+}]}_{\mathrm{Cyt}}$ and ${{[\mathrm{Na}}^{+}]}_{\mathsf{\mu }\mathrm{d}}$ are the concentration of ${\mathrm{Na}}^{+}$ in Cyt and μd, respectively. $V_{\mathrm{mNCX}}$ is the maximal flux through the mNCX, and $k_{\mathrm{Na}}$ and $k_{\mathrm{mNCX}}$ are ${\mathrm{Na}}^{+}$ and ${\mathrm{Ca}}^{2+}$ activation constants for mNCX, respectively. The connectivity coefficient $C_{\mathrm{mNCX}}$ is the proportion of mNCX channels facing the μd.

The MCU channel transports ${\mathrm{Ca}}^{2+}$ into the Mt. Flux through the MCU to the Cyt is defined as follows:

$$J_{\mathrm{MCU}}=(1-C_{\mathrm{MCU}})\cdot V_{\mathrm{MCU}}\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^2}{k_{\mathrm{MCU}}^2+{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^2})$$

(29)

and to the μd, it is as follows:

$$J_{{\mathrm{MCU}}_{\mathsf{\mu }\mathrm{d}}}=C_{\mathrm{MCU}}\cdot V_{\mathrm{MCU}}\cdot (\frac{{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}^2}{k_{\mathrm{MCU}}^2+{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}^2})$$

(30)

where $V_{\mathrm{MCU}}$ $=$ $V_{{\mathrm{MCU}}_0}\mathsf{\Delta }\mathsf{\Phi }$, and $\mathsf{\Delta }\mathsf{\Phi }=\frac{bF(\mathsf{\Psi }{ - \mathsf{\Psi }}_0)}{RT}{\mathrm{e}}^{\frac{bF{(\mathsf{\Psi } - \mathsf{\Psi }}_0)}{RT}}\sinh \frac{bF{(\mathsf{\Psi } - \mathsf{\Psi }}_0)}{RT}$. $V_{{\mathrm{MCU}}_0}$ represents the maximal flux through the MCU, and $\mathsf{\Delta }\mathsf{\Phi }$ is the voltage driving force. $\mathsf{\Psi }$ is the inner mitochondrial membrane voltage (150~180 mV, negative inside). $b$ and ${\mathsf{\Psi }}_0$ are the fitting parameters obtained from Qi [33]. During the simulations, we assume a constant $\mathsf{\Psi }$ of 170 mV, as experimental evidence suggests that $\mathsf{\Psi }$ does not change significantly in response to transient cytosolic ${[\mathrm{Ca}}^{2+}]$ increase, produced by ${\mathrm{IP}}_3$-generating agonists [37,38,39]. $k_{\mathrm{MCU}}$ is the ${\mathrm{Ca}}^{2+}$ activation constant for MCU, and the connectivity coefficient $C_{\mathrm{MCU}}$ is the proportion of MCU channels facing the μd.

2.4. Effective Cytosol

Since the μd is a part of the cytosol, we defined the ${[\mathrm{Ca}}^{2+}]$ of the effective cytosolic compartment ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ as the volume-weighted average of ${[\mathrm{Ca}}^{2+}]$ within the combined Cyt and μd compartments.

$${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}=\frac{Vo{l}_{\mathrm{Cyt}}\cdot {{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}+Vo{l}_{\mathsf{\mu }\mathrm{d}}\cdot {{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}}{Vo{l}_{\mathrm{Cyt}}+Vo{l}_{\mathsf{\mu }\mathrm{d}}}$$

(31)

2.5. μd Volume

We assumed that each mitochondrion is a sphere, and 20% of its surface area closes to the ER [18,40]. Some experimental data suggest that the diameters of mitochondria are 0.5~1.5 μm [41]; here, we choose 0.58 μm. There are about two hundred (N) mitochondria in each cell. Thus, we calculate the volume of the μd compartment as follows:

$$Vo{l}_{\mathsf{\mu }\mathrm{d}}=0.2\cdot SA\cdot N\cdot D$$

(32)

where $SA$ is the surface area, $N$ is the mitochondria total number, and $D$ is the ER–Mt distance.

2.6. Temporal Changes in ${[\mathrm{Ca}}^{2+}]$ in Each Compartment

Temporal changes in ${[\mathrm{Ca}}^{2+}]$ in each compartment are represented as the following ordinary differential equations:

in cytosol, it is the following:

$$\frac{d{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}}{dt}=\frac{(J_{\mathrm{IP}3\mathrm{R}}+J_{\mathrm{mNCX}}+J_{{\mathrm{leak}}_{\mathrm{Cyt}}^{\mathsf{\mu }\mathrm{d}}}+J_{{\mathrm{leak}}_{\mathrm{Cyt}}^{\mathrm{ER}}}-J_{\mathrm{SERCA}}-J_{\mathrm{MCU}})}{1+{\theta }_{\mathrm{Cyt}}}+\frac{S_c}{zFVo{l}_{\mathrm{Cyt}}}(I_{\mathrm{CRAC}}-I_{\mathrm{PMCA}})$$

(33)

in the ER, it is the following:

$$\frac{d{{[\mathrm{Ca}}^{2+}]}_{\mathrm{ER}}}{dt}=\frac{Vo{l}_{\mathrm{Cyt}}}{Vo{l}_{\mathrm{ER}}}(J_{\mathrm{SERCA}}+J_{{\mathrm{SERCA}}_{\mathsf{\mu }\mathrm{d}}}-J_{\mathrm{IP}3\mathrm{R}}-J_{{\mathrm{IP}3\mathrm{R}}_{\mathsf{\mu }\mathrm{d}}}-J_{{\mathrm{leak}}_{\mathsf{\mu }\mathrm{d}}^{\mathrm{ER}}}-J_{{\mathrm{leak}}_{\mathrm{Cyt}}^{\mathrm{ER}}})/(1+{\theta }_{\mathrm{ER}})$$

(34)

in the Mt, it is the following:

$$\frac{d{{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}}{dt}=\frac{Vo{l}_{\mathrm{Cyt}}}{Vo{l}_{\mathrm{Mt}}}(J_{\mathrm{MCU}}+J_{{\mathrm{MCU}}_{\mathsf{\mu }\mathrm{d}}}-J_{\mathrm{mNCX}}-J_{{\mathrm{mNCX}}_{\mathsf{\mu }\mathrm{d}}})/(1+{\theta }_{\mathrm{Mt}})$$

(35)

in μd, it is the following:

$$\frac{d{{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}}{dt}=\frac{Vo{l}_{\mathrm{Cyt}}}{Vo{l}_{\mathsf{\mu }\mathrm{d}}}(J_{{\mathrm{IP}3\mathrm{R}}_{\mathsf{\mu }\mathrm{d}}}+J_{{\mathrm{mNCX}}_{\mathsf{\mu }\mathrm{d}}}+J_{{\mathrm{leak}}_{\mathsf{\mu }\mathrm{d}}^{\mathrm{ER}}}-J_{{\mathrm{SERCA}}_{\mathsf{\mu }\mathrm{d}}}-J_{{\mathrm{MCU}}_{\mathsf{\mu }\mathrm{d}}}-J_{{\mathrm{leak}}_{\mathrm{Cyt}}^{\mathsf{\mu }\mathrm{d}}})/(1+{\theta }_{\mathsf{\mu }\mathrm{d}})$$

(36)

where ${\theta }_{\mathrm{i}}$ (i = Cyt, ER, Mt, μd) is the buffer factor of each compartment, which is defined as follows [36]:

$${\theta }_{\mathrm{i}}=\frac{B{P}_{\mathrm{i}}{K}_{\mathrm{i}}}{{({{[\mathrm{Ca}}^{2+}]}_{\mathrm{i}}+K_{\mathrm{i}})}^2}$$

(37)

The parameter values related to our model are given in

Table 1. Parameters of the model.

Parameter	Value	Description
$g_{\mathrm{CRAC}}$	0.3 ${\mathsf{\Omega }}^{-1}{\mathrm{m}}^{-2}$	the conductance [19]
$E_{\mathrm{m}}$	$-60\mathrm{mV}$	the membrane potential [19]
$R$	$8{.34\mathrm{J}\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}$	the universal gas constant [19]
$T$	$293\mathrm{K}$	the absolute temperature [19]
$F$	${\mathrm{96,485}\mathrm{C}\mathrm{mol}}^{-1}$	Faraday constant [19]
${{[\mathrm{Ca}}^{2+}]}_{\mathrm{act}1/2}$	$5\times {10}^{-4}{\mathrm{mol}\mathrm{L}}^{-1}$	the Ca2+ concentration for half activation of SOC [19]
$I_{\mathrm{PMCA},\mathrm{M}}$	${89.9\mathsf{\mu }\mathrm{M}\mathrm{s}}^{-1}$	the maximum PMCA current [36]
$K_{\mathrm{PMCA}}$	$0.26\mathsf{\mu }\mathrm{M}$	the Ca2+ concentration for half activation of PMCA channels [36]
${{[\mathrm{Ca}}^{2+}]}_{\mathrm{e}}$	$2000\mathsf{\mu }\mathrm{M}$	the extracellular Ca2+ concentration [19]
$Vo{l}_{\mathrm{ER}}$	$0.1\mathrm{pL}$	volume of ER [41]
$Vo{l}_{\mathrm{Mt}}$	$0.05\mathrm{pL}$	volume of Mt [41]
$Vo{l}_{\mathrm{Cyt}}$	$0.85\mathrm{pL}$	volume of Cyt [34]
$S_c$	$0{.28\mathrm{pm}}^2$	cell surface [34]
$V_{\mathrm{IP}3\mathrm{R}}$	$1{.59\mathrm{s}}^{-1}$	max flux of IP3R [34]
$V_{\mathrm{SERCA}}$	$29{.1\mathsf{\mu }\mathrm{M}\mathrm{s}}^{-1}$	max flux of SERCA pump [34]
$k_{\mathrm{SERCA}}$	$0.193\mathsf{\mu }\mathrm{M}$	activation constant for SERCA pump [34]
$a_2$	$0{.0605\mathsf{\mu }\mathrm{M}}^{-1}{\mathrm{s}}^{-1}$	IP3R binding rate at Ca2+ inhibition sites [34]
$d_1$	$0.0377\mathsf{\mu }\mathrm{M}$	IP3R dissociation constant for IP3 sites [34]
$d_2$	$1.33\mathsf{\mu }\mathrm{M}$	IP3R dissociation constant for Ca2+ inhibition sites [34]
$d_3$	$1.74\mathsf{\mu }\mathrm{M}$	IP3R dissociation constant for IP3 sites [34]
$d_5$	$0.239\mathsf{\mu }\mathrm{M}$	IP3R dissociation constant for Ca2+ activation sites [34]
$V_{\mathrm{MCU}}$	$7{.53\mathsf{\mu }\mathrm{M}\mathrm{s}}^{-1}$	max rate of Ca2+ uptake by MCU [34]
$k_{\mathrm{MCU}}$	$1.23\mathsf{\mu }\mathrm{M}$	half-max rate of Ca2+ pumping from Cyt to Mt [34]
$V_{\mathrm{NCX}}$	${119\mathsf{\mu }\mathrm{M}\mathrm{s}}^{-1}$	max rate of Ca2+ release through NCX [34]
$k_{\mathrm{NCX}}$	$43.2\mathsf{\mu }\mathrm{M}$	activation constant for NCX [34]
$k_{\mathrm{Na}}$	$9.4\mathrm{mM}$	Na+ activation constant for MCU [34]
${[\mathrm{Na}]}_{\mathrm{Cyt}}$	$10\mathrm{mM}$	Na+ in Cyt [34]
${[\mathrm{Na}]}_{\mathsf{\mu }\mathrm{d}}$	$10\mathrm{mM}$	Na+ in μd [34]
$k_{\mathsf{\mu }\mathrm{d}}^{\mathrm{ER}}$	$0{.0433\mathrm{s}}^{-1}$	leak constant from ER to μd [34]
$k_{\mathrm{Cyt}}^{\mathrm{ER}}$	$0{.0107\mathrm{s}}^{-1}$	leak constant from ER to Cyt [34]
$k_{\mathrm{Cyt}}^{\mathsf{\mu }\mathrm{d}}$	$0{.0332\mathrm{s}}^{-1}$	leak constant from μd to Cyt [34]
$C_{\mathrm{IP}3\mathrm{R}}$	0.486	fraction of IP3R facing microdomain [34]
$C_{\mathrm{SERCA}}$	0.603	fraction of SERCA facing microdomain [34]
$C_{\mathrm{MCU}}$	0.894	fraction of MCU facing microdomain [34]
$C_{\mathrm{mNCX}}$	0.569	fraction of mNCX facing microdomain [34]
$B{P}_{\mathrm{Cyt}}$	$154\mathsf{\mu }\mathrm{M}$	total buffer concentration in Cyt [42]
$K_{\mathrm{Cyt}}$	11.1	buffer rate constant ratio [42]
$B{P}_{\mathrm{ER}}$	$11,100\mathsf{\mu }\mathrm{M}$	total buffer concentration in ER [42]
$K_{\mathrm{ER}}$	967	buffer rate constant ratio [42]
$B{P}_{\mathrm{Mt}}$	$285,000\mathsf{\mu }\mathrm{M}$	total buffer concentration in Mt [42]
$K_{\mathrm{Mt}}$	698	buffer rate constant ratio [42]
$B{P}_{\mathsf{\mu }\mathrm{d}}$	$191\mathsf{\mu }\mathrm{M}$	total buffer concentration in μd [42]
$K_{\mathsf{\mu }\mathrm{d}}$	12	buffer rate constant ratio [42]

.

3. Results

3.1. Effect of the Degradation and Production of ${\mathit{IP}}_{\mathit{3}}$ on ${\mathit{Ca}}^{\mathit{2}\mathit{+}}$ Oscillations

Based on the ODE45 solver of MATLAB, we regard all formulas as a system to program and solve. Therefore, the results shown below involve all the formulas. Except for the parameter values listed in the caption of each figure, all the other parameter values and meanings are in Table 1.

What we study is biological signals. Due to the compensatory effect of biological systems, there are a few signals with stable periodic oscillations. After a period of time, the oscillation signals often return to the original equilibrium value, or reach a new equilibrium value. Numerical oscillation analysis diagrams imitate the bifurcation diagram of the dynamic system, to analyze the conditions that can produce periodic oscillation solutions for a certain length of time (1000 s in our calculation). We make numerical oscillation analysis diagrams by finding the maximum and minimum values of ${\mathrm{Ca}}^{2+}$ concentration at different $k_{\mathrm{deg}}$, $V_{\mathrm{PLC}}$, and ${\mathrm{IP}}_3$ values in the last 100 s (900–1000 s). If the maximum value and minimum value are not the same, it indicates that the ${\mathrm{Ca}}^{2+}$ concentration is in the state of oscillation, and numerical oscillation analysis diagrams show that one $k_{\mathrm{deg}}$ or $V_{\mathrm{PLC}}$ value corresponds to two ${\mathrm{Ca}}^{2+}$ concentrations. If the maximum value and minimum value are the same, it indicates that the ${\mathrm{Ca}}^{2+}$ concentration is in equilibrium state, and numerical oscillation analysis diagrams show that one $k_{\mathrm{deg}}$ or $V_{\mathrm{PLC}}$ value corresponds to one ${\mathrm{Ca}}^{2+}$ concentration.

The numerical oscillation analysis diagrams show the effect of Mt and μd compartments on the ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillatory dynamics, as functions of $k_{\mathrm{deg}}$ and $V_{\mathrm{PLC}}$, respectively (Figure 2). They also show that the Mt compartment results in a slight decrease in the predicted oscillatory amplitude, and the μd compartment leads to an obvious decrease and the emergence of an oscillatory region at high levels of $k_{\mathrm{deg}}$ (range 4 to 9 s⁻¹) (Figure 2a). These results are consistent with those found by Arash Moshkforoush [34], which showed that the addition of the μd compartment resulted in the appearance of the low-energy ${\mathrm{IP}}_3$ oscillations region in numerical oscillation analysis diagrams. $k_{\mathrm{deg}}$ is positively correlated with ${\mathrm{IP}}_3$ degradation rate, and high levels of $k_{\mathrm{deg}}$ mean that ${\mathrm{IP}}_3$ degrades fast; therefore, ${\mathrm{IP}}_3$ should be kept at low levels. The Mt compartment and μd compartment can cause a decrease in the oscillatory amplitude (Figure 2b). Without μd, ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations appear when $V_{\mathrm{PLC}}$ is in the range of 0~0.004 μMs⁻¹. While with μd, ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations appear when $V_{\mathrm{PLC}}$ is in the range of 0 to 0.003 μMs⁻¹. Figure 2a or Figure 2b shows there are little differences between the three curves at the equilibrium of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$. This indicated that the Mt and μd have little effect on the intracellular equilibrium of ${\mathrm{Ca}}^{2+}$ concentration.

Consider the numerical oscillation analysis diagrams in Figure 2a, which show that with Mem, Mt, and μd, the oscillation range of $k_{\mathrm{deg}}$ is 0.64 to 2.85 ${\mathrm{s}}^{-1}$ and 4.01 to 8.93 ${\mathrm{s}}^{-1}$. Therefore, four values of $k_{\mathrm{deg}}$ (0.1, 0.5, 1.5, and 5.0 ${\mathrm{s}}^{-1}$) are chosen to simulate the temporal traces of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$. According to Figure 3a,b, ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ rises first, due to the stimulation, then fluctuates for a period of time and returns to the equilibrium resting state, when $k_{\mathrm{deg}}$ is 0.1 ${\mathrm{s}}^{-1}$ and 0.5 ${\mathrm{s}}^{-1}$,which are both out of the oscillation range. While, when $k_{\mathrm{deg}}$ is 1.5 ${\mathrm{s}}^{-1}$ or 5.0 ${\mathrm{s}}^{-1}$, which both are in the oscillation range, ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ will maintain steady oscillations after stimulation, as shown in Figure 3c,d. The amplitude of oscillations is higher (~2 μM) and the frequency is lower (~2 oscillations/min) at lower levels of $k_{\mathrm{deg}}$ (Figure 3c, 1.5 ${\mathrm{s}}^{-1}$), compared to those at higher levels of $k_{\mathrm{deg}}$ (amplitudes: ~0.2 μM, frequencies: ~4 oscillations/min, Figure 3d, 5.0 ${\mathrm{s}}^{-1}$).

When $V_{\mathrm{PLC}}$ is 0.001 and 0.003 ${\mathsf{\mu }\mathrm{M}\mathrm{s}}^{-1}$, ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ does not rise immediately, but fluctuates to a value and forms oscillations subsequently. However, the calcium oscillations in Figure 4b are not stable, they return to an equilibrium state over time. In Figure 4c,d, ${\mathrm{Ca}}^{2+}$ concentrations only oscillate for a short period of time, to return to equilibrium. Comparing the four graphs (Figure 4a–d) shows that, with the increase in $V_{\mathrm{PLC}}$, the time required for forming oscillations becomes shorter and shorter until it disappears, the frequency of oscillations decreases, but the amplitude and final equilibrium value increase. When $V_{\mathrm{PLC}}$ increases to a certain value, the oscillations disappear. This means that $V_{\mathrm{PLC}}$ affects the process of calcium oscillations from one equilibrium state to another equilibrium state. These results are consistent with most of the biological signals. Similar to the neural electrical signal, after a period of time, the oscillation signals often return to the original equilibrium value or reach a new equilibrium value.

The temporal profiles of ${[\mathrm{Ca}}^{2+}]$ in the following cellular compartments: Cyt, ER, Mt, μd, after stimulation, are shown in Figure 5a. When ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{ER}}$ is at the valley value (163 μM), ${[\mathrm{Ca}}^{2+}]$ in other cellular compartments are at peak value, as shown in Figure 5a. This indicates that the main ${\mathrm{Ca}}^{2+}$ filling of the Mt and μd come from the ER. The numerical simulation results also prove that high ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}$ is observed when global ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ is lower, and ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ is 20 times that of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$, when ${\mathrm{Ca}}^{2+}$ outflows from the ER. The peak value of ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ (37.3 μM) appears slightly earlier than those of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ (1.6 μM) and ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Mt}}$ (3.5 μM), shown in Figure 5a. This indicates that ${\mathrm{Ca}}^{2+}$ oscillations in μd are not completely synchronized with those in Cyt and Mt. The temporal profiles of ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ and ${{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}$ are shown in Figure 5b. After stimulation, ${{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}$ decrease appears about 10 s earlier than that of ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$, which is in accordance with the fact that the ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ increase happens earlier than the ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ increase, shown in Figure 5a. The valley value of ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ (0.56 μM) is lower than that of ${{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}$ (0.61 μM) in Figure 5b.

3.2. Effect of the ER–Mt Distance (D) on ${\mathit{Ca}}^{\mathit{2}\mathit{+}}$ Oscillations

Figure 6a,c show that with the increase in the ER–Mt distance, the amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations increased slightly. The numerical simulation results of Qi [33] show that with increasing ER–Mt distance at D < 20 nm, the ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ amplitudes decrease, while at D > 20 nm, the ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ amplitudes increase. A result similar to Qi appears in our model, at D = 40 nm, the ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ amplitude is the highest in Figure 6b. Moreover, Figure 6d shows that with ER–Mt distance increases at D < 65 nm, the ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ amplitudes increase, while at D > 65 nm, the ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ amplitudes decrease with ER–Mt distance increases. In Figure 6f, there is also an obvious inflection point of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ amplitudes at D = 65 nm. When D > 65 nm, the ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ amplitudes increase faster than D < 65 nm, of which the amplitudes maintain around 1.6 μM. Combining Figure 6d,f, we find that in this model, when D is larger than 65 nm, the influence of μd on ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations is weak. This is a mathematical explanation of why the μd functional region should have a small distance. Figure 6e shows that ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{ER}}$ amplitudes increase with ER–Mt distance increases, meaning the ER releases more ${\mathrm{Ca}}^{2+}$ under larger ER–Mt distances. This explains why the amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ and ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations increase with ER–Mt distance increases.

Figure 7 shows that the period of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ and ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ oscillations increases with ER–Mt distance increases. In combination with Figure 5a, we find that ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ is higher than ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ and their oscillations are not synchronous, but the periods as well as the frequency of oscillations are the same. When D = 10, 200 nm, the period of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ and ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ is 23.2 and 23.3 s, and 25.5 and 25.5 s, respectively. Figure 7 shows that the period and distance have an approximate linear relationship; therefore, we can obtain the slopes of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ (0.01244) and ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ (0.01232) by fitting. From a mathematical point, this result is beyond our expectation, because in this model, Cyt and μd are calculated as two rooms, and there is only a diffusion relationship between the two rooms. While from a cellular physiological point, this result is reasonable. In the model, the μd is a region that we hypothesize from the cytoplasm. However, in the actual cell, it is a part of the Cyt, hence both the frequencies should be the same. Meanwhile, the frequency of calcium oscillations is one of the ways that cellular calcium signaling transmits information. When the whole intracellular calcium signaling is formed, the same frequency can transmit the same information. Therefore, the calcium signaling must be consistent to prevent cells from receiving different information at the same time, causing functional disorders.

3.3. Effect of the ${\mathit{[}{\mathit{IP}}_{\mathit{3}}\mathit{]}}_{\mathit{Cyt}}$ on ${\mathit{Ca}}^{\mathit{2}\mathit{+}}$ Oscillations

The numerical oscillation analysis diagrams in Figure 8a show that the Mt, μd, and Mem compartments can each reduce the amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations. Contrasting the dot curve with the solid curve, we can find that the addition of the μd compartment causes ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillation regions at low-level ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ of 0.08~0.27 μM. Comparing the dot curve and dot solid curve, we find that the Mem compartment makes the ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations at low-level ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ disappear. The Mt compartment has a slight effect on the left and right bifurcation point value of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations. The ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations region of the model, with Mt, μd, and Mem, shrink obviously, meaning that the μd and Mem compartments limit the range of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations. From Figure 8a, we can also draw a similar conclusion as Figure 2, which is that the presence or absence of Mt, μd, and Mem have little effect on the equilibrium calcium concentration.

4. Discussion and Conclusions

In Qi’s model [33], the ${\mathrm{IP}}_3\mathrm{R}$–MCU distance is regarded as the main factor by which μd affects calcium oscillations. He links the ER with Mt by one-dimensional diffusion of ${\mathrm{Ca}}^{2+}$ between ${\mathrm{IP}}_3\mathrm{R}$ and MCU. However, with the distance increases, the influence of ${\mathrm{Ca}}^{2+}$ diffusion in other directions will be more significant. Hence, the one-dimensional diffusion assumption will make the calculation error of ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ larger, and the error of calcium oscillations larger. Therefore, in our model, we assume that the μd is a separate chamber with volume. As shown in Figure 6d, ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ increases a little before 65 nm and decrease subsequently. From Figure 6e, we find that the ER–Mt distance increase makes the ER release more ${\mathrm{Ca}}^{2+}$. This is because ${\mathrm{IP}}_3\mathrm{R}$ has a ${\mathrm{Ca}}^{2+}$ inhibition binding site, so when ${\mathrm{Ca}}^{2+}$ binds to this site, ${\mathrm{IP}}_3$ activity is inhibited. As the ER–Mt distance increases, ${\mathrm{Ca}}^{2+}$ can spread faster; therefore, the probability of ${\mathrm{Ca}}^{2+}$ binding to the inhibition site decreases, so ${\mathrm{Ca}}^{2+}$ released by ${\mathrm{IP}}_3$ increases. Hence, there is a little upward trend of ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$, and, later, the effect of the increased volume of μd is greater than the release of ${\mathrm{Ca}}^{2+}$, and ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ begins to decrease. However, in Qi’s result, ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ was decreasing all the time. Arash Moshkforoush’s model [34] does not consider the degradation and production of ${\mathrm{IP}}_3$; therefore, the concentration of ${\mathrm{IP}}_3$ in the cells is constant, which is not a physiological reality. ${\mathrm{IP}}_3$ dynamic behaviors have a significant effect on the range of parameter values and the oscillation patterns of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations [43]. Inhibition of protein kinase C eliminates ${\mathrm{Ca}}^{2+}$ oscillations, while ${\mathrm{IP}}_3$ formation is still maintained [44]. This is also shown in Figure 2a,b, which shows that at different degradation and production levels of ${\mathrm{IP}}_3$, there are different amplitudes and frequencies of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations.

In this paper, the dynamic model of calcium oscillations in MCs is developed, which considers the major cellular compartments (Cyt, Mem, ER, and Mt), ${\mathrm{Ca}}^{2+}$ channels and buffer in these compartments, and the μd composed of the ER and Mt. In our simulations, the Mt and μd compartments can reduce the amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations. With the addition of the μd compartment, an oscillatory region will appear at high levels of $k_{\mathrm{deg}}$ (4 to 9 s⁻¹) and at low levels of ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ (0.08~0.27 μM), shown in Figure 2a and Figure 8a. Our model also shows that different concentrations of ${\mathrm{IP}}_3$ stimulation will change the amplitude and frequency of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations. The amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations increases, and the frequency of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations decreases with ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ increases. Figure 4b shows that the ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}^{\mathrm{eff}}$ oscillations process is in line with the actual law of cell calcium signal generation. The calcium signal can be divided into the following three levels: (i) at first, being the most fundamental event, a very low level of stimulation will cause a brief opening of a single channel and the release of calcium, which is called calcium blips; (ii) then, there is the basic event, which results from a small group of channels opening and the release of calcium, to form calcium sparks; (iii) finally, the synchronization of a large number of fundamental events produces the global calcium signal, and subsequently restores the resting state. According to our results, the amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations increases with the increase in the ER–Mt distance. Moreover, the ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ amplitude also increases with the increase in the ER–Mt distance at D < 65 nm, but decreases with the increase in the ER–Mt distance at D > 65 nm. Therefore, we believe that μd has a better regulation effect on ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations when the ER–Mt distance is less than about 65 nm, which also provides reference for determining the distance of μd in the subsequent studies. The periods of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ and ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ oscillations are the same at different ER–Mt distances, and they increase with the ER–Mt distance. Meanwhile, from Figure 5a and Figure 7, we can understand that ${{[\mathrm{Ca}}^{2+}]}_{\mathsf{\mu }\mathrm{d}}$ is 20 times higher than ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$, and their oscillations are not synchronous. The presence of μd causes the ER to release less ${\mathrm{Ca}}^{2+}$, and the effect of μd decreases with ER–Mt distance increases. This proves that μd acts as a buffer against the release of ${\mathrm{Ca}}^{2+}$ from the ER. All these results suggest that μd plays an important role in controlling ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations at D < 65 nm. The degradation and production of ${\mathrm{IP}}_3$ can also regulate ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations by maintaining different levels of ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$. As shown in Figure 5b, before stimulation, ${{[\mathrm{IP}}_3]}_{\mathrm{Cyt}}$ and ${{[\mathrm{IP}}_3]}_{\mathsf{\mu }\mathrm{d}}$ get closer, due to the diffusion between Cyt and μd. This is because the ${\mathrm{IP}}_3$ production and degradation of μd and Cyt maintain dynamic equilibrium at the initial moment, but the concentration of ${\mathrm{IP}}_3$ will affect ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations; therefore, we assume that the production and degradation of ${\mathrm{IP}}_3$ are both zero before the stimulation (50 s), to reduce the impact of ${\mathrm{IP}}_3$ on calcium oscillations in the study of $V_{\mathrm{PLC}}$ and $k_{\mathrm{deg}}$. Then, at this time, the dynamic equilibrium is destroyed, and the concentrations of μd and Cyt are close to each other. We already know that the concentration of ${\mathrm{IP}}_3$ also has an effect on oscillations, hence the given concentration of ${\mathrm{IP}}_3$ is not used as the stimulation, but the system is deviated from the equilibrium state through diffusion after the production and degradation of ${\mathrm{IP}}_3$ is assumed to be zero, and then the values of $V_{\mathrm{PLC}}$ and $k_{\mathrm{deg}}$ are restored, so as to reduce the impact of ${\mathrm{IP}}_3$ on calcium oscillations in the study of $V_{\mathrm{PLC}}$ and $k_{\mathrm{deg}}$.

From the biological point of view, due to the compensatory effect of biological systems, there are a few signals with stable periodic oscillations. After a period of time, the oscillation signals often return to the original equilibrium value, or reach a new equilibrium value. This is consistent with Figure 4b–d. ${\mathrm{IP}}_3\mathrm{R}$ is regulated by ${\mathrm{Ca}}^{2+}$ in a biphasic manner; therefore, ${\mathrm{IP}}_3\mathrm{R}$ activity is inhibited at higher ${\mathrm{Ca}}^{2+}$ concentrations. When the ER–Mt distance increases in our model, the volume of μd increases and the ${\mathrm{Ca}}^{2+}$ concentration of μd gradually decreases, and this causes the activity of ${\mathrm{IP}}_3\mathrm{R}$ to increase; therefore, more ${\mathrm{Ca}}^{2+}$ will outflow from the ER through the ${\mathrm{IP}}_3\mathrm{R}$ channels, causing the ER calcium oscillations amplitude increases shown in Figure 6e, and the cytoplasmic ${\mathrm{Ca}}^{2+}$ concentration increases shown in Figure 6c. Calcium oscillations have been widely accepted as a universal signal mode in cells. With the in-depth study of calcium oscillations, the basic theory of calcium oscillations regulating downstream biological effects through its frequency has been established [45,46]. In Figure 7, the calcium oscillation frequencies of Cyt and μd calculated by our model are basically the same, which also implies that the frequency of calcium oscillations is one of the ways of signal transmission.

In summary, the study provides a dynamic model that simulates calcium oscillations in mast cells and provides a theoretical basis for the mast cell calcium signal observed in the experiment. This enabled us to consolidate previous theoretical and experimental findings. The model results showed that Mem, Mt, and μd can all reduce the amplitude of ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations. Moreover, μd can play a critical role in ${\mathrm{Ca}}^{2+}$ dynamics at appropriate ER–Mt distances (less than 65 nm). In future work, we will continue to study the influence of mitochondrial ATP and important calcium channel parameters on ${{[\mathrm{Ca}}^{2+}]}_{\mathrm{Cyt}}$ oscillations. Additionally, we will improve our model by adding other intracellular calcium pools (such as nucleus, Golgi apparatus, etc.), to make our model more accurate. Although we have oversimplified some details, we believe that this model is still useful and can provide us with some insights into the mechanism of mast cell calcium signaling regulation.