Optimal microdomain crosstalk between endoplasmic reticulum and mitochondria for Ca²⁺ oscillations

Abstract

A Ca²⁺ signaling model is proposed to consider the crosstalk of Ca²⁺ ions between endoplasmic reticulum (ER) and mitochondria within microdomains around inositol 1, 4, 5-trisphosphate receptors (IP₃R) and the mitochondrial Ca²⁺ uniporter (MCU). Our model predicts that there is a critical IP₃R-MCU distance at which 50% of the ER-released Ca²⁺ is taken up by mitochondria and that mitochondria modulate Ca²⁺ signals differently when outside of this critical distance. This study highlights the importance of the IP₃R-MCU distance on Ca²⁺ signaling dynamics. The model predicts that when MCU are too closely associated with IP₃Rs, the enhanced mitochondrial Ca²⁺ uptake will produce an increase of cytosolic Ca²⁺ spike amplitude. Notably, the model demonstrates the existence of an optimal IP₃R-MCU distance (30–85 nm) for effective Ca²⁺ transfer and the successful generation of Ca²⁺ signals in healthy cells. We suggest that the space between the inner and outer mitochondria membranes provides a defense mechanism against occurrences of high [Ca²⁺]_Cyt. Our results also hint at a possible pathological mechanism in which abnormally high [Ca²⁺]_Cyt arises when the IP₃R-MCU distance is in excess of the optimal range.

Introduction

The calcium ion (Ca²⁺) is a ubiquitous intracellular signal controlling diverse cellular functions, such as muscle cell contraction, neurotransmitter release from neurons and astrocytes, metabolic processes, egg activation and cell maturation, differentiation and death¹. Ca²⁺ signals commonly appear as repetitive spikes of cytosolic Ca²⁺ concentration ([Ca²⁺]_Cyt), with signaling information encoded in the frequency, amplitude and duration of these oscillations². The complex spatial-temporal profiles of [Ca²⁺]_Cyt depend heavily on the endoplasmic reticulum (ER) and mitochondria³ which act as the two major intracellular Ca²⁺ stores.

While it has been widely accepted that the inositol 1, 4, 5-trisphosphate receptors (IP₃R) function as Ca²⁺ release channels on the ER membrane¹, mitochondria were previously considered only to function as Ca²⁺-sinks under pathological conditions when there is an abnormally high cellular Ca²⁺¹. This view on mitochondria began to change in the 1990s with two landmark studies demonstrating that non-pathological elevations of [Ca²⁺]_Cyt are accompanied by a marked rise of mitochondrial matrix Ca²⁺ concentrations ([Ca²⁺]_Mt)⁴ and that [Ca²⁺]_Mt responds dynamically to physiological oscillations of [Ca²⁺]_Cyt⁵. Experimental observations indicated that mitochondria can accumulate Ca²⁺ through the mitochondrial Ca²⁺ uniporter (MCU) which is located on the inner mitochondrial membrane (IMM)⁶. In despite of its low Ca²⁺ affinity, the MCU is exposed to a high [Ca²⁺]_Cyt generated in microdomain around the channel pore of open IP₃Rs^7,8, thus even when global [Ca²⁺]_Cyt is low a high value of [Ca²⁺]_Mt may be observed. As a result, mitochondria can serve as buffers of cytoplasmic Ca²⁺ under both pathological and non-pathological conditions.

Although it is now generally believed that mitochondria function as a dynamic Ca²⁺ sequestration system to shape cytosolic Ca²⁺ oscillations and therefore regulate many physiological processes⁶, the detailed mechanism underlying how mitochondria influence global Ca²⁺ signals remains elusive. Contradictory mechanisms have been suggested on how mitochondrial Ca²⁺ uptake regulates IP₃R activity, including up-regulation^9,10,11, down-regulation^12,13 and complete independence¹⁴. Furthermore, although mitochondrial Ca²⁺ uptake via MCU depends primarily on the spacing between the ER and mitochondria¹⁵, the working distance between the two organelles remains unclear due to the complexity of the intracellular environment and the lack of spatial resolution in the experimental observations. Indeed, the existing estimates of the distances between the ER and mitochondria vary hugely, from less than 10 nm to more than 200 nm¹⁶. A recent paper published in PNAS¹⁷ even argued that the mitochondria may not act as a significant dynamic buffer of cytosolic Ca²⁺ under physiological conditions.

Computational models provide a valuable tool for understanding of the mechanisms underlying Ca²⁺ signal transduction¹⁸ and the dynamics of Ca²⁺ release from the ER has been extensively modeled^{19,20,21,22,23}. Some models only focused on mitochondrial Ca²⁺ dynamics^24,25. Dash and co-workers developed a detailed kinetic model for MCU to match published data sets on mitochondrial Ca²⁺ uptake²⁵. However, the involvement of both the ER and mitochondria in Ca²⁺ dynamics has only been considered in a limited number of models^{18,26,27,28,29,30}. In the earlier models^18,26, it was postulated that the MCU substantially sequesters Ca²⁺ when [Ca²⁺]_Cyt < 1 μM, in contrast to experimental evidence that the MCU actually sequesters Ca²⁺ in the range of 10–20 μM^3,31. Dupont et al.²⁸ modeled the effect of Hint2, a mitochondrial protein, on cytosolic Ca²⁺ dynamics in hepatocytes. However, in these models the maximum value of ER luminal [Ca²⁺] ([Ca²⁺]_ER) is in the order of several to several tens of μM, deviating drastically from most of the measured figures of 100–900 μM³² and the value of [Ca²⁺]_Mt was also significantly underestimated. More importantly, except the model proposed by Szopa et al.²⁹, none of them take into account the high [Ca²⁺]_Cyt microdomain, which is critical for mitochondrial Ca²⁺ signaling. The model considered by Szopa et al.²⁹ is a modification of the model from ref. 26, analyzing the influence of microdomain on the period and shape of calcium oscillations. They supposed that MCU sense elevated Ca²⁺ concentration which is directly equal to that in ER²⁹. This assumption is incoherent with the experimental observation that [Ca²⁺]_Cyt in microdomain is about several tens of μM, which is about 10-fold higher than that in the bulk cytosol⁷ and about 10-fold lower than that in ER³². Thus, there is a lack of model to quantitatively investigate how the crosstalk between the ER and mitochondria controls Ca²⁺ signaling.

