A Dynamic Model of Cytosolic Calcium Concentration Oscillations in Mast Cells

: In this paper, a dynamic model of cytosolic calcium concentration ( [ Ca 2 + ] 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 Ca 2 + channels in these cellular compartments. By this model, we calculate [ Ca 2 + ] Cyt oscillations that are driven by distinct mechanisms at varying k deg (degradation coefﬁcient of inositol 1,4,5-trisphosphate, IP 3 and production coefﬁcient of 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 [ Ca 2 + ] Cyt oscillations, and cause the ER to release less Ca 2 + during oscillations; (ii) with increasing cytosolic IP 3 concentration ( [ IP 3 ] Cyt ), the amplitude of oscillations increases (from 0.1 µ M to several µ M), but the frequency decreases; (iii) the frequency of [ Ca 2 + ] 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 [ Ca 2 + ] Cyt oscillations. These results suggest that Mt, µ d, and IP 3 can all affect the amplitude and frequency of [ Ca 2 + ] Cyt oscillations, but the mechanism is different. The model provides a comprehensive mechanism for predicting cytosolic 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.


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 Ca 2+ concentration ([Ca 2+ ] Cyt ) can be divided into the following two pathways: (i) the release of Ca 2+ from the intracellular Ca 2+ stores, mainly the endoplasmic reticulum (ER, the largest Ca 2+ store), or (ii) extracellular Ca 2+ influx to cytosol (Cyt), through the opening of plasma membrane Ca 2+ channels. It has been widely accepted that Ca 2+ release-activated Ca 2+ (CARC) channels are the main mode of 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 Ca 2+ influx and endoplasmic reticulum Ca 2+ supplementation. IP 3 interacts with Ca 2+ channels in the ER, causing the release of stored Ca 2+ , and the depletion of Ca 2+ in the ER triggers Ca 2+ entry through CRAC channels. CRAC channels in MCs are non-voltage-gated and show a characteristic inward rectification [7]. Therefore, the Ca 2+ flow of CRAC channels is related to the Ca 2+ concentration of the ER and Cyt. In contrast to the CRAC channels, the plasma membrane Ca 2+ -ATPase (PMCA) channels extrude Ca 2+ to the extracellular space, to maintain calcium concentration balance. They can transfer 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 Ca 2+ concentration increasing in the mitochondrial matrix ([Ca 2+ ] Mt ) [10]. Moreover, recent studies have shown that the Mt regulate Ca 2+ oscillations by firstly uptaking 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 Ca 2+ concentration will increase, then leukotriene C 4 (LTC 4 ) will be produced, which can activate phospholipase C (PLC), to promote PIP 2 decomposition to IP 3 [19]. Subsequently, IP 3 binds to the inositol 1,4,5trisphosphate receptor (IP 3 R) on the ER membrane, which induces Ca 2+ release from the ER. IP 3 R is regulated by Ca 2+ in a biphasic manner (stimulatory at low levels/inhibitory at high levels) [19][20][21][22]. The mitochondrial Ca 2+ uniporter (MCU) in mitochondria uptakes Ca 2+ quickly when exposed to high Ca 2+ concentration environments around the opened IP 3 R channel pore. MCs are activated by mechanical stimulation, and the 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 Ca 2+ uptake [23][24][25]. This explains why high Ca 2+ is observed when global 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 Ca 2+ dynamic model proposed by Othmer-Tang et al., Falcke et al. [31] added the mitochondrial 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 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 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 [Ca 2+ ] 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 Ca 2+ channels and Ca 2+ buffering in these compartments are considered. The functional region formed by the ER and Mt (µd) is explicitly assumed as a Ca 2+ pool. The degradation and production of IP 3 is added to this model, to investigate the effect of IP 3 dynamics on [Ca 2+ ] Cyt oscillations.

Mathematical Model
The full model includes plasma membrane channels, and degradation and production of IP 3 , Cyt, ER, Mt, and µd. MCs will release IP 3 after they are activated by mechanical stimuli. Then, IP 3 bines to IP 3 R to trigger the intracellular Ca 2+ signal. The whole progress is shown in Figure 1. Calcium dynamics in each compartment are governed by a balance of Ca 2+ fluxes, leaks, and buffering processes. [Ca ] oscillations.

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

Cross-Membrane Ca 2+ Current
According to previous researches, we accept that CRAC channels and plasma membrane Ca 2+ -ATPase (PMCA) channels are the main Ca 2+ channels in MCs [33]. CRAC channels are Ca 2+ influx channels, and PMCA channels are Ca 2+ outflux channels. The CRAC current is given by the Hodgkin-Huxley (HH) model [35], as follows: where g CRAC is the conductance, E m is the membrane potential, and is the Nernst potential for Ca 2+ , φ = RT zF , where R is the universal gas constant, T is the absolute temperature, z = 2 is the valence of Ca 2+ , and F is the Faraday constant. P CRAC is the proportion of CRAC channels in open state, and it is assumed as follows [35]: The PMCA current is given by the following [36]: where I PMCA,M is the maximum PMCA current, and K PMCA is the Ca 2+ concentration for the half activation of PMCA channels.

Ca 2+ Outflows from ER
The calcium outflow from the ER to Cyt or µd, through IP 3 R channels, is defined as follows: where V IP3R is the maximum total flux through IP 3 R channels [33], and P oIP3R and P oIP3R µd are the open probabilities of IP 3 R channels facing the Cyt and µd, respectively, and they are defined as follows [33]: where S act and S act µd express the probability of the activated subunit, respectively, and are defined by sigmoidal functions of [IP 3 ] and [Ca 2+ ] Cyt , and h is the slow inactivation gating variable, defined as follows: for µd, it is the following: where a 2 , d 1 , d 2 , d 3 , and d 5 are parameters. For IP 3 dynamics, the production speed of IP 3 is related to PLC, and the production of phospholipase C isoforms depends on [Ca 2+ ] Cyt , so the production speed of IP 3 is defined as follows: for µd, it is as follows: where V PLC is the maximal production rate of PLC isoforms, and K PLC is the sensitivity of PLC to Ca 2+ . IP 3 is degraded through phosphorylation by IP 3 kinases. The kinetic equation can be written as follows: for µd, it as follows: where k deg represents the phosphorylation rate constants, and K deg is the half-saturation constant of IP 3 kinases. IP 3 leaks from µd to Cyt can be defined as follows: Therefore, the change in [IP 3 ] Cyt and [IP 3 ] µd can be written as follows: SERCA pumps transport Ca 2+ into the ER, and the flux from Cyt to the ER is defined as follows: the flux from µd to the ER is defined as follows: where V SERCA is the maximal flux through SERCA, and k SERCA is the Ca 2+ activation constant for SERCA. 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: and the leak from the ER into the µd is defined as follows: Although the µd is not a membrane-bound compartment, we similarly define the leak from the µd into the Cyt as follows:

Ca 2+ Outflows from Mt
In the Mt, the mNCX channels exchange one Ca 2+ for three Na + . Flux through mNCX channels facing the Cyt is defined as follows: and for mNCX channels facing the µd, it is as follows: where [Na + ] Cyt and [Na + ] µd are the concentration of Na + in Cyt and µd, respectively. V mNCX is the maximal flux through the mNCX, and k Na and k mNCX are Na + and Ca 2+ activation constants for mNCX, respectively. The connectivity coefficient C mNCX is the proportion of mNCX channels facing the µd.
The MCU channel transports Ca 2+ into the Mt. Flux through the MCU to the Cyt is defined as follows: (29) and to the µd, it is as follows: where . V MCU 0 represents the maximal flux through the MCU, and ∆Φ is the voltage driving force. Ψ is the inner mitochondrial membrane voltage (150~180 mV, negative inside). b and Ψ 0 are the fitting parameters obtained from Qi [33]. During the simulations, we assume a constant Ψ of 170 mV, as experimental evidence suggests that Ψ does not change significantly in response to transient cytosolic [Ca 2+ ] increase, produced by IP 3 -generating agonists [37][38][39]. k MCU is the Ca 2+ activation constant for MCU, and the connectivity coefficient C MCU is the proportion of MCU channels facing the µd.

Effective Cytosol
Since the µd is a part of the cytosol, we defined the [Ca 2+ ] of the effective cytosolic Cyt as the volume-weighted average of [Ca 2+ ] within the combined Cyt and µd compartments.

µ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: where SA is the surface area, N is the mitochondria total number, and D is the ER-Mt distance.

Temporal Changes in [Ca 2+ ] in Each Compartment
Temporal changes in [Ca 2+ ] in each compartment are represented as the following ordinary differential equations: in cytosol, it is the following: in the ER, it is the following: in the Mt, it is the following: in µd, it is the following: where θ i (i = Cyt, ER, Mt, µd) is the buffer factor of each compartment, which is defined as follows [36]: The parameter values related to our model are given in Table 1.

Effect of the Degradation and Production of IP 3 on Ca 2+ 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 Ca 2+ concentration at different k deg , V PLC , and 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 Ca 2+ concentration is in the state of oscillation, and numerical oscillation analysis diagrams show that one k deg or V PLC value corresponds to two Ca 2+ concentrations. If the maximum value and minimum value are the same, it indicates that the Ca 2+ concentration is in equilibrium state, and numerical oscillation analysis diagrams show that one k deg or V PLC value corresponds to one Ca 2+ concentration.
The numerical oscillation analysis diagrams show the effect of Mt and µd compart- eff Cyt oscillatory dynamics, as functions of k deg and V 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 deg (range 4 to 9 s −1 ) (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 IP 3 oscillations region in numerical oscillation analysis diagrams. k deg is positively correlated with IP 3 degradation rate, and high levels of k deg mean that IP 3   [Ca ] 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,  eff Cyt 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, 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 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 PLC increases to a certain value, the oscillations disappear. This means that V 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.  [Ca ] in other cellular compartments are at peak value, as shown in Figure   5a. This indicates that the main 2+ Ca filling of the Mt and μd come from the ER. The numerical simulation results also prove that high     . (a,b) show that the oscillations are out of sync, but the frequencies are pretty much the same. [Ca ] oscillations increased slightly. The numerical simulation results of Qi [33] show that with increasing ER-Mt distance at D < 20 nm, the [Ca ] amplitude is the highest in Figure 6b.

Effect of the ER-Mt Distance (D) on Ca 2+ Oscillations
Moreover, Figure 6d shows that with ER-Mt distance increases at D < 65 nm, the   . (a,b) show that the oscillations are out of sync, but the frequencies are pretty much the same. oscillations is weak. This is a mathematical explanation of why the μd functional region should have a small distance. Figure 6e shows that (e) (f)    Figure 7 shows that the period and distance have an approximate linear relationship; therefore, we can obtain the slopes of [Ca 2+ ] Cyt (0.01244) and [Ca 2+ ] µ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.

Effect of the [IP ] Cyt on Ca 2+ Oscillations
The numerical oscillation analysis diagrams in Figure 8a show that the Mt, μd, and

Effect of the [IP 3 ] Cyt on Ca 2+ Oscillations
The numerical oscillation analysis diagrams in Figure 8a show that the Mt, µd, and   oscillations region of the model, with Mt, μd, and Mem, shrink obviously, meaning that the μd and Mem compartments limit the range of 2+ eff Cyt [Ca ] 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.

Discussion and Conclusions
In Qi's model [33], the IP 3 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 Ca 2+ between IP 3 R and MCU. However, with the distance increases, the influence of Ca 2+ diffusion in other directions will be more significant. Hence, the one-dimensional diffusion assumption will make the calculation error of [Ca 2+ ] µ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, [Ca 2+ ] µ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 Ca 2+ . This is because IP 3 R has a Ca 2+ inhibition binding site, so when Ca 2+ binds to this site, IP 3 activity is inhibited. As the ER-Mt distance increases, Ca 2+ can spread faster; therefore, the probability of Ca 2+ binding to the inhibition site decreases, so Ca 2+ released by IP 3 increases. Hence, there is a little upward trend of [Ca 2+ ] µd , and, later, the effect of the increased volume of µd is greater than the release of Ca 2+ , and [Ca 2+ ] µd begins to decrease. However, in Qi's result, [Ca 2+ ] µd was decreasing all the time. Arash Moshkforoush's model [34] does not consider the degradation and production of IP 3 ; therefore, the concentration of IP 3 in the cells is constant, which is not a physiological reality. IP 3 dynamic behaviors have a significant effect on the range of parameter values and the oscillation patterns of [Ca 2+ ] Cyt oscillations [43]. Inhibition of protein kinase C eliminates Ca 2+ oscillations, while 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 IP 3 , there are different amplitudes and frequencies of [Ca 2+ ] 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), Ca 2+ channels and buffer in these compartments, and the µd composed of the ER and Mt. In our simulations,  Figure 4b shows that the [Ca 2+ ] eff Cyt 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 [Ca 2+ ] Cyt oscillations increases with the increase in the ER-Mt distance. Moreover, the [Ca 2+ ] µ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 [Ca 2+ ] 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 [Ca 2+ ] Cyt and [Ca 2+ ] µd oscillations are the same at different ER-Mt distances, and they increase with the ER-Mt distance. Meanwhile, from Figures 5a and 7, we can understand that [Ca 2+ ] µd is 20 times higher than [Ca 2+ ] Cyt , and their oscillations are not synchronous. The presence of µd causes the ER to release less 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 Ca 2+ from the ER. All these results suggest that µd plays an important role in controlling [Ca 2+ ] Cyt oscillations at D < 65 nm. The degradation and production of IP 3 can also regulate [Ca 2+ ] Cyt oscillations by maintaining different levels of [IP 3 ] Cyt . As shown in Figure 5b, before stimulation, [IP 3 ] Cyt and [IP 3 ] µd get closer, due to the diffusion between Cyt and µd. This is because the IP 3 production and degradation of µd and Cyt maintain dynamic equilibrium at the initial moment, but the concentration of IP 3 will affect [Ca 2+ ] Cyt oscillations; therefore, we assume that the production and degradation of IP 3 are both zero before the stimulation (50 s), to reduce the impact of IP 3 on calcium oscillations in the study of V PLC and k 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 IP 3 also has an effect on oscillations, hence the given concentration of 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 IP 3 is assumed to be zero, and then the values of V PLC and k deg are restored, so as to reduce the impact of IP 3 on calcium oscillations in the study of V PLC and k 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. IP 3 R is regulated by Ca 2+ in a biphasic manner; therefore, IP 3 R activity is inhibited at higher Ca 2+ concentrations. When the ER-Mt distance increases in our model, the volume of µd increases and the Ca 2+ concentration of µd gradually decreases, and this causes the activity of IP 3 R to increase; therefore, more Ca 2+ will outflow from the ER through the IP 3 R channels, causing the ER calcium oscillations amplitude increases shown in Figure 6e, and the cytoplasmic 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 [Ca 2+ ] Cyt oscillations. Moreover, µd can play a critical role in 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 [Ca 2+ ] 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.