The effects of calcium on mitochondrial energetics span from the beneficial in the low range (10’s of nmol/mg) to the catastrophic in the high range (greater than 500 nmol/mg). In the low range, matrix calcium activates TCA cycle dehydrogenases and other matrix enzymes that increase the rate of proton pumping and ATP synthesis. In the high range, calcium overload leads to mitochondrial permeability transition and membrane rupture. In this state, mitochondria are unable to maintain a proton electrochemical gradient and thus unable to generate ATP. However, in the intermediate range, calcium leads to a depression in ATP synthesis rates, despite maintaing membrane integrity and being able to maintain a proton electrochemical gradient. Fig 1A shows the prototypical response of isolated mitochondrial respiration to varying amounts of calcium exposure. As the amount of calcium taken up by the mitochondria increases, the ADP-stimulated respiration rates are lowered in a titratable manner. In these experiments, ATP synthesis remains coupled to oxygen consumption. The exact cause of this reduction in ATP synthesis rates in calcium loaded mitochondria is not precisely known. Herein, we test a combination of experimentally and computationally derived hypotheses to determine why ATP synthesis is inhibited in the calcium overloaded state.
Based on the data shown in Fig 1, the following mechanisms explain the calcium inhibition data: i) lower membrane potential caused by sodium/calcium cycling lowers the driving force for ATP production; ii) a subpopulation of mitochondria permeabilize and reduce the total number of mitochondria capable of synthesizing ATP; iii) calpain activation leads to proteolytic loss of electron transport and/or ATP synthesis function; iv) available ADP and/or phosphate for phosphorylation is depleted or reduced by calcium binding, incorporation into calcium phosphate granules, or net loss; and iv) direct inhibition of calcium on mitochondrial processes critical for ATP production. We employed both experimental and computational strategies to test these hypotheses.
The first hypotheses states that the lower membrane potential caused by sodium/calcium cycling lower the thermodynamic driving force, and hence, the rate for ATP synthesis. And while the membrane potential, as shown in Fig 1B, is lower after the CaCl2 bolus was given due to sodium/calcium cycling, the membrane potential during ATP synthesis is the same for all calcium conditions. This rules out that calcium-dependent reduction in respiration rates during ADP-stimulated respiration is due to a lowering driving force for ATP production. Thus, hypothesis i) is disproved.
The second hypothesis involves subpopulations of mitochondria responding differently to the calcium challenges. Isolated mitochondria in suspension do not undergo calcium-dependent permeabilization all at once. It is believed that the individual calcium tolerance of each mitochondrion in a population of suspended mitochondria falls within a distribution , so that the more susceptible mitochondria undergo permeability transition first. Upon permeabilization, they release their calcium content which is taken up by more calcium-resistant mitochondria until they reach their calcium limit, and so on. This cascade of events might explain the data given in Fig 1A. However, both the membrane potential data in Fig 1B and the calcium uptake data in Fig 1C dispute this notion. The membrane potential data reveal that energized status of the population of mitochondria is stable, and the calcium uptake data show that calcium uptake is robust with no detectable permeabilization. This disproves hypothesis ii).
The next hypothesis argues that calpain activation is the reason why calcium lowers ADP-stimulated respiration. Calpains are intracellular calcium-activated cysteine proteases present in nearly all vertebrate cells . They are often, but not always, associated with subcellular organelles. The physiological function of calpains are not fully understood but they likely play major roles during autolysis. Calpains remain inactive until calcium increases to sufficient levels. In the ischemic myocyte, the dysregulation of cytosolic and mitochondrial calcium is thought to be central to calpain activation. Studies have shown that mitochondrial calpain activation does cleave specific sites on complex I and ATP synthase which has been linked to poor substrate oxidation and ATP production [33, 34, 37]. So we decided to test whether or not calpains played any sort of role in the calcium-induced depression of ADP-stimulated respiration. We experimentally tested this hypothesis by incubating isolated mitochondria with calpain inhibitors prior to calcium exposure. If calpains were responsible for this phenomenon, then inhibiting their proteolytic activity should rescue the observed calcium triggered mitochondrial dysfunction. Our results shown in Fig 2 demonstrate that calpain inhibitors do not relieve this inhibition and thus are not involved in the observed phenomenon. Longer incubation times did not improve mitochondrial function (S1 File). Hypothesis iii) is disproved.
Fig 2. Calcium inhibition on ADP-stimulated respiration is independent of calpain activity.
Conditions are like those given in Fig 1. Briefly, 0.1 mg/ml mitochondria are incubated in the presence of 10 μM calpain inhibitor and substrates for five minutes before a bolus of CaCl2 is added followed by bolus of 500 μM ADP five minutes after the CaCl2 bolus. Calpain inhibitors show no effect on the calcium-dependent inhibition of ADP-stimulated respiration. In all experiments, data are presented as mean +/- standard deviation from three or more biological replicates. Individual data are presented as gray dots. There were no statistically significant differences within calcium treatment groups at a 0.05 alpha level. Control data for each calpain inhibitor were combined.
In order to test the remaining hypotheses, we resorted to utilizing a computational approach. We started with the bioenergetics model from Bazil et al.  and added calcium handling, matrix calcium buffering, and updated a few other reactions to improve fits to the original data . See the Supporting Information (S1 File) for details. In brief, very simple mitochondrial calcium uniporter (MCU) and sodium calcium exchanger (NCLX) rate expressions were used with the sodium hydrogen exchanger rate expression from Bazil et al. . We first started with more biophysically detailed MCU  and NCLX  rate expressions; however, we encountered several issues which precluded us from incorporating these more detailed expressions into the current model. First, the detailed MCU rate expression showed saturation with respect to respiration. This aberrant behavior was due to a calcium dissociation constant for the MCU that was too low causing the calcium current to saturate and produce a shark fin like respiratory dynamic. From Fig 1A, it can clearly be seen no such saturation-like respiratory response is evident. As such, we used a much simpler rate expression with a higher Vmax and KD for calcium that is more in line with the known characteristics of the channel . In addition to the MCU, the detailed NCLX flux expression was overly complex and prevented adequate fits to the data. In this case, the steep proton dependency was determined to be the main problem. As with the MCU problem, a simpler rate expression was adequate to fit the experimental data. We also included the model of the mitochondrial calcium sequestration system based on data from Blomeyer et al.  described in detail below. The last major changes to the model was to the lumped TCA cycle rate expression and updating the complex III rate expression. We noticed the rate expression published with the original model failed to adequately simulate the respiratory behavior after calcium loading. This was due to calcium-dependent alkalization of the matrix not being properly compensated for in the expression and was alleviated by slightly altering the NADH feedback component. In addition to changing the TCA cycle rate expression, we replaced the complex III rate expression with a new one that includes dimer functionality from Bazil et al. . The new rate expressions are given in the Supporting Information (S1 File).
Fig 3 shows the updated schematic of the bioenergetics model with new rate expressions and components, the new updated model fits to the experimental data it was originally parameterized on, and the parameter correlations shown as a heat map for the adjustable parameters given in Table 1. As shown in the center four panels in Fig 3, the model fits the data equally well as the original model, and in several cases even better. In particular, the updated model fits are more faithful to the data for the NADH and ADP data. The parameter correlation heatmap shown in the right panel of Fig 3 reveals the model parameters are relatively uncorrelated except for strong correlations between the MitoDH and electron transport related parameters.
Fig 3. Mitochondrial bioenergetics model diagram, calibration, and parameter correlations.
The upper panel shows an updated model diagram of the Bazil et al. model  that includes calcium handling reactions. The updated model was fit to the data from Vinnakota et al. . The lower left four panels show model simulations faithfully reproducing the experimental data. The lower right panel shows the parameter correlations represented as a heat map.
Local sensitivity coefficients are normalized and averaged with all model outputs aligned with the experimental data. Only non-zero sensitivity coefficients were averaged. The leak activity for the model simulations of the guinea pig data (calcium-dependent inhibition data) was lowered by 40% to account for the tighter respiratory coupling that these mitochondria possess.
There is lot of evidence that calcium is reversibly sequestrated as calcium phosphate granules to maintain energy homeostasis [12, 13]. Although we have implicitly modeled granule formation and calcium storage using an empirical model , this approach required an artificially high proton buffering capacity to prevent excess matrix alkalization. This over alkalization was corrected by explicitly modeling calcium phosphate granule formation and including proton release as a part of this process . A general reaction formula for calcium phosphate formation is depicted in Eq 1. This equation is charge and atom balanced and works for many types of calcium phosphate stoichiometries including tricalcium phosphate, hydroxyapatite, and octacalcium phosphate. In the case that 3m-2n < 0, the H in the calcium phosphate complex is understood to be OH with a stoichiometry of 2n-3m.
Assuming an additional prototypical buffering component, the total matrix calcium buffering power that includes the formation of calcium phosphate is shown in Eq 2.
In Eq 2, we use the concentration of hydrogen phosphate, as opposed to the phosphate ion, for maximum compatibility and simplicity. The bioenergetics model does not include the phosphate ion as a state variable because this ion only exists in appreciable quantities under very basic, non-physiological conditions. Including the phosphate ion would not change the results but require addition state variables and parameters to be added to the model.
Fig 4 shows the relationship between mitochondrial calcium buffering power and matrix free calcium with the explicit representation of calcium phosphate granule formation compared to the experimental data from Blomeyer et al. . The model consists of two essential components. The first component characterizes prototypical calcium buffering caused by binding to proteins and lipids not explicitly accounted for in the model. The second component represents the formation calcium phosphate granules. Table 2 gives the parameters used to simulate the data given in Fig 4. The calcium phosphate complexation constant is not a solubility product constant. It is a parameter that governs the matrix calcium buffering power. For this data, either tricalcium phosphate, hydroxyapatite, or octacalcium phosphate formation fit the data with near equality. As shown in the inset of Fig 4, the major difference between the three types of calcium phosphates is that higher buffering powers for matrix calcium concentrations exceeding 10 μM are obtained as calcium stoichiometric coefficient increases. For simplicity, we chose to model tricalcium phosphate formation. In addition, we assume calcium phosphate precipitate formation is rapid relative to calcium uptake and efflux. Further experiments are required to determine which calcium phosphate species constitute the main precipitate.
Fig 4. Mitochondrial calcium buffering model.
The mitochondrial calcium buffering data from Blomeyer et al.  for the conditions most similar to the conditions in this study was used to fit the new calcium phosphate buffering system using an explicit representation of calcium phosphate granule formation. In the inset, the buffering power for the tricalcium phosphate, hydroxyapatite, and octacalcium phosphate are shown for higher matrix calcium concentrations. These simulations were performed with the following matrix conditions: pH 7.4, [K+] = 150 M, [Na+] = 5 mM, [Mg2+] = 0.5 mM, and [Pi]tot = 10 mM.
We then proceeded to simultaneously fit the model to the data shown in Fig 1 by allowing some of the calcium related parameters to change. For these simulations, we ignored the calcium-dependent reduction in ADP-stimulated respiration to get a baseline model capable of simulating the calcium handling data. To minimize the number of adjustable parameters, we selected the most important parameters with respect to calcium handling using sensitivity analysis conducted with estimated parameter values obtained from the literature. These parameters are given in Table 3 and were determined to be the MCU activity, MCU calcium dissociation constant, NCLX activity, and the NCLX calcium dissociation constant. All other model parameters used for the simulations are given in the Supporting Information (S1 File). The most sensitive parameter is the MCU activity. This parameter has a high, normalized local sensitivity coefficient relative to the others. The next highest ranked sensitivity is for the MCU calcium dissociation constant. The remaining two parameters, the NCLX activity and the NCLX calcium dissociation constant have near equal sensitivity values. The parameters were relatively correlated (between 0.7 and 1) which indicates additional information is needed to uniquely identify these values. Also, the NCLX parameters are dependent on the matrix calcium buffering parameters given in Table 2. However, these correlations do not impact the model results in this study.
Local sensitivity coefficients are normalized and averaged with all model outputs are aligned with the experimental data. Only non-zero sensitivity coefficients were averaged.
With the model calibrated, we began to computationally test the remaining hypotheses. Fig 5 shows that the reduction in ADP-stimulated respiration is not due to a reduction in available ADP and/or phosphate for oxidative phosphorylation. Even the highest calcium bolus only leads to less than 1% and 2% reduction in available ADP and phosphate, respectively. This makes it highly unlikely that this is the cause of the calcium-dependent reduction in respiration shown in Fig 1. However, these simulations do not rule out the possibility that ADP is physically incorporated into the calcium phosphate granules. That remains a possibility to be tested, thus the hypothesis is not fully disproven. We will revisit this hypothesis below.
Fig 5. Calcium complexation with ATP synthase substrates does not explain the calcium-induced inhibition of ADP-stimulated respiration.
A) Simulations of the percent of matrix ADP during the calcium loading and ADP bolus experiments given in Fig 1. The times when the CaCl2 and ADP boluses, as well as, when the environment becomes anoxic are marked. These times are the same for all panels. B) The fraction of matrix ADP chelated with calcium for the same simulation. C) The simulated concentrations of matrix Pi during the experiments. D) The corresponding simulations of solvated calcium phosphate species (not calcium phosphate precipitate) during the experiments. The simulation protocols were as faithful to the experimental protocols as reasonably possible.
We next tested the remaining hypothesis by adding a calcium-dependent Michaelis Menten-like inhibitory function to the main reactions involved in ADP-stimulated respiration. We opted to use a simple inhibitory function as opposed to more non-linear Hill type functions for simplicity. For this, we fit the adjustable parameters given in Table 3 to simulate the calcium handling components of the mitochondrial bioenergetics model for the experimental conditions described in Fig 1. All other parameters were fixed and given in detail in the Supporting Information (S1 File). Table 4 summarizes the results and shows that the best mechanism that matches the data is the inhibition of complex I as a function of the concentration calcium phosphate granules. The next best mechanism is the inhibition of complex III as a function of calcium phosphate granule concentration. All other conditions did not lead to adequately good fits to the data based on the normalized least squares values and visual inspection. The relatively small inhibition constants for the inorganic phosphate carrier and F1FO ATP synthase are expected, as these rate expressions have high activities to put them in a near equilibrium state. More complex inhibition functions involving multiple reactions were not considered. These results show that the most likely site of inhibition is at complex I or III and the calcium phosphate granule concentration is the likely cause.
Fits were quantified using a normalized least squares value which was computed by summing the square of the difference between model and data relative to the data and divided by the total number of data points.
The model results using the complex I inhibition mechanism compared to the experimental data is shown in Fig 6. As shown, the model faithfully reproduces the experimental data. The model captures the calcium-dependent stimulation of respiration caused by sodium/calcium cycling and the inhibitory effect of calcium phosphate precipitates on ADP-stimulated respiration. In addition, the model reveals why the mitochondrial membrane potential is more polarized during leak state respiration when 1 mM EGTA is present. With the model, we can mechanistically explain this observation. The reason why the membrane potential is higher in this respiratory state when 1 mM EGTA is present, the matrix pH is closer to the buffer pH. This lowers the ΔpH component of the proton motive force and leads to an increase in the membrane potential. The matrix alkalization when 1 mM EGTA is not present is due to uptake of the residual calcium in the buffer. This residual calcium contaminate comes from reagent impurities mainly coming from potassium chloride and potassium phosphate salts. This is a small, but important, validation of the model ability to accurately simulate mitochondrial physiology.
Fig 6. Model simulations compared to the experimental data of the calcium-induced inhibition of ADP-stimulated respiration.
A) Model simulations of respiration for the different calcium conditions for the protocol detailed in Fig 1. The inset shows the model simulations for leak and sodium/calcium cycling respiratory states in more detail. B) Model simulations of membrane potential dynamics. The model output for membrane potential was converted to a percent basis from a theoretical 200 mV maximum so that it was comparable to the normalized TMRM data. C) Simulations of the stimulatory effect of calcium uptake on respiration. D) Simulations of the buffer calcium dynamics. E) Simulations of the matrix calcium dynamics. F) Simulations of the matrix pH dynamics. These simulations were carried out assuming that Complex I was inhibited by the concentration of calcium phosphate granules using the inhibitory constant given in Table 4. To remain faithful to the experimental conditions, we added an ATPase reaction that represents the typical ATPase contaminate that is present in all isolated mitochondrial preps. For details, see the Supporting Information (S1 File). Error bars are given as 95% confidence intervals. Individual data are represented as dots. The colors for each panel represent the conditions given in the legend of panel B.