Here we consider a more realistic model to investigate the role of mitochondria in Ca²⁺ signaling to couple the ER and mitochondria based on the Ca²⁺ microdomain. A cytosolic microdomain is specifically considered in order to discuss the MCU dynamics and so each MCU specifically responds to the local high Ca²⁺ concentration in microdomain generated by a cluster of IP₃Rs. Our results demonstrate the critical role of mitochondrial Ca²⁺ uptake in modulating IP₃R-released Ca²⁺ signaling. We show that the location of the MCU relative to IP₃Rs is a key determinant for modulating Ca²⁺ signals. There is a critical IP₃R-MCU distance at which 50% of the ER-released Ca²⁺ is taken up by mitochondria and that mitochondria modulate Ca²⁺ signals differently when outside of this critical distance. When the IP₃R-MCU distance is greater than the critical distance mitochondrial Ca²⁺ uptake stimulates ER-Ca²⁺ release by rapidly reducing the amplitude of the [Ca²⁺]_Cyt signal within the microdomain and suppressing the inhibition dynamics of high [Ca²⁺]_Cyt on IP₃R. However, when the distance is less than the critical distance, mitochondrial Ca²⁺ uptake demonstrates a different effect on IP₃R dynamics. Initially strong mitochondrial uptake forces [Ca²⁺]_Cyt to a low level, thus preventing the strong activation dynamics of [Ca²⁺]_Cyt on IP₃R. This leads to strong IP₃R inhibition on the decay of subsequent [Ca²⁺]_Cyt spikes. Thus, our model predicts that when MCU are too closely associated with IP₃Rs, the enhanced mitochondrial Ca²⁺ uptake will produce an increase of cytosolic Ca²⁺ spike amplitude. As a result, there is an optimal spacing about 30–85 nm between the ER and mitochondria for generation of experimentally observed Ca²⁺ oscillations in living cells.

Modelling coupled ER and mitochondria Ca²⁺ dynamics

A “closed” cell model, which does not include Ca²⁺ exchange between cytosol and external medium through the plasma membrane, is considered to investigate the coupling effects between the ER and mitochondria on Ca²⁺ signaling (Fig. 1(a)). As a result, the model has three compartments: cytosol (Cyt), ER and mitochondria (Mt). As shown in Fig. 1(a), the dynamics of the free [Ca²⁺] in these compartments is determined by the Ca²⁺ fluxes from various channels and pumps and by buffering processes with various Ca²⁺ binding proteins (BP), given by,

Figure 1

The model design.

(a) The schematic diagram of the components and fluxes included in the model. When a cell is stimulated to generate IP₃, IP₃ binds to IP₃R and triggers Ca²⁺ release from the ER through IP₃Rs. The Ca²⁺ in the cytsol can be pumped back into ER by SERCA. (b) IP₃Rs and MCUs in close proximity form a microdomain, a half 3-dimensional cytosolic space separated by ER membrane. A high microdomain [Ca²⁺] ([Ca²⁺]_Mic) can be generated upon opening of IP₃Rs which in turn leads to Ca²⁺ uptake into mitochondria through MCUs. Ca²⁺ is extruded from mitochondria back into the cytosol by NCX. (c) The four-state model of MCU monomer with binding and unbinding rates.

Here [Ca²⁺] represents the free [Ca²⁺], V the volume of three compartments, jⁱⁿ the Ca²⁺ flux from outside to inside and j^out the opposite flux. The three terms j_Ca represent Ca²⁺ dissociated from the BP.

In each compartment, Ca²⁺ can be buffered by various BPs. The kinetics between free Ca²⁺ and BP is described by a simple chemical reaction

giving

where [BP]_{Tot_i} and [CaBP]_i represent the total concentration of BP and the concentration of Ca²⁺-bound BP in each compartment, k_on and k_off denote the on and off rate constants of Ca²⁺ with the BP, respectively.

Ca²⁺ fluxes in ER component

The Ca²⁺ effluxes from the ER to cytosol are given by the IP₃R release and a passive leakage with the concentration gradient across the ER membrane as their driving force, i.e.,

with v_IP3R being the maximal efflux from the IP₃Rs in the ER membrane, P_oIP3R the open fraction of IP₃Rs and c_Leak the Ca²⁺ leak constant.

The IP₃R channel is an assembly of four equivalent subunits, each of which is mutually gated by IP₃ and Ca²⁺. The gating dynamics of each IP₃R subunit is described by the Li-Rinzel model¹⁹ and we assumed that the channel is open if at least three of its subunits are in the active state²⁰. Thus, we get the following expression for P_oIP3R:

Here s_act is the probability in the active state of the subunit, [IP₃] is the IP₃ concentration and the slowly gating variable h is given by

where a_i and d_i are constants.

The cytosolic Ca²⁺ can be pumped back into ER by active SERCA which is given by,

where v_SERCA and k_SERCA are maximum SERCA pump flux and SERCA activation constants.

Ca²⁺ fluxes in mitochondrial component

Mitochondrial Ca²⁺ uptake through MCU is driven primarily by the large voltage across IMM. Besides the IMM voltage, it has been shown that MCU flux is also regulated by MICU³³. For simplicity, in our model we only consider the role of MICU1 which has been well investigated in experiment. Then the MCU influx can be expressed as

Here v_MCU is the maximal flux of MCU, ΔΦ is the voltage driving force, R_MICU represents the MICU1 regulatory term and P_oMCU describes the open fraction of MCU.

The IMM voltage Ψ is 150 ~ 180 mV (negative inside). Although the net Ca²⁺ flux into the mitochondria might associate with a reduction in voltage, some evidence suggests that mitochondrial Ca²⁺ uptake does not appear to alter Ψ significantly during physiological condition³⁴. With the assumption of constant Ψ, the driving force for MCU can be given by³⁵

Here F, R and T are the Faraday constant, the gas constant and the Kelvin temperature, respectively; b and Ψ₀ are fitting parameters³⁵.

As a regulatory subunit, MICU1 senses [Ca²⁺] in mitochondrial intermembrane space to set a threshold for MCU³³. The influx rate of MCU only becomes substantial when the extra-mitochondrial [Ca²⁺] reaches values above 3 μM³⁶. As a result, the regulation of MICU1 to MCU activity has a sharp activation threshold (k_MICU) with an apparent Hill coefficient equal to 4³⁶:

where [Ca²⁺]_Mic is the microdomain [Ca²⁺]_Cyt that MCU is exposed to, which will be discussed later.

