The Function of Mitochondrial Calcium Uniporter at the Whole-Cell and Single Mitochondrion Levels in WT, MICU1 KO, and MICU2 KO Cells

Mitochondrial Ca2+ ([Ca2+]M) uptake through its Ca2+ uniporter (MCU) is central to many cell functions such as bioenergetics, spatiotemporal organization of Ca2+ signals, and apoptosis. MCU activity is regulated by several intrinsic proteins including MICU1, MICU2, and EMRE. While significant details about the role of MICU1, MICU2, and EMRE in MCU function have emerged recently, a key challenge for the future experiments is to investigate how these regulatory proteins modulate mitochondrial Ca2+ influx through MCU in intact cells under pathophysiological conditions. This is further complicated by the fact that several variables affecting MCU function change dynamically as cell functions. To overcome this void, we develop a data-driven model that closely replicates the behavior of MCU under a wide range of cytosolic Ca2+ ([Ca2+]C), [Ca2+]M, and mitochondrial membrane potential values in WT, MICU1 knockout (KO), and MICU2 KO cells at the single mitochondrion and whole-cell levels. The model is extended to investigate how MICU1 or MICU2 KO affect mitochondrial function. Moreover, we show how Ca2+ buffering proteins, the separation between mitochondrion and Ca2+-releasing stores, and the duration of opening of Ca2+-releasing channels affect mitochondrial function under different conditions. Finally, we demonstrate an easy extension of the model to single channel function of MCU.

Significant information about the role of regulatory proteins described above in the function of MCU has emerged over the last few years [10][11][12][13][14][15][18][19][20][21][22]. In particular, MICU1 has been shown to mediate the observed suppression of MCU activity in the low cytosolic Ca 2+ concentration ([Ca 2+ ] C ) regime, called "gatekeeping" [12]. This gatekeeping of MCU activity is relieved when [Ca 2+ ] C reaches approximately 1.3 µM. Below this critical concentration, Ca 2+ uptake through MCU is essentially negligible. MICU1 is also responsible for the observed highly cooperative activation of MCU activity where the gatekeeping is cooperatively relieved when [Ca 2+ ]C exceeds 1.3 μM [12] ( Figure 1B). Similarly, MICU2 interacts with MICU1 to increase the [Ca 2+ ]C threshold for the gatekeeping and reduce the gain of cooperative activation of MCU activity [12] ( Figure 1B). EMRE, on the other hand, plays a role in regulating MCU activity on the matrix side of the inner mitochondrial membrane (IMM) such that the uptake is minimum when [Ca 2+ ]M is approximately 400 nM [18,19]. Thus, MCU exhibits a biphasic behavior where the current through the channel increases as [Ca 2+ ]M drops below or rises above 400 nM ( Figure 1C). Furthermore, MICU1 and MICU2 are required for the matrix Ca 2+ regulation of MCU such that the inhibition at the intermediate [Ca 2+ ]M disappears when either one of the two proteins are knocked out [19].  [12] (B) and [19] (C-E).
Despite this wealth of information, many key issues about the activity of MCU and the role of different regulatory proteins in mitochondrial function remain unresolved and are beyond the scope of current experimental techniques. For instance, while experimental tools can be used to study the activity of MCU or how it is affected by a given regulatory protein at discrete fixed values of one or two variables, several variables affecting mitochondrial Ca 2+ uptake in real cells vary continuously in mutually dependent manner. Furthermore, the observations about MCU activity are based on experiments performed at different spatiotemporal scales ranging from patch clamp electrophysiology of single Ca 2+ uniporter channels and individual mitoplasts to the imaging of Ca 2+ uptake at the cell culture level using fluorescence microscopy (e.g., see [6,11,12,19]). Linking these  [12] (B) and [19] (C-E).
Despite this wealth of information, many key issues about the activity of MCU and the role of different regulatory proteins in mitochondrial function remain unresolved and are beyond the scope of current experimental techniques. For instance, while experimental tools can be used to study the activity of MCU or how it is affected by a given regulatory protein at discrete fixed values of one or two variables, several variables affecting mitochondrial Ca 2+ uptake in real cells vary continuously in mutually dependent manner. Furthermore, the observations about MCU activity are based on experiments performed at different spatiotemporal scales ranging from patch clamp electrophysiology of single Ca 2+ uniporter channels and individual mitoplasts to the imaging of Ca 2+ uptake at the cell culture level using fluorescence microscopy (e.g., see [6,11,12,19]). Linking these observations at different scales is key to understanding the role of Ca 2+ in mitochondrial and cell function. For example, while significant data about mitochondria-dependent variables such as ATP and reactive oxygen species as a function of Ca 2+ at the whole-cell or cell culture levels exists, mitochondrial Ca 2+ uptake and function are mainly regulated by local Ca 2+ concentration in a few tens of nanometers wide microdomain formed by the close apposition of plasma membrane or intracellular organelles with individual mitochondrion [23][24][25][26]. These and a range of other observations underscore the importance of incorporating the key observations about the role of regulatory proteins in the function of MCU in a comprehensive computational framework.
The importance of a computational model incorporating the role of regulatory proteins in the function of MCU is further highlighted by the role of these proteins in various pathologies. For example, the human loss-of-function mutations in MICU1 results in mitochondrial Ca 2+ overload, impaired bioenergetics, and mitochondrial fragmentation [27]. These MICU1-induced defects in mitochondrial Ca 2+ signaling and function are believed to result in early-onset neuromuscular weakness, impaired cognition, and extrapyramidal motor disorder [27][28][29]. A whole-body knockout of MICU1 has also been reported to result in perinatal lethality in mouse models [30,31]. In addition, the expression levels of EMRE were also altered in patients with MICU1 mutations compared to controls [27]. Recently, it has been shown that mutations in mitochondrial m-AAA proteases associated with spinocerebellar ataxia and hereditary spastic paraplegia inflict their cytotoxicity by affecting the biogenesis of EMRE. This leads to the accumulation of active MCU-EMRE channels lacking gatekeeping, which facilitates mitochondrial Ca 2+ overload, opening of mitochondrial permeability transition pore, and neurodegeneration [32].
In this paper, we (1) develop a comprehensive data-driven model that reproduces the key observations about the function of MCU under a range of [Ca 2+ ] C , [Ca 2+ ] M , and mitochondrial membrane potential (MMP or ∆ψ) values in WT, MICU1 KO, and MICU2 KO cells; (2) combine the observations at the cell culture (whole-cell) and single mitoplasts levels into a single mathematical framework; (3) investigate how the Ca 2+ signaling in the microdomain between the endoplasmic reticulum (ER) and mitochondrion compares in WT, MICU1 KO, and MICU2 KO mitoplasts by modeling the interaction between the ER-bound inositol 1,4,5-trisphosphate (IP 3 ) receptor (IP 3 R) channel and MCU; (4) how does MICU1 or MICU2 KO affect mitochondrial function both in terms of Ca 2+ uptake and ATP production; and (5) investigate the effect of Ca 2+ buffers, width of the microdomain, and open duration of IP 3 R on mitochondrial Ca 2+ uptake and ATP production in WT, MICU1 KO, and MICU2 KO cells.