Whereas MICU1 acts as a regulator, MCU represents the core pore-forming component responsible for mitochondrial Ca²⁺ uptake^37,38. Recent experimental and bioinformatic analysis reveals that MCU is an oligomer³⁷ and most likely a tetramer³⁹. Considering the biphasic dependence of MCU on [Ca²⁺]_Cyt, i.e., fast Ca²⁺-dependent activation and slow Ca²⁺-dependent inactivation of the MCU³¹, we simply propose that each monomer has two Ca²⁺ binding sites, one for activation and another for inhibition. Therefore a monomer can exist in 4 states of s₀₀, s₁₀, s₀₁ and s₁₁ with the first index for activation binding site and the second index for inhibition binding site and with 1 for binding and 0 unbinding (Fig. 1(c)). Similar as the IP₃R, we posit that the MCU is open if either 3 or 4 of the 4 monomers are in the active state (s₁₀). Thus, the MCU open fraction is given by,

where x_ij denotes the probability of a monomer in states s_ij. The dynamical equations governing the state of a monomer are

with x₀₀ + x₀₁ + x₁₀ + x₁₁ = 1.

Ca²⁺ is extruded from mitochondria primarily by NCX which swaps 3 Na⁺ ions with 1 Ca²⁺ ion during each exchange cycle, Thus, the kinetic of outward calcium flux can be described as³⁵

where [Na⁺]_Cyt is Na⁺ concentration in the cytosol, v_NCX is maximal NCX activity, both k_Na and k_NCX are activation constants for NCX.

Microdomain crosstalk between the ER and mitochondria

Now we describe a microdomain model that represents the crosstalk occurring between closely situated ER and mitochondria (Fig. 1(b)). The interacting players are IP₃Rs and MCUs in the microdomain. A microdomain of local high [Ca²⁺]_Cyt produced by nearby open IP₃Rs is essential for high capacity Ca²⁺ transport into mitochondria via MCU. It has been reported that 90% of IP₃Rs are located close to mitochondria and the majority of MCU sense the high [Ca²⁺]_Mic in the vicinity of IP₃Rs²⁹. Here we make a specific assumption that MCU senses the local high [Ca²⁺]_Mic generated by the opening of clustered IP₃Rs, but other Ca²⁺ handling components, such as IP₃Rs, SERCA, NCX, etc., are determined by the global averaged [Ca²⁺]_Cyt. The pros and cons of our method will be discussed in the last section of the paper.

It has been suggested that each mitochondrial Ca²⁺ uptake site must face multiple IP₃Rs for effective mitochondrial Ca²⁺ uptake⁴⁰. Recent experiment evidence indicates that IP₃Rs are typically distributed in clusters on the ER membrane, with each cluster containing a handful of IP₃Rs^41,42. Patch clamp experiments also show that the Ca²⁺ current (I_sIP3R) passing through a single open IP₃R is linearly correlated with the concentration gradient ([Ca²⁺]_ER - [Ca²⁺]_Cyt) across the ER membrane with a slope of S_IP3R = 0.3 pA/mM⁴³. Considering the closely distributed IP₃Rs, we treat the Ca²⁺ flux from clustered IP₃Rs as a Ca²⁺ current from a point source at cluster center. As a result, an effective Ca²⁺ current from clustered IP₃Rs is expressed as

Here, P_oIP3R is the IP₃R open probability given by Eq. 6 and n_IP3R is the IP₃R number in a cluster. According to Ref. 41, we consider a typical cluster with n_IP3R = 4. These 4 IP₃Rs are in tight association, organized in a square pattern.

Around the center of IP₃R cluster, the open IP₃Rs will generate a very sharp distribution of [Ca²⁺]_Cyt due to rapid Ca²⁺ diffusion producing a microdomain of high [Ca²⁺]_Cyt. Here we also assume that the ER membrane acts as a flat boundary to limit the Ca²⁺ diffusion in a half 3-dimensional space, i.e. the cytosolic microdomian (Fig. 1(b)). Using hemispherical symmetry of Ca²⁺ diffusion from a point release source, the sharp distribution around the release point can be established with a signal equation analogous to the linear cable equation^44,45. As a solution of the linearized reaction-diffusion equation, the microdomain [Ca²⁺] is given by⁴⁵

with

Here, D_c is the diffusion coefficient for free Ca²⁺. The parameter λ is an important factor that determines the sharp decay of [Ca²⁺] within a microdomain, which accounts for the Ca²⁺ diffusion and binding to BPs. [BP]_Mic is the total concentration of BP in the microdomain and K_d is its dissociation constant^44,45. One can see that the [Ca²⁺]_Mic within a microdomain becomes the cellular average [Ca²⁺]_Cyt when there are no open IP₃Rs.

Eq. 16 describes the [Ca²⁺] distribution at distance r around the point source, i.e., the IP₃R cluster center. In our model the parameter r is the distance between the MCU and the IP₃R cluster center, which we will call the IP₃R-MCU distance in the paper. Thus, [Ca²⁺]_Mic is the microdomain Ca²⁺ concentration that the MCU sees. Considering MCUs located on IMM, here we also simply assume that the outer mitochondrial membrane (OMM) has little effect on the Ca²⁺ diffusion as it has high permeability to Ca^{2+ 6}. If we suppose the IMM and ER membrane are parallel, parameter r also represents the distance between IMM and ER membrane. As a result, Eq. 16 represents the microdomain coupling between the ER and mitochondria in our model.

Parameter values

The parameters used in this model are listed in Table 1 (see supplementary material). They are categorized on the basis of their source to: (a) those extracted directly from experimental data; (b) those obtained from previously published models; or (c) the free parameters which were fit to experimental results. The model presented in this study permits simultaneous calculation of [Ca²⁺] changes in the 3 compartments. Currently, it is difficult to measure experimentally Ca²⁺ concentrations in the cytosol and ER or mitochondria simultaneously because of the highly dynamic nature of the Ca²⁺ signals that occur within organelles⁴⁶. Only a few existing experimental data sets^46,47 consist of measurements from all 3 compartments. The free parameters used in our model were extrapolated from fits to these data sets.

Results

The Ca²⁺ oscillations

Our first priority was to determine whether the model realistically simulated cellular Ca²⁺ dynamics. As shown in Fig. 2a, when setting [IP₃] = 0.5 μM our model generates baseline Ca²⁺ oscillations occurring as a series of discrete spikes. The values for basal and stimulated [Ca²⁺] for the cytosol, mitochondria and ER in the model are in close accord with experimental data^38,48. The oscillations period is in the order of several ten seconds, as seen in many experimental observations with non-exitable cells^2,48. The Ca²⁺ spikes show relatively rapid rising phases and slow falling phases, mimicking those seen in experimental recordings⁴⁹. The cytosolic and mitochondrial [Ca²⁺] oscillate almost in phase, with [Ca²⁺]_Cyt as the leader, while the ER [Ca²⁺] is in anti-phase with them, also as observed experimentally^13,47.

Figure 2

Mitochondria serve as Ca²⁺ reservoirs.

(a) Oscillations of: [Ca²⁺]_Mic in the IP₃R-MCU microdomain, [Ca²⁺]_Cyt in the cytosol, [Ca²⁺]_Mt in mitochondria and [Ca²⁺]_ER in the ER in the presence of mitochondria with a IP₃R-MCU distance of 35 nm. The width at half maximum, the rising duration and decay duration of [Ca²⁺]_Cyt spike are marked in the second panel of (a). (b) Oscillations of: [Ca²⁺]_Cyt in the cytosol and [Ca²⁺]_ER in the ER in the absence of mitochondria. [IP₃] = 0.5 μM. The unit of [Ca²⁺] is μM.

Having established that the salient features of the Ca²⁺ dynamics of our model are in good agreement with existing experimental observations, we examined the effect of manipulating the mechanisms modulated by the crosstalk of the ER and mitochondria. Thus, we also performed simulations in which mitochondria were absent from the model (Fig. 2b). The minimal values of [Ca²⁺]_ER are 311 μM and 276 μM in the absence and in the presence of mitochondria, respectively (Fig. 2), indicating that more Ca²⁺ ions are released from the ER during each spiking cycle in the presence of mitochondria. The maximal values of [Ca²⁺]_Cyt are 5.6 μM and 2.5 μM in the absence and in the presence of mitochondria, respectively (Fig. 2), showing that mitochondria can significantly decrease [Ca²⁺]_Cyt oscillation amplitude.

Ca²⁺ oscillation modulated by IP₃R-MCU distance

Next we looked at the effect of IP₃R-MCU distance on Ca²⁺ oscillation in detail. Fig. 3 compares oscillating [Ca²⁺] trajectories in microdomain, ER, mitochondria and cytosol during a single cycle at varying IP₃R-MCU distances. By decreasing the IP₃R-MCU distance from infinity (we used simulations in which mitochondria were absent from the model as the infinite IP₃R-MCU distance), the [Ca²⁺]_Mic in microdomain increases according to Eq. 16 against distance (Fig. 3a), giving an increase in mitochondrial Ca²⁺ uptake. As a result, the minimal [Ca²⁺]_ER steadily decreases (Fig. 3b) and the [Ca²⁺]_Mt spike amplitude steadily increases until saturating at about 65 μM (Fig. 3c). However, an unexpected finding is that the [Ca²⁺]_Cyt spike amplitude first decreases but then increases as the distance decreases (Fig. 3d).

Figure 3

The time evolution of [Ca²⁺]_Mic (a), [Ca²⁺]_ER (b), [Ca²⁺]_Mt (c) and [Ca²⁺]_Cyt (d) in an oscillating cycle at distance r = 200 (green), 100 (purple), 40 (red), 20 (dark) and 10 nm (grey). The brown lines in (b) and (d) are in the absence of mitochondria. In (a) and (b) the duration of an abrupt increase in Ca²⁺ release is marked by arrows for r = 10 nm. [IP₃] = 0.5 μM.

The critical IP₃R-MCU distance

The amplitudes of [Ca²⁺]_Mic, [Ca²⁺]_Mt and [Ca²⁺]_Cyt are plotted in detail as a function of distance from 10 to 150 nm at [IP₃] = 0.5 μM in Fig. 4a and 4b. A sharp increase is observed for [Ca²⁺]_Mic and [Ca²⁺]_Cyt at r < 30 nm and [Ca²⁺]_Cyt reaches its minimum at r = 30 nm. Fig. 4b also shows that the [Ca²⁺]_Cyt spiking period increases monotonically from 30 to 50 sec, following the increased ability of mitochondria to take up Ca²⁺ as the IP₃R-MCU distance diminishes. The [Ca²⁺]_Cyt spiking amplitude is determined by the number of Ca²⁺ ions that diffuse into the cytosol, which is related to the total amount of Ca²⁺ ions released from ER through IP₃Rs and the fraction which are taken by mitochondria. In the model, the total amount of Ca²⁺ released from the ER through IP₃R and the fraction accumulated by mitochondria through MCU during each oscillating cycle are plotted in Fig. 4c as a function of IP₃R-MCU distance. The fraction of mitochondrial Ca²⁺ uptake increases monotonically with decreasing distance indicating an increase in the ability of the mitochondria to take up Ca²⁺ at closer distances. Our simulation shows that at 30 nm distance, about 50% of the ER-released Ca²⁺ ions are taken up by mitochondria (see green line with stars in Fig. 4c). As a consequence, r_crt = 30 nm acts as a critical distance for mitochondria to modulate cytosolic Ca²⁺ signaling: At r > r_crt, more ER-released Ca²⁺ ions will diffuse into the cytosol than be taken up by mitochondria; whilst at r < r_crt, mitochondria will capture the majority of Ca²⁺ ions released from ER.

Figure 4

Ca²⁺ dynamics modulated by the IP₃R-MCU distance.

(a) The spike amplitudes of [Ca²⁺]_Mic (triangles) and [Ca²⁺]_Mt (circles). (b) The spike amplitudes of [Ca²⁺]_Cyt (diamonds) and spike period (spheres). (c) The total amount of Ca²⁺ released from the ER through IP₃Rs during each cycle (squares), the amount of Ca²⁺ released from the ER through IP₃Rs during the rising time (circles) and the fraction delivered to mitochondria through MCU (stars). (d) The spike width (triangles), the rising time (squares) and the decay time (circles) of [Ca²⁺]_Cyt oscillations. [IP₃] = 0.5 μM.

Ca²⁺ dynamics modulated by IP₃R-MCU distance

Surprisingly, Fig. 4c also shows another unexpected finding that the total amount of Ca²⁺ released from the ER first increases with distance down to r = 20 nm and then decreases at smaller distance (see black line with squares in Fig. 4c). This raises a critical question: with decreasing distance at r < 20 nm, how can [Ca²⁺]_Cyt spiking amplitudes increase (violet line with diamonds in Fig. 4b), while the total amount of Ca²⁺ released from the ER decreases (black line with squares in Fig. 4c) but mitochondrial Ca²⁺ uptake fraction increases (green line with stars in Fig. 4c)?