Materials and Methods
The experimental data used in this paper is adapted from Riley et al. [12] and Vais et al. [19] with permission. These experiments were performed on HEK293 cells. We refer the interested reader to these two papers for details about experimental methods.

Kinetic Model for MCU Function
Our kinetic model for MCU implements the scheme proposed in Payne et al. [12] and Vais et al. [19].  (50-800 nM), Ca 2+ binds to the inhibitory sensor that closes the channel (states C 31 and C 41 in the model shown in Figure 1A). At low [Ca 2+ ] M , the matrix Ca 2+ -mediated gatekeeping is released (state O 40 ). At high [Ca 2+ ] M , Ca 2+ binds to the activating senor, relieving the gatekeeping on the matrix side of the IMM. Our best fit criterion (BIC score) warrants the binding of three Ca 2+ to the activating sensor in addition to the one Ca 2+ that is bound in state C 41 to open the channel again (state O 44 ). More recently, it was shown that the affinities of activating and inhibitory sensors on the matrix side are regulated by another Ca 2+ sensor termed "flux sensor" [18]. The flux sensor is believed to represent a site on the matrix side of IMM that binds the Ca 2+ as it enters through the channel. This sensor is sensitive to the Ca 2+ buffering capacity of mitochondria, especially fast buffers that uptake Ca 2+ immediately after entering the matrix. Although our study does not concern the role of mitochondrial Ca 2+ buffers, our model can reproduce the observations about the role of "flux sensor" in mitochondrial Ca 2+ uptake without changing the topology of the model (results not shown). We also remark that EMRE plays a key role in the matrix Ca 2+ -mediated regulation of MCU activity, and expressing mutant EMRE abolished this regulation. EMRE is also key to the strong coupling between the mechanisms regulating the channel activity on the matrix and the opposite side of IMM.
Note that, in our model, there are transitions that involve binding of more than one Ca 2+ . For example, the MCU has 0 and 2 Ca 2+ bound to the sensors on the cytosolic side of IMM in states C 00 and C M1 20 , respectively. However, the direct link between C 00 and C M1 20 does not necessarily mean that the channel binds 2 Ca 2+ simultaneously. It rather means that the state with 1 Ca 2+ bound have very low occupancy and is not required in the model. However, these transition states act as speed-bumps for the probability flux, and their effect can be incorporated in the mean transition times or transition rates. Similarly, for the connectivity of the model, we apply Occam's razor and include the transitions that are warranted by Bayesian Information Criterion [33]. These issues have been discussed in detail elsewhere [34].
The above considerations lead to the kinetic scheme with six close (C XY ) and two open (O XY ) states ( Figure 1A), where subscripts X and Y represent the number of Ca 2+ bound to the domains on the cytosolic and matrix sides of the IMM, respectively. Relative to the reference unliganded close state C 00 , the occupancies (unnormalized probabilities) of the close C XY and open O XY states are proportional to ([Ca 2+ ] C ) X ([Ca 2+ ] M ) Y with occupancy parameters K C XY and K O XY , respectively. The occupancy parameter of a state is the product of equilibrium association constants along any path connecting C 00 to that state (K C 00 = 1) [34,35]. In other words, the product of all forward rates starting from C 00 to a given state X divided by the product of the rates from state X back to C 00 gives the occupancy of state X. This applies to all states and is independent of the path taken from C 00 to X. For example, the ratio of the transition rate from C 00 to C M1 20 to the rate from C M1 20 to C 00 gives the occupancy of C M1 20 . This will become clear further when we derive the transition rates for the single MCU channel in the Discussion section. We refer the interested reader to [34,35]  /Z, respectively, where Z is the total occupancy of all states and is given as where Z C and Z O are the unnormalized occupancies of all close and all open states, respectively. The P O of the channel can be written as The Ca 2+ uptake rate of MCU (V MCU ) as a function of [Ca 2+ ] C and [Ca 2+ ] M is given by where V MCU,max = 600 nM/s is the maximum uptake rate observed in the cell culture experiments [12]. We remark that the initial Ca 2+ uptake rate for each experiment was determined by fitting an exponential function to [Ca 2+ ] C beginning from Ca 2+ addition until a new steady-state is reached 300 s later to obtain parameters A (extent of uptake) and τ (time constant). The instantaneous rate of uptake at t = 0 was taken as equal to A/τ. In these experiments, only mitochondrial Ca 2+ uptake through MCU was active, whereas all other Ca 2+ pathways were either pharmacologically blocked or were inactive (see Ref. [12] for details).