In order to further examine this conflicting effect of distance on [Ca²⁺]_Cyt spike amplitude and the amount of ER-released Ca²⁺, we investigated how mitochondrial Ca²⁺ uptake modulates the spike width, the spike rising time and decay time of [Ca²⁺]_Cyt at varying IP₃R-MCU distance. Here the spike width is defined as the duration at half maximum of [Ca²⁺]_Cyt spike amplitude, the rising time is defined as the duration from twice of the minimal [Ca²⁺]_Cyt to the peak of the response and the decay time as the duration from the peak of the response to twice of the minimal [Ca²⁺]_Cyt. An example of definition of these parameters is marked in Fig. 2a for spike curve at distance r = 35 nm.

We first discuss the simple situation at r > r_crt. As shown in Fig. 4d when r > r_crt, the spike width increases with decreasing distance, which primarily arises from an increase in the duration of the decay time, whereas the rising time is rapid and shows little change with distance. As r decreases, the ability of mitochondrial Ca²⁺ uptake increases, leading to a decrease in [Ca²⁺]_Cyt amplitude. The decrease in [Ca²⁺]_Cyt amplitude in turn leads to a weakening of the Ca²⁺ inhibition of IP₃R dynamics, resulting in the prolonged spike decay time with decreasing distance.

Differently, at r < r_crt, as the IP₃R-MCU distance decreases, even though the rising time increases (black line in Fig. 4d), the spike width actually becomes smaller due to a decrease in the decay time (blue and red lines in Fig. 4d). At this small distance, the rising time consists of both a slow rising phase and an abrupt rising phase. As an example, the slow rising and abrupt rising phases for a spike at distance r = 10 nm are indicated by arrows in Fig. 3a&b. The mitochondrial Ca²⁺ uptake is so strong that, at the beginning of each oscillating cycle, most of the ER-released Ca²⁺ ions are directly driven into mitochondria, leading to a rapid increase of [Ca²⁺]_Mt (Fig. 3c), but a very slow increase of [Ca²⁺]_Cyt (Fig. 3d). With low [Ca²⁺]_Cyt, there are a small fraction of activating IP₃Rs which remain open for a long time, resulting in a slow rise phase for [Ca²⁺]_Cyt (e.g. the gray dash line for r = 10 nm in Fig. 3d with time from 0 to 30 sec). Such a long duration with low [Ca²⁺]_Cyt actually contributes little to the spike width. When [Ca²⁺]_Mt increases to approach to its saturation state, most of the ER-released Ca²⁺ ions will then diffuse into cytosol. Because of the Ca²⁺ induced Ca²⁺ release, there occurs an abrupt increase in the numbers of open IP₃Rs, generating a spike (e.g. the gray dash line for r = 10 nm in Fig. 3d with time from 30 to 35 sec).

With shorter IP₃R-MCU distance, the mitochondria become more saturated prior to the abrupt rising phase and thus during the abrupt rising phase the Ca²⁺ ions released from the ER diffuse more into cytosol, generating a higher amplitude spike in [Ca²⁺]_Cyt (violet line with diamonds in Fig. 4b at r < 30 nm). In turn, the elevated [Ca²⁺]_Cyt peak induces a strong inhibitory effect on IP₃R dynamics, such that spikes decay faster as the IP₃R-MCU distance decreases (red line with circles in Fig. 4d).

The changes with distance seen in the amount of Ca²⁺ ions released from the ER during the rise and decay phases can then explain why the total amount of ER-released Ca²⁺ reaches a maximum around r = 20 nm. At r > r_crt, as the rising time is fast and little affected by distance (close to 6 s as shown in Fig. 4d), the elevation in the total amount of ER-released Ca²⁺ with decreasing distance mainly related to the increase of the decay phase. By decreasing r at r_crt > r > 20 nm, although the amount of ER-released Ca²⁺ ions during the decay phase decreases because of the shortened decay duration (red line in Fig. 4d), the amount of ER-released Ca²⁺ ions during the rising phase increases because of the prolonged rising duration (black line in Fig. 4d). As a result, the total amount of Ca²⁺ released by ER keeps increasing with decreasing r for r_crt > r > 20 nm. However, at very short distance of r < 20 nm, during the rising phase the ER-released Ca²⁺ ions almost saturate both mitochondria and the cytosol, giving a maximal Ca²⁺ amount released by ER during the rising phase (see the open circles in Fig. 4c at r < 20 nm). Then, the change of the total amount of ER-released Ca²⁺ is typically determined by the change of decay time, leading to a decrease in the total amount of ER-released Ca²⁺ with decreasing distance at r < 20 nm.

Thus we can see how the apparently paradoxical relationship between [Ca²⁺]_Cyt spike amplitude and total amount of ER-released Ca²⁺ at IP₃R-MCU distances <20 nm develops. At r < 20 nm, the [Ca²⁺]_Cyt spike amplitude is mainly determined by the amount of Ca²⁺ released during the abrupt rising phase of a Ca²⁺ response, giving an increasing relationship with decreasing distance. In contrast to this, the change in the total amount of ER-released Ca²⁺ mainly follows the behavior of the spike decay time, giving a decreasing relationship with decreasing distance at r < 20 nm.

The synchronization of the oscillation phases between [Ca²⁺]_Cyt, [Ca²⁺]_ER and [Ca²⁺]_Mt is also modulated by the IP₃R-MCU distance. Our data (not shown here) indicate that at r > r_crt, [Ca²⁺]_Cyt peaks first, followed by [Ca²⁺]_ER minimum and then [Ca²⁺]_Mt peak, due to an increase of the amount of ER-released Ca²⁺ ions diffusing into the cytosol at greater distances. However, when r < r_crt, more ER-released Ca²⁺ ions are sequestered into mitochondria, so [Ca²⁺]_Mt peaks first, followed by [Ca²⁺]_ER minimum and then [Ca²⁺]_Cyt peak.

Ca²⁺ dynamics modulated by [IP₃]

In the following discussion, we investigate how [IP₃] controlled [Ca²⁺]_Cyt oscillations can be modulated by the IP₃R-MCU distance. The bifurcation diagram for [Ca²⁺]_Cyt oscillation (i.e. the steady state of [Ca²⁺]_Cyt or the maximum and minimum of [Ca²⁺]_Cyt oscillation) with [IP₃] as the bifurcation parameter is plotted at different IP₃R-MCU distances in Fig. 5a. When IP₃ levels are low, [Ca²⁺]_Cyt can be seen to remain in a low stable steady state. Beyond a certain threshold value of [IP₃], [Ca²⁺]_Cyt undergoes periodical oscillations. We also found that, at distance range of 27–122 nm, the system only demonstrates sustained Ca²⁺ oscillations when [IP₃] values are in an intermediate range, evolving into a stable steady state at higher [IP₃] (Fig. 5b). Such a bifurcation diagram for r = 60 nm is depicted in the inset of Fig. 5b. Similar bifurcation behavior has been observed in the experimental Ca²⁺ waves in Xenopus laevis oocytes. In oocytes it has been shown that repetitive Ca²⁺ waves occur only at intermediate [IP₃] between 0.1 and 1 μM⁵⁰. Thus, as shown in the inset of Fig. 5b, r = 60 nm is the optimum IP₃R-MCU distance in which our model replicates these experimental observations. According to the experiment⁵⁰, for a biologically meaningful Ca²⁺ oscillation model there should be a right bifurcation point for high [IP₃], beyond which there occurs no Ca²⁺ oscillation. As plotted in Fig. 5b, only in the range of r of 27–122 nm the model can produce Ca²⁺ oscillations in an intermediate [IP₃] range. Outside of such distance range the Ca²⁺ oscillations always occur even at infinitely high [IP₃] (Fig. 5b).