Whole-Cell Ca 2+ Model
In the cell culture experiments, HEK293 cells were treated with 0.004% digitonin at 50 s after the beginning of the experiment to permeabilize plasma membrane in intracellular-like medium lacking free Ca 2+ buffers, equilibrating Ca 2+ concentration in bath solution and cytoplasm. Ca 2+ uptake into the ER and mitochondrial efflux were inhibited by applying 2 µM thapsigargin (Sarco/ER Ca 2+ ATPase (SERCA) blocker) and 20 µM CGP37157 (Na + /Ca 2+ exchange blocker) at 100 and 400 s, respectively. Similarly, no inositol 1,4,5-trisphosphate (IP 3 ) was used in these experiments, implying no Ca 2+ efflux from the ER through IP 3 receptor (IP 3 R) channels. After [Ca 2+ ] C reached a steady-state at 700 s, MCU-mediated Ca 2+ uptake was initiated by adding boluses of CaCl 2 to achieve increases in [Ca 2+ ] C between 100 nM and 10 µM. All these considerations lead to the following rate equations for [Ca 2+ ] C and [Ca 2+ ] M , respectively.
d Ca 2+ where V Stim represents Ca 2+ stimulus (bolus of CaCl 2 ). The value of V Stim is set so that the peak [Ca 2+ ] C given by the model matches the observed value. Parameters f m = 3 × 10 −4 and δ = 0.067 are the Ca 2+ buffering of mitochondria and the ratio of mitochondrial to cytosolic volume [36,37]. As we will see later, both of these parameters significantly affect free [Ca 2+ ] M . We add V Leak to the rate equation for [Ca 2+ ] C to incorporate the effect of potential Ca 2+ leak from the ER to the cytoplasm: where Ca 2+ ER = 200 µM is Ca 2+ concentration in the ER.

The Effect of MICU1 or MICU2 KO on Mitochondrial Function at the Whole-Cell Level
To assess how MICU1 or MICU2 KO affect mitochondrial function, we extend the model to incorporate the dynamics of mitochondrial NADH concentration ([NADH] M ), mitochondrial ADP concentration ([ADP] M ), ∆ψ, and cytosolic ADP concentration ([ADP] C ). The rate equations, relevant fluxes (with the exception of V MCU and V NaCa ), and parameters are adopted from Ref. [36] and are given in the Supplementary Information Text. Furthermore, Equation (5) is modified slightly to include Ca 2+ efflux from the mitochondria through Na + /Ca 2+ exchanger (V NaCa ) [38].
[Na] i = 11 × 10 3 µM, F, R, and T = 310.16 o K is the cytosolic Na + concentration, Faraday's constant, gas constant, and temperature, respectively. Note that ∆ψ 0 = 91 mV does not represent the resting MMP. It was originally selected so that V NaCa is maximum when ∆ψ = ∆ψ 0 , i.e., under uncoupled conditions [38]. As we discuss later, V NaCa,max is estimated from the observed values of resting [Ca 2+ ] M in WT, MICU1, and MICU2 KO conditions. The full model used for mitochondrial function at the whole-cell level consists of Equations (7) (with Equations (3) and (8)) and Equations (S1)-(S4) (Section "Bioenergetics model" in the Supplementary Information Text). Equation (3) is multiplied by Equation (9) given below to incorporate the voltage dependence of V MCU .

The Model Reproduces the Function of MCU both at Cell Culture and Single Mitoplast Levels
Model fits to the mitochondrial Ca 2+ uptake rate as a function of [Ca 2+ ] C observed in cell cultures using fluorescence microscopy on WT (spheres, solid line), MICU1 KO (squares, dashed line), and MICU2 KO (diamonds, dashed-dotted line) cells are shown in Figure 1B. The values of occupancy parameters for the three conditions are given in Table S2. The model closely reproduces the general features of the MCU activity including the suppressed activity in the low [Ca 2+ ] C regime (0-1.3 µM), the cooperative opening once the gatekeeping is relieved, relief of gatekeeping and reduced cooperative activation in MICU1 KO cells, and reduced threshold for gatekeeping and increased gain of cooperative activation in MICU2 KO cells. We remark that our fitting criterion searches for the parameters that result in the best fit to all observations simultaneously. For example, in case of MICU1 KO experiments, we do not restrict our search to states that are associated with MICU1 only. In other words, we make no a priori assumption about the independent binding of Ca 2+ to the MICU1 EF hands. This is motivated by the complex structure of MCU [5,8,9,39] and the understanding that "the protein is not a rigid system in which a ligand moves in a fixed potential. Instead, there is a strong mutual interaction between ligand and protein" [40] that affects several aspects of the channel gating. This approach leads to an interesting observation in our model: MICU1 (or MICU2) KO not only affects the probability of states associated with MICU1 but also change the occupancy of other binding sites (  Figure 3J,K in Ref. [19] Figure 1B go beyond 8 µM. We believe that this information gap could be causing the discrepancy between the model and observations in the MICU2 KO cells ( Figure 1B). We noticed that another set of parameters give a closer fit to the uptake rate observed in the MICU2 KO cell culture, but it deteriorates the fit to the MCU current (I MCU ) observed at the single mitoplast level [19] ( Figure S1). Since both [Ca 2+ ] C and [Ca 2+ ] M are known in the patch clamp experiments, we choose the first set of parameters (Table S2) over the second for the rest of this paper.
Payne et al. [12] suggested that the state with Ca 2+ bound to both EF hands of MICU1 and one EF hand of MICU2 is an open state. To test this hypothesis, we change the state C 30 in the model to an open state (O 30 ). However, we found that such change is not supported as the model predicts higher Ca 2+ uptake rate in WT cells ( Figure S2A). Furthermore, the overall fit to the observed uptake rate in MICU1 KO cells significantly deteriorates ( Figure S2B). We therefore leave this state as a close state from now on.
In our formalism, the transition from whole-cell model (based on cell culture experiments) to the single mitoplast model is straightforward. The only parameter that changes while fitting to the single mitoplast data is that V MCU,max in Equation (3) Figure 1D). The discrepancy in case of MICU1 KO could be due to the fact that the model represents MICU1 KO mitoplasts, whereas the experimental data are from MICU1 knockdown mitoplasts and is shown only for a crude comparison due to the lack of such data on MICU1 KO mitoplasts. Indeed, the overall theoretical trend is consistent with the observations in Figure 1C (Figure 2A). In other words, the range of [Ca 2+ ] M values where gatekeeping is effective gets narrower as [Ca 2+ ] C increases. This, together with the cooperative opening of MCU once the gatekeeping on the cytosolic side is relieved, could potentially lead to a vicious cycle, leaving mitochondria increasingly vulnerable to Ca 2+ overload and the opening of mitochondrial permeability transition pore [41] as [Ca 2+ ] C increases beyond a certain limit.