Figure 5

Ca²⁺ dynamics modulated by [IP₃].

(a) The bifurcation diagram (the stable steady states and the maxima and minima of Ca²⁺ oscillations) of the model with [IP₃] as the bifurcation parameter at distance r = 200 (green), 100 (purple), 40 (red), 20 (dark) and 10 nm (grey). (b) The value of [IP₃] at the right bifurcation point plotted against distance. The inset is the bifurcation diagram at distance r = 60 nm, in which the arrow marks the right bifurcation point. (c) The maximal [Ca²⁺]_Cyt against the IP₃R-MCU distance for [IP₃] = 0.1 (lozenges), 0.5 (squares), 1.0 (circles), 2.0 (triangles) and 5.0 μM (stars). The dashed and dotted lines represent [Ca²⁺]_Cyt = 4 and 5 μM, respectively.

[Ca²⁺]_Cyt oscillations with limited amplitudes regulate a host of vital cell functions¹. A global [Ca²⁺]_Cyt less than 3.5 μM is usually observed in living cell^38,48. However, high [Ca²⁺]_Cyt will trigger and modulate apoptosis⁵¹. Fig. 5c shows the maximal [Ca²⁺]_Cyt values seen in our model plotted against IP₃R-MCU distance for different [IP₃]. One can see that a rapid elevation in [Ca²⁺]_Cyt amplitude is observed at IP₃R-MCU distances around 10 nm. A global [Ca²⁺]_Cyt which is smaller than 4.0 or 5.0 μM (i.e. dashed or dotted line in Fig. 5c) is typically obtained with IP₃R-MCU distances of 20 to 85 nm or distances of 15 to 150 nm, respectively.

Robustness with varying clustered IP₃R number

In above study, we only consider a cluster in which n_IP3R = 4 IP₃Rs correspond to a single MCU. Experimental evidence suggests that IP₃R clusters typically contain a handful (~4) of IP₃Rs^41,42. Next we discuss the robustness of the model's behavior upon slight changes in n_IP3R (from 4 to 3 and 5). By calculating the fraction of Ca²⁺ ions accumulated by mitochondria compared to the total amount of Ca²⁺ released from the ER, we see that the critical distance is changed to 20 nm for n_IP3R = 3 and 35 nm for n_IP3R = 5 (Fig. 6(a)). Fig. 6(b) gives the value of [IP₃] at the right bifurcation point against distance with varying n_IP3R. A common region to find the right bifurcation [IP₃] is obtained with distances of 33 to 93 nm. Fig. 6(c) shows the maximal [Ca²⁺]_Cyt plotted against IP₃R-MCU distances for varying n_IP3R. A global [Ca²⁺]_Cyt which is smaller than 4.0 or 5.0 μM (i.e. dashed or dotted line in Fig. 6(c)) is typically obtained with IP₃R-MCU distances of 25 to 85 nm or distances of 20 to 110 nm, respectively.

Figure 6

Robustness of Ca²⁺ dynamics on n_IP3R.

(a) The fraction of Ca²⁺ accumulated by mitochondria of the total amount of Ca²⁺ released from the ER during each oscillating cycle plotted against IP₃R-MCU distance. (b) The value of [IP₃] at the right bifurcation point against distance. (c) The maximal [Ca²⁺]_Cyt against the IP₃R-MCU distance at saturating [IP₃] = 5.0 μM. The dashed and dotted lines represent [Ca²⁺]_Cyt = 4 and 5 μM, respectively. Here n_IP3R = 3 (stars), 4 (circles) and 5 (triangles).

Taken together, compared with experimental observations^38,48,50,51, the simulation results shown in Fig. 5 and Fig. 6 suggest that the optimal IP₃R-MCU distance for physiological Ca²⁺ signaling should be around 30 to 85 nm. When the IP₃R-MCU distance is between 30 and 85 nm, as illustrated in Fig. 6(a) approximately 15–55% of ER-released Ca²⁺ ions can be taken up by mitochondria, compatible with experimental measurements made in various cells^52,53. Our model thus suggests that mitochondria, typically working at a short distance from the ER (but larger than r_crt), serve as Ca²⁺ sinks, sequestering Ca²⁺ ions released from the ER.

Discussion

Since the 1990s, the manner in which mitochondrial Ca²⁺ uptake through MCU shapes intracellular Ca²⁺ signaling has attracted much attention⁶. The key questions presently debated are: to what extent mitochondria acquire Ca²⁺; what impact mitochondria have on cytosolic Ca²⁺ signals; and what dimensions the ER/mitochondria Ca²⁺ microdomain might have. To address these questions, we constructed a Ca²⁺ signaling model to consider the Ca²⁺ crosstalk within microdomains between IP₃Rs and MCU. A 4-states MCU model has been suggested based on the latest experimental results^31,37,38,39. The MCU is exposed to the microdomain Ca²⁺ concentration which is released from the nearby clustered IP₃Rs^34,36, rather than the bulk cytosolic Ca²⁺ concentration. As a result, the model replicates the shape, amplitude, period and kinetics of the Ca²⁺ spikes observed in experimental data, allowing us to examine the effects of mitochondrial Ca²⁺ uptake on cytosolic Ca²⁺ signals.