Membrane Potential Dependence of MCU Current Density
Membrane potential (Δψ) is another key variable influencing current through MCU by providing a driving force for Ca 2+ influx. Figure 2C shows that the inward current density through MCU decreases as we make mitochondrial membrane more depolarized. While the amplitude of IMCU/Cm changes as we change [Ca 2+ ]C, [Ca 2+ ]M, and knockout MICU1 or MICU2, normalizing all traces with respect to their respective peak values shows that they follow the same functional form with respect to Δψ. Thus, we use the following function to represent the Δψ dependence of Ca 2+ flux through MCU: where A = 1.62349 and b = 2.69 × 10 -4 . Fit to IMCU/Cm given by the model is shown by the thick solid line in Figure 2D. Multiplying Eqs. (3 and 9) leads to a complete model for MCU Ca 2+ uptake rate as a function of [Ca 2+ ]C, [Ca 2+ ]M, and Δψ. In Figure 2E, [44,45]. Such concentrations would cause the gatekeeping on both sides of the IMM to relieve, resulting in a maximum flux through MCU [12,19]. Furthermore, normalized TMRE fluorescence representing membrane potential (Δψ) increases by more than 0.1 on a 0.5 to 1 scale (where 1 represents the complete dissipation of the potential difference across IMM) when [Ca 2+ ]C is

Membrane Potential Dependence of MCU Current Density
Membrane potential (∆ψ) is another key variable influencing current through MCU by providing a driving force for Ca 2+ influx. Figure 2C shows that the inward current density through MCU decreases as we make mitochondrial membrane more depolarized. While the amplitude of I MCU /C m changes as we change [Ca 2+ ] C , [Ca 2+ ] M , and knockout MICU1 or MICU2, normalizing all traces with respect to their respective peak values shows that they follow the same functional form with respect to ∆ψ. Thus, we use the following function to represent the ∆ψ dependence of Ca 2+ flux through MCU: where A = 1.62349 and b = 2.69 × 10 −4 . Fit to I MCU /C m given by the model is shown by the thick solid line in Figure 2D. Multiplying Equations (3) and (9) [44,45]. Such concentrations would cause the gatekeeping on both sides of the IMM to relieve, resulting in a maximum flux through MCU [12,19]. Furthermore, normalized TMRE fluorescence representing membrane potential (∆ψ) increases by more than 0.1 on a 0.5 to 1 scale (where 1 represents the complete dissipation of the potential difference across IMM) when [Ca 2+ ] C is raised by 5 µM in patch clamp experiments on WT mitoplasts (see Figure S3 in Ref. [12]). The drop in ∆ψ would cause a decrease in Ca 2+ uptake by mitochondria, giving another crucial control to these key variables to regulate mitochondrial function. We remark that all analysis in Figure 2 is applicable to MCU function at the cell culture level except that the peak I MCU /C m (206 pA/pF) gets replaced by the maximum uptake rate observed at the whole-cell level (600 nM/s).  Figure 3B).
The peak [Ca 2+ ] M (peak free Ca 2+ in the mitochondria) given by the model is an order of magnitude lower than that observed experimentally ( Figure 3C,D). For example, raising [Ca 2+ ] C close to 800 nM in WT cells increases [Ca 2+ ] M to more than 1 µM ( Figure 3J,K in Ref. [12]). However, the Ca 2+ buffering capacity of mitochondria (f m ) and the relative volume of mitochondria with respect to cytosol (δ) both play key roles in the free [Ca 2+ ] M . We repeat the simulation in Figure 3A for WT cells using Ca 2+ bolus of different strengths and record the peak [Ca 2+ ] C during the simulation and [Ca 2+ ] M 300 s after the Ca 2+ bolus is added to the cytoplasm. As shown in Figure 3E,F, the peak free mitochondrial Ca 2+ increases dramatically as we increase f m (decreasing mitochondrial Ca 2+ buffering capacity) or decrease δ (decreasing the relative mitochondrial volume). Using the simultaneously measured values of [Ca 2+ ] C and [Ca 2+ ] M from Ref. [12] as a reference, we found that f m~0 .09 results in a close correspondence between [Ca 2+ ] C and [Ca 2+ ] M when δ = 0.067 is considered. The f m value giving the correct correspondence between model and observations would decrease (higher mitochondrial Ca 2+ buffering capacity) as we decrease δ.

Mitochondrial Bioenergetics in WT, MICU1 KO, and MICU2 KO Cells
In simulations shown in Figure 3, Δψ is kept fixed at -160 mV. However, Δψ changes as mitochondria buffers Ca 2+ . We therefore extend our model for mitochondrial Ca 2+ uptake at the whole-cell level to incorporate the dynamic behavior of Δψ and how mitochondrial NADH ([NADH]M) and ATP ([ATP]M) change as mitochondria buffer Ca 2+ in WT, MICU1 KO, and MICU2 KO cells. A key parameter missing in our whole-cell Ca 2+ model is the maximum flux through Na + /Ca 2+ exchanger (VNaCa,max in Eq. 8). Although we do not have direct access to this parameter in our experimental results, we use the observed resting value of [Ca 2+ ]M to estimate VNaCa,max. Knowing that the Ca 2+ fluxes through MCU and Na + /Ca 2+ exchanger balance each other out in resting conditions [27], and the MCU flux is constrained by the observations discussed above, we vary VNaCa,max until

Mitochondrial Bioenergetics in WT, MICU1 KO, and MICU2 KO Cells
In simulations shown in Figure 3, ∆ψ is kept fixed at −160 mV. However, ∆ψ changes as mitochondria buffers Ca 2+ . We therefore extend our model for mitochondrial Ca 2+ uptake at the whole-cell level to incorporate the dynamic behavior of ∆ψ and how mitochondrial NADH ([NADH] M ) and ATP ([ATP] M ) change as mitochondria buffer Ca 2+ in WT, MICU1 KO, and MICU2 KO cells. A key parameter missing in our whole-cell Ca 2+ model is the maximum flux through Na + /Ca 2+ exchanger (V NaCa,max in Equaiton (8)). Although we do not have direct access to this parameter in our experimental results, we use the observed resting value of [Ca 2+ ] M to estimate V NaCa,max . Knowing that the Ca 2+ fluxes through MCU and Na + /Ca 2+ exchanger balance each other out in resting conditions [27], and the MCU flux is constrained by the observations discussed above, we vary V NaCa,max until [Ca 2+ ] M given by the model in resting state matches the observed value ( Figure S3). The values given by the fit in WT, MICU1 KO, and MICU2 KO conditions at the whole-cell level are given in the legends of Figure S3. In line with the observations where the deregulated mitochondrial Ca 2+ uptake due to the loss of MICU1 function is shown to initiate a futile cycle whereby continuous Ca 2+ influx is balanced by efflux through Na + /Ca 2+ exchange [27], our model predicts higher V NaCa,max in MICU1 KO cells (2.144 nM/s) when compared to WT cells (0.958 nM/s).
After reaching steady state, we apply a square pulse of Ca 2+ to the cytoplasm, raising [Ca 2+ ] C from its resting value of 0.1 µM to 1.  Figure S3C) that also leads to higher [ATP] M in resting conditions ( Figure S3D). The small depolarization in IMM as soon as [Ca 2+ ] C is increased (note the very small increase in ∆ψ before the relatively large decrease) results from the Ca 2+ and Na + influx through MCU and Na + /Ca 2+ exchanger, respectively. The small depolarization is followed by a relatively large hyperpolarization (although still very small in absolute terms), mainly due to the increase in [ [27] showing that the futile Ca 2+ cycle due to uninhibited MCU influx balanced by larger efflux through Na + /Ca 2+ exchanger leads to a lower ATP production in cells with MICU1 mutations.

Ca 2+ Uptake and Bioenergetics at the Single Mitoplast Level
There are two main discrepancies between ∆ψ given by the model at the whole-cell level and that observed experimentally in the cell culture experiments (see Figure S3 in Ref. [12]). First, the observed IMM potential in MICU1 KO cells is more depolarized than the WT and MICU2 KO cells. Second, the observed depolarization of IMM after the application of Ca 2+ bolus is significantly larger when compared to the model predictions where the depolarization due to Ca 2+ influx is almost unnoticeable. This is mainly due to the fact that the estimated maximum flux through MCU at the whole-cell level is small as it represents the average uptake of the entire mitochondrial network in the cell culture. That is, it represents the rate at which the average extra-mitochondrial Ca 2+ in the cell culture decreases due to the net uptake by the entire mitochondrial network. Similarly, [Ca 2+ ] C and [Ca 2+ ] M also represent the average concentrations and ignore the large concentration gradients that could occur in the microdomain formed by the close apposition of a mitochondrion and the ER. In reality, a mitochondrion experiences Ca 2+ on the cytosolic side that could be significantly higher than the bulk [Ca 2+ ] C . Furthermore, as shown below, the maximum flux through MCU estimated from the patch-clamp experiments on single mitoplasts comes out to be significantly larger. We would like to point out that the experiments in [19] investigated the function of MCU at the single mitoplast level, not the bioenergetics of the mitochondrion. Thus, we model bioenergetics at the single mitochondrion level as if it is residing inside a living cell, but the MCU functions the same way as it would in the mitoplast.
The maximum value of MCU current density (I MCU /C m ) observed in the patch-clamp experiments is equal to 206, 237, and 217 pA/pF in the case of WT, MICU1 KO, and MICU2 KO mitoplasts respectively at ∆ψ = −160 mV [19]. These values can be used to estimate the maximum flux through MCU (V MCU,max in Equation (3)) as follows: where Vol m is the volume of the mitochondrion. We consider a spherical mitochondrion of radius 0.5 µm. Vais et al. [19] observed the mitoplast capacitance (C m ) to be in the range 0. In intact cells, mitochondrial Ca 2+ uptake at the single mitochondrion level also depends on the spatial distance between the ER and mitochondrion. Ca 2+ is released from the ER through IP 3 R that diffuses away from the channel in the microdomain ( Figure S4). We consider a single IP 3 R, allowing it to open for 100 ms, and model Ca 2+ concentration in the microdomain ([Ca 2+ ] mic ) as a function of distance from the channel (r) and current through IP 3 R (I IP3R ) using the following equation [46,47]: where I IP3R = 0.05 pA [48], [Ca 2+ ] rest = 100 nM, and D Ca = 223 µm 2 /s [49] is the current flowing through an open IP 3 R, [Ca 2+ ] C at rest, and diffusion coefficient of Ca 2+ . λ is given by B T is the total concentration of Ca 2+ binding protein with K d = 2.0 µM [50]. Figure S5 shows examples of Ca 2+ concentration in the microdomain given by Equation (11) as a function of distance from the IP 3 R at different buffer concentrations.
Thus, the full model used for mitochondrial function at the single mitoplast level consists of Equations (7) (with Equations (3), (8) and (10)), (11), Equations (S1)-(S4). Equation (3) is multiplied by Equation (9) to incorporate the voltage dependence of V MCU . We remove the parameter δ in Equation (7) for simulating the single mitoplast because V MCU in the whole-cell experiments was estimated based on the dynamics of extra-mitochondrial Ca 2+ concentration, and we needed to factor-in the relative size of mitochondria with respect to the cytoplasm. In case of single mitoplast experiments, V MCU is estimated directly from the current flowing into the mitoplast through MCU and such factoring is not required.
Like the whole-cell model, the maximum flux through Na + /Ca 2+ exchanger at the single mitoplast level is determined by fitting the resting value of [Ca 2+ ] M given by the model to the observed values ( Figure S6). The values for V NaCa,max in WT, MICU1 KO, and MICU2 KO are given in the legends of Figure S6. The fit results in a significantly large V NaCa,max in MICU1 KO mitoplasts as compared to WT mitoplasts. This is in line with recent observations from whole-genome data from Alzheimer's disease patients where significant upregulation of genes encoding Na + /Ca 2+ exchangers is observed when compared to age-matched healthy individuals (see also the Discussion section below) [51].  Figure S6, where we consider a 24 nm wide microdomain [52]. A comparison between mitochondrial dynamics at the whole-cell ( Figure S3) and single mitoplast ( Figure S6) levels reveals four key differences. First, instead of hyperpolarization as in the case of whole-cell model, the Ca 2+ uptake leads to the depolarization of IMM at the single mitoplast level ( Figure S6B)-in line with the observations in Ref. [12]. Second, in the resting conditions, IMM in the MICU1 KO cells is more depolarized than WT and MICU2 KO cells, which is also consistent with observations [12]. Third, [ATP] M under resting conditions is lowest in MICU1 KO cells followed by MICU2 KO and WT cells ( Figure S6D). Finally, the mitochondrial Ca 2+ in MICU1 KO cells results in [ATP] M depletion due to the large depolarization that overshadows the rise in [NADH] M ( Figure S6C). Furthermore, the Ca 2+ influx in MICU1 KO mitoplasts results in significantly large depolarization of IMM as compared to the other two conditions. As we discuss in the next section, this is also in line with recent observations [18].
Csordás et al. found that the average spatial width of the microdomain is~24 nm in RBL-2H3 cells [52]. However, significant variability exists in the microdomain width where its value can vary from a few nanometers to several hundred nanometers [53,54].   We repeat the simulations in Figure S6 varying the duration for which the IP3R channel is open and analyze the changes in the mitochondrial variables as [Ca 2+ ]mic rises. As is clear from Figure 5  We repeat the simulations in Figure S6 varying the duration for which the IP 3 R channel is open and analyze the changes in the mitochondrial variables as [Ca 2+ ] mic rises. As is clear from Figure 5  To investigate this counterintuitive result further, we look at the three fluxes that control [NADH]M in the model (Eq. S1). As expected, the pyruvate dehydrogenase-catalyzed reaction (Eq. S8) leading to NADH production increases with [Ca 2+ ]M ( Figure S7A) [55,56]. However, the resting [Ca 2+ ]M in MICU1 KO cells (~0.258 μM) is in the range where the reaction is almost saturated already and any further increase in [Ca 2+ ]M would lead to little change in this reaction. The second contribution to [NADH]M in the model comes from the malate-aspartate shuttle (MAS) [57,58]. The model takes into account the activation of two mammalian carriers aralar and citrin involved in MAS by moderate cytosolic Ca 2+ increases (Eq. S9). By activating the Krebs cycle, the rise in [Ca 2+ ]M would lead to a decrease in the amount of α-ketoglutarate, a key-metabolite of the MAS [57,58]. This limiting step is taken into account in the model [36] and the flux due to MAS shuttle decreases with [Ca 2+ ]M and is almost at the minimum value in MICU1 KO cells ( Figure S7B) where the resting [Ca 2+ ]M ~ 0.26 μM ( Figure S6A). Finally, in the model, the oxidation of NADH in the Electron Transport Chain and the coupled extrusion of protons from the mitochondria are incorporated in a single equation (Eq. S10) [36,38,[59][60][61][62][63]. This flux declines for more hyperpolarized IMM, since it is difficult to pump protons against a large potential gradient and vice versa ( Figure S7C  To investigate this counterintuitive result further, we look at the three fluxes that control [NADH] M in the model (Equation (S1)). As expected, the pyruvate dehydrogenase-catalyzed reaction (Equation (S8)) leading to NADH production increases with [Ca 2+ ] M ( Figure S7A) [55,56]. However, the resting [Ca 2+ ] M in MICU1 KO cells (~0.258 µM) is in the range where the reaction is almost saturated already and any further increase in [Ca 2+ ] M would lead to little change in this reaction. The second contribution to [NADH] M in the model comes from the malate-aspartate shuttle (MAS) [57,58]. The model takes into account the activation of two mammalian carriers aralar and citrin involved in MAS by moderate cytosolic Ca 2+ increases (Equation (S9)). By activating the Krebs cycle, the rise in [Ca 2+ ] M would lead to a decrease in the amount of α-ketoglutarate, a key-metabolite of the MAS [57,58]. This limiting step is taken into account in the model [36] and the flux due to MAS shuttle decreases with [Ca 2+ ] M and is almost at the minimum value in MICU1 KO cells ( Figure S7B) where the resting [Ca 2+ ] M~0 .26 µM ( Figure S6A). Finally, in the model, the oxidation of NADH in the Electron Transport Chain and the coupled extrusion of protons from the mitochondria are incorporated in a single equation (Equation (S10)) [36,38,[59][60][61][62][63]. This flux declines for more hyperpolarized IMM, since it is difficult to pump protons against a large potential gradient and vice versa ( Figure S7C and smaller increase in response to cytosolic Ca 2+ in MICU1 KO cells is consistent with observations where a larger concentration of NADH in resting state and smaller change in response to Histamine-induced rise in intracellular Ca 2+ was observed in primary fibroblasts with MICU1 mutations in individuals with a disease phenotype characterized by proximal myopathy, learning difficulties, and a progressive extrapyramidal movement disorder as compared to those from age-matched healthy individuals [28]. The difference between the changes in [NADH] M and [ATP] M in MICU1 KO mitoplasts with respect to WT cells grow as we increase the open duration of IP 3 R, indicating that, in the face of prolonged [Ca 2+ ] C increases, MICU1 KO cells would be more vulnerable to damage due to the lack of ATP to perform necessary cell functions. To test whether adding Ca 2+ buffer to the cytoplasm would restore mitochondrial function in MICU1 KO cells, we vary the buffer concentration from 0 to 1000 µM. Like the results shown in Figure 4, the peak value of [Ca 2+ ] M decreases significantly in WT and MICU2 KO cells as we increase the buffer concentration ( Figure S8) Figure S8 (dashed lines in Figure S9).

Discussion
While the sigmoidal dependence of mitochondrial Ca 2+ uptake on cytosolic Ca 2+ was observed 60 years ago [64,65], details about how MICU1, MICU2, and EMRE define this relationship by regulating the gating of MCU began to emerge recently [8,66]. The functional importance of these regulatory proteins is further highlighted by their crucial role in various pathologies [27][28][29][30][31][32]. A key challenge for future experiments is to investigate how these regulatory proteins modulate mitochondrial Ca 2+ influx through MCU in intact cells under physiological and pathological conditions, and how the interaction between MCU and regulatory proteins affects the bioenergetics of mitochondria in particular and cell signaling in general. This task is further complicated by the fact that mitochondrial Ca 2+ uptake depends on the mitochondrial organization in the cell, the size of individual mitochondrion, the spatial positioning of a single mitochondrion with respect to other Ca 2+ -releasing organelles and channels, and various Ca 2+ buffering proteins inside and outside the mitochondria. To make matters worse, mitochondria are in a constant state of transporting, fission, fusion, and fragmentation as the local Ca 2+ gradients and metabolic demands change dynamically during cell functioning. Existing experimental techniques are not equipped to address the majority of these and other related issues. Data-driven computational models closely reproducing key observations about MCU function and its dependence on various variables and modulators offer a viable alternative.
Several models for mitochondrial Ca 2+ signaling have been developed over the past (for example, see [36,38,47,59,[61][62][63][67][68][69][70]). However, they are mostly based on Ca 2+ uptake in permeabilized cell cultures or isolated mitochondria suspended in aqueous solutions. These models do not properly incorporate the fine details of MCU's gating kinetics, which could lead to erroneous conclusions [70]. The model in [71] links the biophysical properties of MCU activity to the mitochondrial Ca 2+ uptake at the whole-cell level using Michaelis-Menten type relationship, but it does not incorporate the regulation of MCU by MICU1, MICU2, or EMRE. Similarly, the model in [72] investigates the mitochondrial Ca 2+ dynamics in cells and isolated mitochondria in suspensions, but does not incorporate the effects of these proteins. Our model not only reproduces many key observations about MCU function at the single mitochondrion and whole-cell levels but also incorporates the role of mitochondrial membrane potential, MICU1, MICU2, and EMRE at both scales. Our analysis reaffirms that the three regulatory proteins together with MMP, the spatial positioning of the mitochondria with respect to the Ca 2+ releasing channels or organelles, and the spatiotemporal dynamics of Ca 2+ concentration and buffering proteins make for a robust and sophisticated signaling machinery for regulating mitochondrial metabolic function.
Our model provides an excellent platform for investigating several of the issues outlined above. As an example, we investigate the interaction of a single IP 3 R with MCU in a microdomain and show that the separation between the two channels, the concentration of Ca 2+ buffer in the microdomain, the open duration of IP 3 R, and the regulatory proteins all play key roles in mitochondrial metabolic function. Interestingly, our microdomain model shows that, while Ca 2+ release through IP 3  Increasing the open duration of IP 3 R also has an opposite effect in WT (and MICU2 KO) and MICU1 KO mitochondria where we observe an increasingly positive and negative change in ATP, respectively. These results are in line with reports claiming that the human loss-of-function mutations in MICU1 results in early-onset neuromuscular weakness, impaired cognition, extrapyramidal motor disorder, and perinatal lethality [27][28][29][30][31] through mitochondrial Ca 2+ overload and impaired bioenergetics [27]. The model also leads to two key observations that warrant further discussion. First, the Ca 2+ influx in MICU1 KO mitoplasts results in significantly larger depolarization of IMM. This is in line with recent observations [18]. Specifically, HEK293 cells expressing two mutant MCU channels individually were exposed to repeated 5 to 7 µM boluses of cytosolic Ca 2+ that were 150 s apart. While mitochondria expressing WT MCU rapidly took up Ca 2+ in response to 9-12 successive boluses, those in cells expressing mutant MCU failed to take up Ca 2+ after only two boluses and exhibited a complete loss of IMM potential. Interestingly, both mutants used in these experiments results in the loss of matrix Ca 2+ -mediated gatekeeping (inhibition) like MICU1 KO.
Another key observation from the model is the significantly larger V NaCa,max in MICU1 KO mitochondria as compared to WT. This came out naturally from the fit as the Na + /Ca 2+ exchangers have to compensate for larger Ca 2+ influx in MICU1 KO mitochondria under resting conditions. Remarkably, a recent analysis of mitochondrial Ca 2+ -related genes in 25 publicly available microarray and RNA-Sequencing datasets revealed a significant upregulation of Na + /Ca 2+ exchanger-encoding gene SLC8B1 in Alzheimer's disease patients when compared to age-matched healthy individuals [51]. The authors attributed this to the counteracting effect in the diseased human brain to avoid mitochondrial Ca 2+ overload. Along these lines, upregulating the expression level of SLC8B1 gene rescued mitochondrial Ca 2+ -induced pathology in 3xTg Alzheimer's disease mouse model [73].
Given that kinetic data on the single channel function of MCU is slowly emerging [5,6,11,74], our model is formulated so that such future observations can be incorporated. For example, if time traces representing the gating state of MCU as a function of time is available, one can extract the transition rates between different states by fitting the following log-Likelihood function to the traces (see [34,75] for details): log(f(tc 1 , to 1 , tc 2 , to 2 , . . . . . . tc n , to n )) = log π C exp(Q CC tc 1 ) Q CO exp(Q OO to 1 ) Q OC exp(Q CC tc 2 ) Q CO exp(Q OO to 2 ) . . . .
where π C is a vector of the initial probabilities of all close states being occupied at equilibrium, to i and Alternatively, if the distributions of dwell-times in open or close states are available, one can extract the transition rates by fitting the following functions to such distributions, respectively [34,75]: where π O is a vector of initial probabilities of all open states being occupied at equilibrium, u C is a vector of all ones with length equal to the number of close states in the model, Q CO is a (6 × 2) matrix of the transition rates from all close to all open states, and Q CC is a (6 × 6) matrix of the transition rates from all close to all close states etc. As a toy example, we used the double-exponential behavior of open dwell-time histogram observed in [11]. Specifically, in these experiments, MCU and the regulatory proteins were inserted into a planar bilayer membrane and electrophysiological recording of single channel activity was carried out. In one of these experiments ( Figure S3B in [11]), the open dwell-time distribution could be fitted with double exponential with time constants τ 1 = 12.06 ms and τ 2 = 161.92 ms. We used these time constants to generate experimental dwell-time distribution (green bars in Figure 6A) and fit Equation (14) to it. The model fit to the data is shown by the dashed line in Figure 6A. The transition rates between different states from the fit are given in section "Extracting single channel gating parameters from dwell-time distributions" of the Supplementary Information Text. The rates obtained from the fit are then used to stochastically simulate the gating of single MCU channel at different [Ca 2+ ] C and [Ca 2+ ] M values using the method previously developed [35,76,77]. Consistent with the results at the whole-cell and single mitochondrion levels, the single channel activity of MCU increases progressively as we increase [Ca 2+ ] C ( Figure 6B,C). Furthermore, the frequency of the channel visiting the open states in general, and O 44 state in particular increases significantly as we increase [Ca 2+ ] M from 0.4 µM ( Figure 6B) to 1 µM ( Figure 6C). Cells 2020, 9, x FOR PEER REVIEW 17 of 22 f C t C =π C exp(Q CC t C )Q CO u O . (15) where is a vector of initial probabilities of all open states being occupied at equilibrium, is a vector of all ones with length equal to the number of close states in the model, QCO is a (6 × 2) matrix of the transition rates from all close to all open states, and QCC is a (6 × 6) matrix of the transition rates from all close to all close states etc.
As a toy example, we used the double-exponential behavior of open dwell-time histogram observed in [11]. Specifically, in these experiments, MCU and the regulatory proteins were inserted into a planar bilayer membrane and electrophysiological recording of single channel activity was carried out. In one of these experiments ( Figure S3B in [11]), the open dwell-time distribution could be fitted with double exponential with time constants τ1 = 12.06 ms and τ2 = 161.92 ms. We used these time constants to generate experimental dwell-time distribution (green bars in Figure 6A) and fit Eq. (14) to it. The model fit to the data is shown by the dashed line in Figure 6A. The transition rates between different states from the fit are given in section "Extracting single channel gating parameters from dwell-time distributions" of the Supplementary Information Text. The rates obtained from the fit are then used to stochastically simulate the gating of single MCU channel at different [Ca 2+ ]C and [Ca 2+ ]M values using the method previously developed [35,76,77]. Consistent with the results at the whole-cell and single mitochondrion levels, the single channel activity of MCU increases progressively as we increase [Ca 2+ ]C ( Figure 6B, C). Furthermore, the frequency of the channel visiting the open states in general, and O44 state in particular increases significantly as we increase [Ca 2+ ]M from 0.4 μM ( Figure 6B) to 1 μM ( Figure 6C).

Conclusions
Mitochondrial Ca 2+ overload either due to MICU1 mutations [27][28][29][30][31], MCU upregulation [51], or impaired function of Na + /Ca 2+ exchangers [73] plays a causal role in various diseases. For example, recent observations show that it not only contributes to the progression of disease but also precedes neuropathology in different animal models of Alzheimer's disease [51,73]. Mitochondrial Ca 2+ overload was shown to accelerate memory decline, amyloidosis, and tau pathology by promoting the production of reactive oxygen species, metabolic dysfunction, and upregulation of β secretase expression, and neurodegeneration [51,73]. Furthermore, preventing mitochondrial Ca 2+ overload was proven sufficient to impede Alzheimer's disease-associated pathology and memory loss [51,73].

Conclusions
Mitochondrial Ca 2+ overload either due to MICU1 mutations [27][28][29][30][31], MCU upregulation [51], or impaired function of Na + /Ca 2+ exchangers [73] plays a causal role in various diseases. For example, recent observations show that it not only contributes to the progression of disease but also precedes neuropathology in different animal models of Alzheimer's disease [51,73]. Mitochondrial Ca 2+ overload was shown to accelerate memory decline, amyloidosis, and tau pathology by promoting the production of reactive oxygen species, metabolic dysfunction, and upregulation of β secretase expression, and neurodegeneration [51,73]. Furthermore, preventing mitochondrial Ca 2+ overload was proven sufficient to impede Alzheimer's disease-associated pathology and memory loss [51,73]. A data-driven model such as ours will help understand the role of different mechanisms regulating mitochondrial Ca 2+ uptake in Alzheimer's disease and many other diseases where mitochondrial Ca 2+ overload is believed to play a major role [78,79]. Our model not only incorporates several key observations about MCU function under WT, MICU1 KO, and MICU2 KO conditions at the whole-cell and single mitochondrion levels, but also leaves room for incorporating future observations and providing information about experimentally inaccessible parameters at the single MCU level such as the transition rates between different conducting states of the channel, etc.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4409/9/6/1520/s1. Supplementary Methods (Equations (S1)-(S15)); Figure S1: Using a different set of parameters to fit the model the MCU Ca 2+ uptake rate in MICU2 KO cell cultures, Figure S2: Considering the state with two and one Ca 2+ bound to MICU1 and MICU2 respectively to be an open state results in higher Ca 2+ uptake rate at low [Ca 2+ ] C , Figure S3: Changes in mitochondrial variables in response to a step-like increase in [Ca 2+ ] C at the whole-cell level, Figure S4: Schematic of the mitochondrial Ca 2+ uptake model at the single mitochondrion level, Figure S5: Ca 2+ concentration in the microdomain ([Ca 2+ ] mic ) due to the opening of a single IP 3 R as a function of distance from the channel, Figure S6: Changes in mitochondrial variables in response to the opening of a single IP 3 R channel at the single mitoplast level, Figure S7