From experimental observations, it has been suggested that mitochondrial Ca²⁺ uptake contributes significantly to the physiological regulation of cytosolic Ca²⁺³⁸. Regarding the recent controversy whether mitochondria act as significant buffers to cytosolic Ca²⁺ under physiological conditions¹⁷, our model indicates that mitochondria can accumulate large quantities of Ca²⁺ under physiological conditions (about 15–55% of the total released from the ER) at IP₃R-MCU distances of 30–85 nm. Furthermore, our model indicates that the distance between the ER and mitochondria appears to be an important factor in modulating Ca²⁺ signaling, especially considering that mitochondria can be highly mobile organelles⁵⁴. There is a critical IP₃R-MCU distance r_crt (around 20 to 35 nm depending on n_IP3R) at which half of the Ca²⁺ released by the ER is taken up by mitochondria. At distances smaller than r_crt mitochondria can take up more Ca²⁺, while at distances larger than r_crt, more Ca²⁺ ions diffuse from the ER into the cytosol. This leads to mitochondria displaying different kinds of Ca²⁺ signal modulation depending on their distance from the ER. The effects of mitochondria on the spike amplitude, spike width and rising and decay time of [Ca²⁺]_Cyt spikes, the synchronization phase among [Ca²⁺]_Cyt, [Ca²⁺]_ER and [Ca²⁺]_Mt and the total amount of Ca²⁺ released from the ER are all substantially modulated by how the IP₃R-MCU distance deviates from this critical value.

To explore experimentally the consequences of mitochondrial Ca²⁺ uptake on global Ca²⁺ signals, mitochondrial Ca²⁺ uptake is often prevented directly by applying an MCU inhibitor¹⁴ or modulated by using agents which either depolarize or hyperpolarize the mitochondrial membrane potential and thus dissipate or increase the driving force for Ca²⁺ uptake into the organelle^10,14. How mitochondrial Ca²⁺ uptake affects IP₃R release mechanism is still under investigation. Some researchers have argued that mitochondrial Ca²⁺ uptake may provide a positive modulation of Ca²⁺ release by suppressing the negative feedback of high [Ca²⁺]_Cyt on IP₃R^9,10,11. Others have suggested that mitochondrial Ca²⁺ uptake exerts a negative control on [Ca²⁺]_Cyt by preventing the positive feedback of Ca²⁺ on IP₃R^12,13.

In our simulations, the ability of mitochondria to take up Ca²⁺ is changed by varying the IP₃R-MCU distance and so the mitochondrial modulation on IP₃R dynamics has been examined. Our model indicates that altering mitochondrial Ca²⁺ uptake induces a nontrivial change on [Ca²⁺]_Cyt amplitude, which can reconcile the conflicting results of the change of [Ca²⁺]_Cyt amplitude observed in various experiments^{9,10,11,12,13}. When IP₃R-MCU distance is larger than the critical distance r_crt, mitochondrial Ca²⁺ uptake is weak and more Ca²⁺ ions diffuse into cytosol, the enhancement of mitochondrial Ca²⁺ uptake causes a decreasing [Ca²⁺]_Cyt spike peak, exerting a weakened inhibition effect on IP₃R dynamics, resulting in a prolonged decay time of spikes. Thus, at this region, mitochondrial Ca²⁺ uptake offers a positive modulation of Ca²⁺ release by reducing [Ca²⁺]_Cyt amplitude to suppress the inhibition dynamics (i.e. the negative feedback) of high [Ca²⁺]_Cyt on IP₃R during decay time of spikes (Fig. 7), similar as suggested in refs 9, 10, 11. While, when IP₃R-MCU distance is smaller than the critical distance, mitochondria take more Ca²⁺ ions than those diffusing into cytosol, the strong mitochondrial Ca²⁺ uptake will force [Ca²⁺]_Cyt at a small value to prevent the strong activation dynamics (i.e. the positive feedback) of [Ca²⁺]_Cyt on IP₃R (Fig. 7), similar as proposed in refs 12, 13. Thus, our model shows that this weak activation effect causes a small Ca²⁺ release from ER, giving a slow rising phase. When mitochondria approach to saturation and [Ca²⁺]_Cyt increases beyond the threshold of excitable dynamics, an abrupt rising phase occurs to generate a [Ca²⁺]_Cyt spike, giving a high [Ca²⁺]_Cyt amplitude and then inducing a strong inhibitory effect on IP₃R dynamics to shorten the decay time of spikes.

Figure 7

Schematic representation of how mitochondria modulate [Ca²⁺]_Cyt.

IP₃R channel is regulated by Ca²⁺ in a biphasic manner: a lower [Ca²⁺]_Cyt promotes Ca²⁺ release and causes its elevation whereas higher [Ca²⁺]_Cyt inhibits Ca²⁺ release. There is a critical distance at which 50% of the IP₃R-released Ca²⁺ ions are taken up by mitochondria. When r<r_crt, mitochondrial Ca²⁺ uptake exerts a negative control on IP₃R-released Ca²⁺ by preventing the positive effect of low [Ca²⁺]_Cyt on IP₃R; whilst when r > r_crt, mitochondrial Ca²⁺ uptake provides a positive modulation on IP₃R-released Ca²⁺ by suppressing the negative effect of high [Ca²⁺]_Cyt on IP₃R. Ca²⁺ oscillations merely occur for moderate values of [IP₃] (orange range) with physiological amplitude (gray range), so the optimal distance of IP₃R-MCU should be within the intersection zone between the gray and orange range.

An unexpected prediction with the model is that when the IP₃R-MCU distance is less than the critical distance, even though mitochondrial Ca²⁺ uptake keeps increasing with decreasing distance, the total amount of Ca²⁺ released from the ER decreases while the cytosolic Ca²⁺ spike amplitude increases. The [Ca²⁺]_Cyt peak is mainly determined by the amount of Ca²⁺ released during the abrupt rising phase of spikes. At shorter IP₃R-MCU distances, the mitochondria become almost saturated before the abrupt rising phase and thus during the abrupt rising phase more of the Ca²⁺ ions released from the ER diffuse more into cytosol, generating a higher amplitude of [Ca²⁺]_Cyt. However, the change in the total amount of Ca²⁺ release from the ER seen with decreasing distance is mainly determined by the release of Ca²⁺ occurring during the decay phase of spike. The enhanced inhibitory effect of higher amplitude [Ca²⁺]_Cyt spikes on IP₃R dynamics in turn leads to a shortening of the decay time of the spikes, resulting in a decrease in the total Ca²⁺ amount released from ER with decreased IP₃R-MCU distance.

The existence of high [Ca²⁺]_Mic in microdomains has been put forward as one way to resolve the apparent paradox that even though the Ca²⁺ affinity of the MCU is low (K_d of 10–20 μM) highly efficient mitochondrial Ca²⁺ uptake still occurs under physiological conditions with low micromolar global [Ca²⁺]_Cyt³. However, an intriguing question related to the role of the microdomain is the distance between the ER and mitochondria that would be required to control distinct cellular processes. The distance from the ER membrane to OMM was originally estimated to be in excess of 100 nm¹⁶, recently it has been suggested that this distance could be in the order of 30 nm^15,55. A key prediction of our model is that the optimal distance between the IP₃R and MCU is about 30–85 nm. This arises when parameters for [IP₃] and [Ca²⁺] are derived from the experimental observations that Ca²⁺ oscillations generally occur for intermediate values of IP₃⁵⁰ and that physiological [Ca²⁺]_Cyt amplitudes seldom exceed 3.5 μM in living cells^38,48, which are illustrated in Fig. 7. Considering that the thickness between OMM and IMM is around 20 nm⁵⁶, our model predicts that the optimal gap between the ER membrane and OMM should be approximately 10–65 nm, which is in good agreement with measurements obtained via electron microscopy and tomography^15,31,55.

Interestingly, our simulation also provides a reason for why MCUs need to be located ~20 nm away from the OMM. In experiment, elevations in the magnitude of [Ca²⁺]_Cyt during typical agonist induced responses could lead to a triggering of [Ca²⁺]_Cyt-induced apoptosis⁵¹. Maintaining a distance between IMM and OMM of 20 nm may present a defense mechanism to keep the global [Ca²⁺]_Cyt at a small physiological range. Our model highlights the importance of the IP₃R-MCU distance and suggests that potential pathological mechanisms for generating elevated [Ca²⁺]_Cyt and leading to possibly [Ca²⁺]_Cyt-induced apoptosis^6,51 would be found in disorders in which the optimal distance between the ER and mitochondria is disturbed. Considering that when untethered to the ER, mitochondria can be highly mobile organelles⁵⁴ and that little is known of the physiological stimuli governing the distance between the ER and mitochondria, or how the tethers between the ER and mitochondria function or are disrupted during pathology, future research elucidating how the IP₃R-MCU distance is maintained and controlled is essential.

Finally, we discuss the advantages and limitations of our model. First, the previous published models typically simulated MCU either by a simple Hill equation^26,29 or by the Monod-Wyman-Changeux type allosteric regulation^24,27,28. In contrast to these models, we introduced a mechanistic kinetic model to discuss in detail the binding and unbinding processes of Ca²⁺ to MCU. Second, the numerical simulation about the peak values of [Ca²⁺]_ER or/and [Ca²⁺]_Mt were usually underestimated in most of the models^{24,26,27,28,29,30}. For example, the concentration of free Ca²⁺ in the ER is typically in the range of a few hundred micromolars. Our model is capable of reproducing the observed magnitudes of [Ca²⁺]_Cyt, [Ca²⁺]_ER and [Ca²⁺]_Mt.

In 2013, Szopa et al. suggested a model to explicitly take into account the existence of Ca²⁺ microdomain between the ER and mitochondria to discuss the mitochondrial Ca²⁺ uptake²⁹. However, the assumption in the model that MCU directly senses the [Ca²⁺]_ER may be not very appropriate, particularly in view of the existing sharp Ca²⁺ gradients in the microdomain against the distance between the IP₃R and MCU. Previously, there are mainly two basic ways to handle Ca²⁺ microdomain in the whole-cell models. The first approach involved the separation of a specified small region from the cell. But this subspace is simply treated as another well-mixed compartment^57,58 and thus this approach cannot capture the fact of sharp Ca²⁺ concentration gradients in the microdomain^43,44. Sneyd et al. used an alternative approach to divide the cell into a few separate domains, connected by many boundary fluxes^59,60. Although this is a powerful approach, it greatly increases the complexity of the model, leading to extreme computational difficulties.

Different from these methods, a new approach has been suggested in the present paper. For simplicity, all the IP₃Rs dynamics and SERCA are determined by the bulk [Ca²⁺]_Cyt. As a result, the bulk Ca²⁺ concentrations in cytosol, ER and mitochondria can be examined in deterministic ODE which actually assuming that the cytosolic, ER and mitochondrial compartments are mixed homogenously. Upon such well mixed compartments, a microdomain is specifically considered to discuss the MCU dynamics which is crucial for mitochondrial Ca²⁺ uptake. The microdomain consists of a cluster of IP₃Rs and a MCU separated at a short distance. Treated as a point release source for the clustered IP₃Rs which dynamics is controlled by bulk [Ca²⁺]_Cyt, a sharp distribution of [Ca²⁺]_Mic can be built up in the microdomain. Each MCU responds to the local high [Ca²⁺]_Mic generated by the clustered IP₃Rs, rather than the bulk [Ca²⁺]_Cyt. As a result, this method avoids the simulation of the Ca²⁺ diffusion in microdomain between the ER membrane and IMM.

Peglow et al.³⁰ discussed a three-dimensional model to analyze the influence of the distance between mitochondria and the CRAC channel on the plasma membrane, which is responsible for Ca²⁺ influx from the outside of the cell, on global Ca²⁺ oscillation. They confirmed that the relative position between mitochondria and the CRAC channels controls the global Ca²⁺ signal. Their simulation results showed that close mitochondria will lower the microdomain [Ca²⁺] right at the CRAC channels and prevent the Ca²⁺-dependent inactivation of CRAC channels, leading to higher Ca²⁺ influx and thereby increasing global [Ca²⁺]_Cyt. They found that the global [Ca²⁺]_Cyt monotonously changes with the mitochondria-CRAC distance from 0.1 to 2 μm. A similar result is obtained in our model at r > r_crt.

As a fact, our model does not realistically consider the spatial structure and so many important factors of Ca²⁺ signaling have been ignored. In the model the microdomain is modeled as a radially symmetric hemisphere and so the linearized reaction-diffusion equation has been applied to obtain a simple expression for [Ca²⁺]_Mic. The mitochondrial boundary is actually ignored. One may expect that the mitochondrial boundary effect may modulate the sharp [Ca²⁺]_Mic distribution. The Ca²⁺ release from the clustered IP₃Rs are treated as the point release source, which is correct for far away MCU. But if MCU is too close to IP₃Rs, a better treatment is to consider the spatial distribution of IP₃Rs. The model also assumes that OMM has no effect on the Ca²⁺ diffusion as it has high permeability to Ca²⁺. A better MCU model should also consider the regulation effect of MICU2 on the MCU channel activity. Our model is primarily designed for non-excitable cells, i.e. a ‘closed’ cell model, ignoring the Ca²⁺ exchange between cytosol and extracellular medium through the plasma membrane. Thus, a more realistic Ca²⁺ model should be considered in the future. However, despite oversimplification in some details, we believe that the present model is useful and may provide us some insights on how the crosstalk between the ER and mitochondria modulates Ca²⁺ signals in living cells.