Calcium oscillations optimize the energetic efficiency of mitochondrial metabolism

Summary Energy transduction is central to living organisms, but the impact of enzyme regulation and signaling on its thermodynamic efficiency is generally overlooked. Here, we analyze the efficiency of ATP production by the tricarboxylic acid cycle and oxidative phosphorylation, which generate most of the chemical energy in eukaryotes. Calcium signaling regulates this pathway and can affect its energetic output, but the concrete energetic impact of this cross-talk remains elusive. Calcium enhances ATP production by activating key enzymes of the tricarboxylic acid cycle while calcium homeostasis is ATP-dependent. We propose a detailed kinetic model describing the calcium-mitochondria cross-talk and analyze it using nonequilibrium thermodynamics: after identifying the effective reactions driving mitochondrial metabolism out of equilibrium, we quantify the mitochondrial thermodynamic efficiency for different conditions. Calcium oscillations, triggered by extracellular stimulation or energy deficiency, boost the thermodynamic efficiency of mitochondrial metabolism, suggesting a compensatory role of calcium signaling in mitochondrial bioenergetics.


Highlights
A thermodynamically consistent model of Ca 2+ mitochondria cross-talk is established The efficiency of Ca 2+regulated mitochondrial ATP synthesis is quantified Induction of Ca 2+ spikes by energy deficiency or stimulation increases efficiency Oscillatory and not stationary Ca 2+ concentration boosts mitochondrial efficiency

INTRODUCTION
Life relies on permanent conversions between different forms of energy, a phenomenon referred to as energy transduction.A wide range of cellular processes are fueled by the chemical energy stored in adenosine triphosphate (ATP), but the compartmentalization of eukaryotic cells also enables the storage of potential energy across the membranes of organelles. 1 Energy transduction is mediated by enzymes and pumps driven in a nonequilibrium thermodynamic manner by the hydrolysis of ATP, chemical gradients or membrane potentials.
In optimal scenarios where transduction is fully efficient, the input energy is completely transformed into usable work.However, biological processes are typically accompanied by entropy production, i.e., dissipation of energy in the form of heat and/or chemical waste that is unusable for transduction. 2For example, the action of many transmembrane ionic pumps transporting ions against their concentration gradient is often based on catalyzing the hydrolysis of ATP.The chemical energy released by hydrolysis is partly used to drive ionic transport while another part is dissipated.In the extreme case of pump uncoupling, also known as ''slippage'', all the energy of ATP hydrolysis is dissipated without any ion transport. 3ifferent nonequilibrium kinetic models have been developed to account for energy loss in pumps [4][5][6][7] but have only provided limited insights into energetic costs at the pathway level.New approaches based on metabolic network reconstruction and nonequilibrium thermodynamics are gradually emerging to rationalize the energetic costs of cellular processes 8 including gene regulation, 9 repair mechanisms, 10,11 enzymatic catalysis, 12 information processing 13 or signaling. 14,15A framework to study energy transduction in complex open chemical reaction networks (CRN) has recently been proposed and used to study the efficiency of pathways of the central energy metabolism in the absence of regulations. 16Evaluating the efficiency of tightly coupled transduction processes, i.e., processes whose input and output currents are equal, is straightforward as it does not depend on the net reaction flux.However, when regulations come into play, this tight coupling can be lost and kinetic models become indispensable to evaluate the flux of the different processes contributing to the efficiency.
Here, we resort to such a kinetically-detailed nonequilibrium thermodynamic approach to show and quantify how active signaling can have a beneficial energetic impact on metabolism.In particular, we analyze the efficiency of the mitochondrial production of ATP via the tricarboxylic acid (TCA) cycle and oxidative phosphorylation (OXPHOS), and take into account its regulation by calcium (Ca 2+ ).Ca 2+ is a ubiquitous and

Modeling and theoretical frameworks
[27][30][31][32][33][34][35] Balanced chemical equations, detailed expressions of the reaction rates, thermodynamic forces and reference parameter values are given in Tables 1, 2, 3, and 4, respectively.The model was calibrated to reproduce key features of the Ca 2+ -metabolism cross-talk based on experimental data from cultured astrocytes, 36 to generate physiological concentrations and concentration ratios with a realistic mitochondrial membrane potential (Figures S1 and S2 and Table 5), and to be thermodynamically consistent (details in Model calibration in the STAR Methods).1][22][23] The underlying kinetic models originally aimed at capturing the essential mechanisms of the Ca 2+ -metabolism interplay and at rationalizing experimental data about the response of Ca 2+ signals to changes in mitochondrial activity (and vice versa).][39][40][41][42][43] To compute their metabolic efficiency, we analyzed mitochondria as out-of-equilibrium chemical engines (Figure 1B) satisfying the second law of thermodynamics: 38,41 T _ s = À d t G + _ w nc + _ w driv : (Equation 1) Mitochondrial metabolism constitutes an open CRN that continuously harnesses the free energy stored in buffered species (e.g., AcCoA, CoQ, O 2 , H + m ) to convert cytosolic ATP (ATP c ) from cytosolic ADP (ADP c ) while being influenced by Na + homeostasis and cytosolic processes such as Ca 2+ signaling and ATP c consumption.From a thermodynamic perspective, the synthesis of ATP c and the regulations correspond to free energy exchanges between the mitochondrial engine and its surroundings.They appear in the second law (Equation 1) as the nonconservative work rate, _ w nc , and the driving work rate, _ w driv , respectively._ w nc is the energy current maintaining the CRN out of equilibrium while _ w driv is the energy current resulting from the modification of the underlying equilibrium state by the out-of-equilibrium dynamics (Figure 1C).The difference between their sum and the variation in time of the internal Gibbs free energy of mitochondria, G, equals the free energy dissipated by the mitochondrial reactions, i.e., the entropy production rate (EPR) _ s times the absolute temperature T.

A B C
Figure 1.Representation of the model components, conceptualization of mitochondria as a chemical engine, expected work contributions to the thermodynamic efficiency in equilibrium and nonequilibrium conditions and corresponding abbreviations Balanced chemical equations, detailed expressions of the reaction rates, thermodynamic forces and reference parameter values are given in Tables 1, 2, 3, and 4, respectively.
(A) The upper part depicts the Ca 2+ (red) and ATP (blue) fluxes responsible for the cross-talk between Ca 2+ dynamics and mitochondrial metabolism.The bottom part is a detailed description of the model components.The kinetic rates for TCA cycle fluxes and processes involving exchanges across the mitochondrial membrane are respectively originating from Dudycha 24 and Magnus-Keizer models, 24À26.except for the transformation of MAL into OAA, which is described more realistically by a reversible flux. 25Here, OXPHOS corresponds to the net redox reaction resulting from the electron transport chain (Ox) and the synthesis of ATP by the F1F0-ATPase (F1).A last module, consisting of Ca 2+ exchanges across the ER membrane and cytosolic ATP hydrolysis are taken from the models from Komin et al. 26 and from Wacquier et al. 27 Controlled species (i.e., species whose concentration is assumed to be constant) are shown in grey, dynamical species in black and dashed arrows represent regulations.Processes are annotated in yellow and black boxes for mitochondrial and cytosolic/ER processes, respectively.(B) Mitochondrial metabolism is conceptualized as an open chemical engine that transforms ADP c into ATP c through a set of 2 emergent cycles split in 3 effective reactions (r1 out , r1 in and r2).Some of the controlled species involved in the internal reactions are buffered at a constant concentrations (green), while Na + (brown) and Ca 2+ (red) regulate reaction rates by activating specific enzymes or acting on the mitochondrial membrane potential.(C) Expected work contributions in equilibrium and nonequilibrium conditions.Stars denote the state of the system at different time points in nonequilibrium conditions (violet) and the corresponding underlying equilibrium state (orange).As illustrated, the non-zero driving work contribution in the oscillatory regime can modify the underlying equilibrium state of the system.Abbreviations: AcCoA -acetyl coenzyme The expressions of the thermodynamic quantities in Equation 1 are derived for mitochondrial metabolism using a topological analysis 40,41,43 of the corresponding CRN, which allowed us to identify 13 conservation laws and 2 emergent cycles (see subsection Conservation laws and emergent cycles in the STAR Methods).The conservation laws define parts of molecules that remain intact in all mitochondrial reactions and are instrumental to determine the Gibbs free energy G.The emergent cycles are split into 3 effective reactions that can be intuitively identified as mitochondrial output or input 3) (Equation 4) with the subscripts c and m referring to the location of some species in the cytosol and mitochondria, respectively.Notice that r1 out and r1 in are not independent as they have the same reaction current and that their sum constitutes the first emergent cycle.In contrast, r2 is the exact chemical equation of the second emergent cycle and has a different reaction current.The nonconservative work rate is correspondingly split into the sum of 3 contributions, equal to the product of the free Gibbs energy of reaction and the effective current of each reaction: , where _ w r1out quantifies the mitochondrial free energy output corresponding to the synthesis of ATP (Equation 2), while _ w r1in and _ w r2 quantify the mitochondrial free energy power source in a regime optimizing ATP production -Equation 3 -or the proton driving force -Equation 4-respectively.
The average thermodynamic efficiency h of mitochondria can then be calculated as h = À w output nc w input nc +w driv (Equation 5) where w input nc = w r1in + w r2 and w output nc = w r1out , and the overline denotes either steady-state quantities or averages over one period of Ca 2+ oscillations (notice that d t G = 0).Subscripts c, ER and m refer to the cytosol, the endoplasmic reticulum and the mitochondria, respectively.

Magnus and
Keizer 30 CS

Magnus and
Keizer 30 with Keizer 30 Hyd The EPR, nonconservative and driving work contributions vanish at equilibrium according to the second law of thermodynamics but take finite values in nonequilibrium regimes (Figure 1C).The dynamics of the full system (Equations 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,  21, 22, and 23) was simulated for different stimulation conditions and mitochondrial substrate concentrations (i.e., for different ½IP 3 and ½AcCoA in the simulations), which allowed for the calculation of the corresponding nonconservative and driving work contributions and, ultimately, of the efficiency of mitochondrial metabolism.

Ca 2+ -metabolism cross-talk affects the oscillation period and the production of ATP
To calibrate the kinetic model, we compared our simulation results to experimental and simulation data from the literature.More details about the calibration procedure can be found in the STAR Methods.We considered physiological conditions as well as substrate-limited and overstimulation conditions.To provide some baselines, a normal substrate level corresponds to ½AcCoA = 10 mM 44 and in the absence of external stimulation of GPCRs (i.e., unstimulated conditions), the so-called basal ½IP 3 is close to 0:1 mM. 45½AcCoA and ½IP 3 are respectively decreased and increased with respect to those baseline values to reach substrate-limited and overstimulation conditions.7][48] Overstimulation of IP 3 Rs by IP 3 can lead to energy deficiency but is above all characterized by the saturation of the oscillation period or by a high steady-state ½Ca 2+ c of several hundreds of nanomolars.As described in subsection The relation between mitochondrial efficiency and dissipation is different in fast Ca 2+ spiking regimes triggered in substrate-limited or overstimulation conditions, substratelimited and overstimulation conditions also display different efficiency-dissipation profiles.
Upon increase of ½IP 3 or decrease of ½AcCoA, a transition (bifurcation point) from steady-state to slow oscillations marks the onset of the signaling machinery.A decrease in the oscillation period is observed as ½IP 3 is increased (Figures 2A-2C) or as ½AcCoA is decreased (Figure 2C).0][51][52] They also reproduce the Ca 2+ dynamics of cultured astrocytes 36 and Xenopus oocytes 53 reported for limited availability of mitochondrial substrate.In the oscillatory regime, ½ATP c displays a maximum in dependence on ½IP 3 and ½AcCoA (Figures 2D, 2E, and S3), a feature that is also predicted by the model of Wacquier et al. 27 In our simulations, a cusp in the average of ½ATP c is additionally observed at the critical point (Figures 2D, 2E, and S3).Furthermore, our model predicts, at a more intense stimulation level of GPCRs, a slight decrease in the period as ½AcCoA is increased (Figure 2C, e.g., at ½IP 3 = 0:50 mM).This result contrasts with the results of Duchen and co-workers who reported that energized mitochondria tend to slow down Ca 2+ oscillations due to increased mitochondrial uptake and reduction of CICR. 54,55In our model and for those stimulation levels of GPCRs, the increased mitochondrial uptake is accompanied by a larger mitochondrial efflux through NCX, which maintains the activation of CICR (data not shown).Ca 2+ uptake does not decrease as DJ decreases (Figure S2) indicating that DJ has a minor impact on Ca 2+ uptake in this model, which can probably not be compared with the chemically induced collapse of DJ reported in the experiments of Duchen et al. 54,55 Most of these observations can be rationalized based on the dependence of SERCA pumps on ATP c , which enables the switch between ER and mitochondrial Ca 2+ sequestration and is a key signature of the Ca 2+ -metabolism cross-talk.More precisely, this dependence includes the ATP c consumption associated to Ca 2+ uptake into the ER and the regulation of SERCA flux (J SERCA ) by ATP c , which are expressed mathematically in the evolution equation of ½ATP c (Equation 9) and by the ½ATP c -dependent factor in J SERCA (Table 2), respectively.From a thermodynamic point of view, J SERCA contributes to the effective current of r1 out and r1 in , hence to mitochondrial thermodynamic efficiency.We stress that in physiological conditions, ½ATP c is in the millimolar range, which is about 20 times larger than the ATP c dissociation constant of SERCA
(Table 4).SERCA's activity is thus ATP c -insensitive but nevertheless impacts ½ATP c as long as ½Ca 2+ c is close to or larger than the Ca 2+ c dissociation constant of SERCA (Table 4).At basal ½Ca 2+ c (about 0:1 mM 54, 55 ) SERCA's activity is thus reduced but is the dominant Ca 2+ removal mechanism (Figure S1).
As ½IP 3 increases, more Ca 2+ is released into the cytosol through IP 3 Rs.The steady-state ½ATP c thus decreases due to a more demanding maintenance of the basal ½Ca 2+ c via SERCA pumps.At the critical ½IP 3 corresponding to the onset of oscillations, mitochondrial sequestration of Ca 2+ becomes significant, which not only relieves SERCA pumps but also enables the activation of Ca 2+ -sensitive dehydrogenases of the TCA cycle (Figure S1).These combined effects result in an increase of the average ½ATP c .Upon further increase of ½IP 3 , the activation of mitochondrial dehydrogenases by Ca 2+ reaches saturation (Figure S1) and the ATP c consumption associated to Ca 2+ homeostasis is no longer counterbalanced by the Ca 2+ -enhanced mitochondrial activity, which results in a slow decrease in average ½ATP c .Meanwhile, increasing stimulation by IP 3 favors more frequent opening of the IP 3 Rs, which results in an decrease of the oscillation period.In most mathematical models for Ca 2+ signaling and in agreement with experimental observations, the oscillation period saturates at high ½IP 3 69 and, beyond a critical ½IP 3 , Table 3. Forces of the system Transformed Gibbs free energies of reaction (D r G k ) associated to each process of the system.The indices c and m associated to D r G 0+ k indicate that this thermodynamic quantity is evaluated at cytosolic and mitochondrial pH, that is, pH = 7:2 and pH = 8:0, respectively.The value of D r G 0+ k , which also accounts for physiological ionic strength (I = 0:12 M 28 ) and ½Mg 2+ m (pMg = 3.4), was retrieved for each relevant process via Equilibrator. 29Starting from the corresponding pseudoisomer concentrations, Magnus and Keizer estimate that  oscillations disappear.The cell then exhibits a high-½Ca 2+ c steady-state, 50 as reproduced by our simulations (Figure S1).Further stimulation by IP 3 does not affect the steady-state concentrations reached after termination of the oscillations (see Figure 3B bottom for ½ATP c ; Figure S1 for ½Ca 2+ c and ½Ca 2+ m ), suggesting that IP 3 Rs have reached their maximal release rate and contribute to saturation effect.
The impact of AcCoA level on Ca 2+ oscillations and the presence of the cusp behavior in ½ATP c is visible at low stimulation by IP 3 , i.e., ½IP 3 % 0:24 mM (Figures 2C and S3B).In those conditions, oscillations arise for ½AcCoA % 5 mM (Figure 2C), which correlates with a low ½ATP c (Figure S3B), and SERCA's activity is ATP-sensitive (Figure S2).Starting from physiological conditions (½AcCoA = 10 mM), ½ATP c decreases as ½AcCoA is decreased except at the onset of oscillations, where a cusp in ½ATP c can be observed (Figures 2E and S3D).At this critical point, ½ATP c becomes limiting for Ca 2+ uptake by SERCAs and mitochondrial exchanges take over despite a low DJ (Figure S2).Mitochondrial exchanges intensify as ½AcCoA is further decreased (Figure S2) and the larger Ca 2+ efflux from mitochondria (Figure S2) maintains CICR active, hence the decrease in oscillation period.The amplitude of the cusp is much smaller as compared to the cusp observed with GPCR stimulation because the activation of dehydrogenases by Ca 2+ plays here a minimal role due to the limited activity of the TCA cycle at low AcCoA level.As shown in Figure S2, the average reaction fluxes of the Ca 2+ -sensitive dehydrogenases (J IDH and J KGDH ) and OXPHOS (J Ox and J F1 ) decrease monotonically as ½AcCoA decreases.As explained in more detail in subsection Robustness of the efficiency-rescuing effect of Ca 2+ oscillations, the oscillatory dynamics plays a central role in the rescue of the ATP level because of the phase shift between ER and mitochondrial Ca 2+ sequestration mechanisms.
To assess the impact of mitochondrial uptake on DJ and mitochondrial metabolism, we inspected the bifurcation diagrams of the mitochondrial uptake flux (J UNI ), on the one hand, and of DJ, F1F0-ATPase flux (J F1 ) and respiration flux (J Ox ), on the other hand (Figures S1 and S2).For ½AcCoA close to 10 mM and upon increase of ½IP 3 , the elevation of J UNI is accompanied by a significant increase in J F1 , and J Ox , and a slight increase in average DJ, confirming that the activation of the mitochondrial dehydrogenases and the subsequent enhancement of OXPHOS overcomes the depolarization induced by the Ca 2+ entry (Figure S1).In contrast, a decrease in DJ, J F1 and J Ox is correlated to an increase of J UNI as ½AcCoA is reduced.However, the decrease in DJ cannot be attributed only to the larger J UNI , but also to the reduced TCA cycle flux (J IDH , J KGDH ) feeding the electron transport chain (J Ox ) that maintains the polarization of the mitochondrial membrane (Figure S2).In conclusion, the role of Ca 2+ on metabolism is beneficial upon GPCR stimulation but is potentially inhibitory in substrate-limited conditions.
To summarize, our data suggest that the cusp occurs when ATP c hydrolysis is counterbalanced by the increased ATP production by Ca 2+activated mitochondrial metabolism, which occurs when mitochondrial Ca 2+ uptake becomes significant.The Ca 2+ build-up activating mitochondrial exchanges and metabolism is either due to CICR (high ½AcCoA and ATP-insensitive SERCA's activity) or to insufficient Ca 2+ uptake by SERCA due to ATP c deficiency (low ½AcCoA and ATP-sensitive SERCA's activity).The role of oscillations in the cusp behavior of ATP c is also to be emphasized.For high ½AcCoA and ATP-insensitive SERCA's activity, the oscillatory dynamics allows for higher Ca 2+ levels in the cytosol, which can then be transmitted to mitochondria and sensed by dehydrogenases, without inducing cytotoxic effects.For low ½AcCoA and ATPsensitive SERCA's activity, oscillations enable a desynchronization between ER and mitochondrial Ca 2+ uptake, which intermittently boosts ATP c production while reducing ATP c hydrolysis (see subsection Robustness of the efficiency-rescuing effect of Ca 2+ oscillations).

The efficiency of mitochondrial metabolism displays a maximum in the regime of Ca 2+ spiking
The nonlinear ATP production observed for different ½IP 3 and ½AcCoA (Figures 2D, 2E, and S3B) suggests variations in the output work of mitochondria and, possibly, in the thermodynamic efficiency of their metabolism.As confirmed computationally, the output nonconservative work (w r1out ) displays a minimum (corresponding to maximal export of energy from mitochondria) that coincides with the maximal ½ATP c in the kinetic simulations (Figures 3A top vs. 3B bottom).In the extreme case where oscillations disappear for large stimulation by IP 3 , the efficiency drops and reaches a plateau (Figure 3B top).otably, the maximum in ½ATP c translates into a maximum in the efficiency of mitochondrial metabolism (Figures 3C and S4B).Such maxima are not systematically observed when ½AcCoA is varied at fixed ½IP 3 (Figures S3A and S3B).However, the increase in efficiency at the onset of the oscillatory regime is a robust feature that points to the stabilizing effect of Ca 2+ spikes on mitochondrial energetics.An increase over a similar range of mitochondrial Ca 2+ concentration can be induced in an IP 3 -independent way, for example by increasing the leakage of Ca 2+ from the ER (i.e., increase in V LEAK with a low ½IP 3 ).In that case, the Ca 2+ build-up is not accompanied by the emergence of oscillations and does not boost the associated mitochondrial efficiency (Figure S5), which highlights the special role of the oscillatory Ca 2+ dynamics on mitochondrial metabolism in contrast to steady-state regimes.
The relation between mitochondrial efficiency and dissipation is different in fast Ca 2+ spiking regimes triggered in substrate-limited or overstimulation conditions Like in other biological processes such as the migration of molecular motors along microtubules, kinetic proofreading or the regulation of circadian clocks, 70 the system's efficiency is maximal at intermediate levels of dissipation, corresponding to a limited range of ½IP 3 (Figure 3C).  4. Figure S3B illustrates the behavior of ½ATP c for an extended range of ½AcCoA and ½IP 3 .
Overstimulation of the signaling machinery by IP 3 is counterproductive since it only increases dissipation (Figure 3A bottom).By exploring the behavior of efficiency and EPR at large ½IP 3 , we observed a saturation effect (Figures 3A-3C) leading to limiting values for the efficiency (z 0:349) and total dissipation (z26600 JL À 1 s À 1 ).The dependency of efficiency on the total dissipation is highly nonlinear (Figure 3C).Around the onset of oscillations (and for a limited range of dissipation rates), a given dissipation rate can be associated to different efficiencies, in which case the highest efficiency is always reached for the highest Ca 2+ spiking frequency, while the lowest efficiency corresponds to the steady-state regime.Reversely, different dissipative regimes can yield the same efficiency.In that case, steady-state regimes display the lowest EPR while the fast-spiking regimes are the more dissipative regimes.We hypothesize that in such instances, the selection of the dissipative regime could be guided by constraints imposed by the global energy budget of the cell.
As ½AcCoA is increased, slow-spiking regimes are more dissipative than the fast-spiking ones (Figure S4A, bottom panel) and the efficiency increases almost linearly with dissipation in the oscillatory regime (Figure S4C), which contrasts with the overstimulation by IP 3 .A reason for this difference might be that upon overstimulation by IP 3 , ATP c hydrolysis overbalances the enhanced production of ATP by mitochondrial metabolism, which is not the case at low ½AcCoA.To summarize, fast spiking is less efficient than the slow spiking observed around the bifurcation point, and is also more dissipative when oscillations frequency intensifies due to perturbations in ½IP 3 .

Robustness of the efficiency-rescuing effect of Ca 2+ oscillations
Our proposed mechanism for the maximum in efficiency relies on the dependence of the SERCA flux on the hydrolysis of ATP c and on the resulting modulation of Ca 2+ sequestration mechanisms.If Ca 2+ homeostasis was ATP-independent, Ca 2+ would always exert a positive feedback on the TCA cycle flux and the efficiency of metabolism would increase monotonically with the Ca 2+ release from IP 3 Rs.We validated this hypothesis by performing simulations with a modified SERCA flux that is uncoupled from ATP c hydrolysis.The degradation of ATP c then relies exclusively on other ATP-consuming processes mimicking cellular activity (Hyd reaction in Figure 1A).As expected, the Ca 2+ -enhanced ATP production by mitochondria is not restrained upon more intense stimulation by IP 3 (Figure 4A).This uncoupling also disables the feedback of ATP production by mitochondria on Ca 2+ oscillations: instead of decreasing, the spike period is barely changed as ½AcCoA decreases (Figure 4B).
In the uncoupled case, J SERCA is not limited by the depletion of ATP c and both removal mechanisms proceed synchronously, although Ca 2+ uptake to the ER is predominant (Figure 4C, green dotted curve).While the mitochondrial Ca 2+ influx J UNI slightly increases with J SERCA , the Ca 2+ -dependent currents of the TCA cycle, J IDH and J KGDH (Figure 4C, purple dotted curves), are barely affected.Upon regulation by ATP c and in substrate-limited conditions, J SERCA varies over a more restricted range, is on average smaller and proceeds with a slight phase shift with respect to J UNI , which allows for a larger Ca 2+ influx in mitochondria and a more intense activation of the TCA cycle enzymes (Figure 4C solid curves).This mechanism reinforces the ''efficiency-rescuing'' role of mitochondrial buffering in substrate-limited conditions.
We also explored the robustness of our results against perturbations in the uptake rate of SERCA pumps of the original model.We mimicked the inhibition of SERCA pumps by decreasing the limiting rate V SERCA max .The efficiency-dissipation relation displays the same features as in the non-inhibited case (Figures 3C and S4C vs. Figures 3E and S4E, respectively).The results of an extended perturbation analysis (Table 6) confirming the robustness of our conclusions are detailed in subsection Sensitivity analysis of the STAR Methods.
Together, these results confirm that the cross-talk between Ca 2+ signaling and mitochondrial energy metabolism is a major mechanism underlying the maximum in efficiency arising in the spiking regime, even when the amplitude of this coupling is reduced due to the inhibition of SERCA pumps.

DISCUSSION
Here, we examined the impact of Ca 2+ signaling on the efficiency of mitochondrial metabolism by using tools from the nonequilibrium thermodynamics of CRNs on a detailed and experimentally calibrated kinetic model of the Ca 2+ -metabolism cross-talk.Our results highlight that, despite a usually higher dissipation rate compared to steady-state regimes, Ca 2+ oscillations can enhance the efficiency of mitochondrial metabolism.Increasing stimulation by IP 3 or inducing energy deficiency by decreasing the mitochondrial AcCoA level reduces the steadystate efficiency of metabolism, but at the onset of oscillations, the efficiency raises with a cusp-like transition and reaches a maximum of about 40% before decreasing again at higher stimulation/lower AcCoA level.This value is in the range of the efficiencies of the TCA cycle and OXPHOS (30% and 42%, respectively) estimated in the absence of regulation with a nonequilibrium thermodynamic approach. 16Moreover, slow-spiking is less dissipative than fast-spiking.Thus, we hypothesize that, for a given cell state, there exists an optimal stimulation level leading to slow-spiking/low-dissipation oscillations which maximize the efficiency of metabolism during signaling.For higher stimulation, the Ca 2+ signaling machinery then generates more dissipative regimes of gradually decreasing efficiency.
In the broader context of physical bioenergetics, energetic costs are usually assessed by evaluating the Gibbs free energy of reaction (D r G) dissipated or the equivalent number of ATP molecules produced/consumed along the processes of interest. 8,14,15However, such purely thermodynamic approaches do not account for reaction kinetics and thus cannot quantify the rates of free energy transduction and dissipation.Significant efforts have been made in the direction of adding thermodynamic constraints in flux balance analysis of metabolic networks. 71,72A few attempts have also been made to account for more complex kinetic effects such as enzyme saturation, leading to insights into the trade-offs between energy production and enzyme costs in glycolysis. 12,73,74Nevertheless, all these approaches rely on optimized nonequilibrium steady-states, which may not correspond to physiological conditions and cannot capture the energetic impact of time-dependent behaviors, such as the energy-rescuing effect of Ca 2+ oscillations quantified here.Our approach overcomes these limitations, based on the rigorous thermodynamic analysis of a curated dynamical model.Due to the modular structure of the model, our approach can be extended with additional pathways, such as glycolysis or one-carbon metabolism, with the aim to perform integrative modeling of cell metabolism.
Finally, a key element of the metabolic efficiency management in the presence of Ca 2+ regulation is the possibility for Ca 2+ oscillations, as an increase in the steady-state Ca 2+ level is not sufficient to boost mitochondrial efficiency.Ca 2+ oscillations are shaped by the interplay between SERCA and mitochondrial uptakes.Alterations in Ca 2+ removal mechanisms due to mutations, generation of reactive oxygen species or remodeling of channel and pump expression are ubiquitous in pathological states such as mitochondrial 75 and neurodegenerative diseases, 76,77 cancer 78 or diabetes, 79 and therapeutic strategies targeting Ca 2+ homeostasis and signaling have started to emerge. 80,81Some of these changes can be captured by perturbations in the kinetic parameters of the Ca 2+ fluxes, 82 which would make the use of our approach in the context of disease quite straightforward.Overall, our methodology thus paves the way for a more systematic characterization of the dynamical energetic impact of metabolism regulation, which could improve the current understanding of pathway selection mechanisms in health and disease.

Limitations of the study
A limitation of the study is related to the scarcity of kinetic data available per cell and organism.The computational model was carefully calibrated on astrocyte data but a wide range of kinetic parameters (Table 4) was collected from previous models based on different cell types and mammalian organisms.Ideally, the kinetic consistency of the model could be improved by using parameters characterized in the same cell line and experimental conditions, bearing in mind that in vitro or in vivo data may also differ. 83Another limitation is that the rate of ATP consumption by cellular activity other than SERCA's is a highly coarse-grained estimation, although widely used in Ca 2+ modeling.Ongoing research about the energetic costs in biology might provide more accurate values in the future.Finally, the experimental verification of our modeling predictions on the  2+ release by IP 3 Rs (that is, increasing ½IP 3 ) monotonically increases the efficiency in the uncoupled case, which strongly contrasts with the nonmonotonic dependency on ½IP 3 in the coupled case.While varying over different ranges, the period evolves according to the same trends in both cases.(B) Both systems display similar responses in their efficiency upon variations in ½AcCoA, but the oscillation period of the uncoupled system is decreasing as ½AcCoA increases.Empty and filled dots correspond respectively to steady-state and oscillatory regimes -note the use of linear colorbar schemes for the period.(C) Phase portraits of Ca 2+ -dependent TCA cycle currents (purple) and mitochondrial Ca 2+ uptake (green), namely J IDH , J KGDH and J UNI , vs. ER Ca 2+ uptake, namely J SERCA .Note that J IDH and J KGDH are indistinguishable.Symbols denote steady-state values.Triangles (½IP 3 = 0:24 mM) and dotted curves (½IP 3 = 0:32 mM) correspond to the uncoupled case, while circles (½IP 3 = 0:18 mM) and solid curves (½IP 3 = 0:26 mM) represent the coupled case.(A) ½AcCoA = 10 mM, (B) ½IP 3 = 0:30 mM for the uncoupled case and ½IP 3 = 0:20 mM for the coupled case, (C) ½AcCoA = 1 mM.k Hyd = 1:3310 À 1 mM s À1 for the uncoupled case.The other parameter values are the same as in Figures 2 and 3. Table 6.Parametric conditions for each realization of the sensitivity analysis thermodynamic efficiency is challenging, since it requires the dynamic measurement of fluxes and Gibbs free energies of reactions, independently of competing metabolic processes.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: The Ca 2+ uptake by SERCAs is balanced by the efflux of Ca 2+ from the ER through IP 3 receptors (IP 3 Rs) and passive leak through the membrane of the ER.The flux through IP 3 Rs is modeled by a function depending solely on the concentrations of IP 3 and cytosolic Ca 2+ , as in previous models. 26

Fluxes and evolution equations
Overall, we consider 32 species and 17 chemical reactions taking place in or between compartments corresponding to the mitochondrial matrix, the cytosol and the ER.The chemical equations describing these reactions can be found in Table 1.A current J k (given in  1 in the main text), meaning that the effect of chemical reactions is balanced by additional processes which are not described explicitly in the model.The rate equations for the other species are given in Equations 6, 7, 8, 9, 10, 11, 12, 13,  14, 15, 16, 17, 18, 19, 20, 21, 22, and 23.Stoichiometric and volumetric coefficients are included in Equations 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,  17, 18, 19, 20, 21, 22, and 23 to guarantee mass balance across cell compartments.Specifically, a = V ER =V c and d = V m =V c , where V c , V ER and V m are the volumes of the cytosol, of the ER and of mitochondria, respectively.Ca 2+ buffering in these compartments is accounted for by the coefficients f c , f ER and f m , which correspond to the fraction of free Ca 2+ in the compartment of interest.Finally, we mention that all fluxes are expressed in mM/s and concentration units are mM, except for Ca 2+ concentrations which are in mM, and coefficient g = 10 3 mM/mM is therefore introduced to ensure consistency in units.6) 7) 8) d t ½ATP m = À J ANT + J F1 + J SL (Equation 10) 11) 12) 13) 14) 15) 16) 17) 18) 19) 20) 21) 22) (Equation 23)

Model calibration
Except for the leading constants (V max , r and k f ), most parameter values of our model were taken from previous models.The calibration of our model thus consisted in finding a set of V max , r and k f generating data in agreement with experimental data according to the following criteria: 1. negative Gibbs free energy of reaction for all processes (thermodynamic consistency); 2. realistic concentrations, concentration ratios and mitochondrial membrane potential for physiological conditions (½AcCoA = 10 mM and ½IP 3 = 0:1 mM) based on experimental data from different mammalian cell types (Table 5); 3. reproducing key trends (in ATP:ADP and period of Ca 2+ oscillations) in response to GPCR stimulation and mitochondrial substrate availability, characteristic of the Ca 2+ -metabolism cross-talk, based on data from cultured astrocytes.
Thermodynamic consistency was satisfied for all bifurcation diagrams, with a notable exception for Gibbs free energy associated to SERCA in the uncoupled case.Indeed, since SERCA's activity is modeled as independent from ATP c hydrolysis in the uncoupled case and Ca 2+ physiologically more abundant in the ER than in the cytosol, a positive Gibbs free energy of reaction is to be expected.We emphasize that the uncoupled case is unphysiological.
The bidirectional feedback between Ca 2+ signaling and mitochondrial metabolism has recently been investigated in C8-D1A astrocytes 36 and our model was calibrated to reproduce the following trends: 1. decreasing the glucose concentration in the cell culture medium decreases the ATP:ADP ratio; 2. increasing the stimulation of GPCRs decreases the period of Ca 2+ oscillations; 3. decreasing the glucose concentration in the medium decreases the period of Ca 2+ oscillations; 4. increasing the stimulation of GPCRs increases the ATP:ADP ratio.
The impact of ½IP 3 and ½AcCoA on the period of Ca 2+ oscillations (Figure 2C of the manuscript) is to be compared with the experimental data of Figures 4. 13

Concepts of biothermodynamics Definitions
The entropy production rate (EPR) associated to a chemical reaction r is given by where T is the absolute temperature while J r and D r G r are the current and Gibbs free energy of reaction r, respectively.In nonequilibrium thermodynamics, À D r G r is the force driving the reaction while J r is the reaction flux resulting from this force.The equilibrium state is characterized by zero forces and hence zero fluxes.The Gibbs free energy of reaction r is defined by 91 where S r i is the net stoichiometric coefficient of species i in reaction r and m i is the chemical potential of species i.Under the hypothesis of local equilibrium, i.e., state variables such as temperature and pressure relax to equilibrium on a much faster timescale than chemical reactions, the expressions for the chemical potentials derived at equilibrium still hold locally out-of-equilibrium. 91The chemical potential m i is thus given by where m + i and a i denote the standard chemical potential and the activity of species i, respectively.The activity accounts for the interactions between chemical species present in solution and is related to the concentration by the coefficient of activity g i , which depends on the ionic strength, 92 I, such that a i = g i ½i c + , where ½i and c + are the concentration of species i and standard concentration, respectively.In ideal solutions, g i = 1.Standard conditions correspond to atmospheric pressure p + = 1 bar and molar concentrations c + = 1 M.
Standard chemical potentials are directly related to the standard Gibbs free energy of reaction D r G + r = i i;eq is the equilibrium constant of reaction r.Hence, D r G r can then be rewritten as 27) From a practical point of view, standard Gibbs free energies of reaction and activity coefficients are usually available in thermodynamic tables.As described in the next subsection, further adaptations can be done to describe more adequately the physiological environment in which biochemical reactions take place.

Physiological conditions
Cells are compartmentalized into specialized organelles whose composition can widely differ.For example, mitochondrial and cytosolic pH are 8 60 and 7.2, 59 respectively.A plethora of buffering mechanisms regulate their internal concentration of metallic ions and pH, and thereby ensures the maintenance of homeostasis.Some chemical species can exist in different forms, that is, bound to metallic cations or at different levels of protonation (for example, ''ATP'' can be ATP 4 À , HATP 3 À , MgATP 2 À , etc.), and their relative abundance depends on the internal environment of the organelle.For the sake of compactness, biochemical reactions are thus usually written in terms of pseudoisomers, that is, without explicitly mentioning the state of the species and without detailing the consumption or production of protons/metallic ions by the reaction (e.g.ATP + H 2 O#ADP + P i ). 92o describe biochemical reactions from a thermodynamic point of view, their associated standard Gibbs free energy can be rescaled to match the equilibrium corresponding the physiological pH and metallic ions concentrations, but also incorporates the activity coefficients corresponding to a physiological ionic strength (I = 0:120 M 28,92 ).The resulting D r G 0+ k is subsequently used to calculate the transformed Gibbs free energy of reaction 28) where ½j is the concentration of pseudoisomer j and S k j is the net stoichiometric coefficient of pseudoisomer j in reaction k.Complementary approaches 93,94 have led to the development of databases 29,95 from which we retrieved D r G 0+ k for different pH, ionic strength and metallic ion concentrations.
Electrogenic processes, such as Ca 2+ exchanges, the exchange of ATP 4 À and ADP 3 À via the antiporter and the transfer of protons accompanying oxidative phosphorylation, constitute notable exceptions where the charges of species need to be explicitly accounted for.Indeed, electrostatic interactions affect the Gibbs free energy and in that case, the right-hand side of Equation 28 must also comprise the term related to the electric potential in the compartment of interest.We thus distinguish the charged species fig from the pseudoisomers fjg.More precisely, D r G 0 k = 29) where z i is the charge of species i (for example, z = + 2 for Ca 2+ ), F is the Faraday constant and V rðiÞ is the electric potential in the compartment r where species i is considered.Overall, this leads to S k i z i FV rðiÞ : (Equation 30) When a charged species is exchanged between two compartments, the D r G 0 k associated to this transport process depends on the difference of potential between the two compartments.For example, if we consider the transport of Ca 2+ from cytosol to mitochondria, D r G 0 k = RT ln Although pH (and hence proton concentrations) in mitochondria and cytosol are assumed to be constant due to strong buffering, proton transfer across mitochondrial membrane still affects the membrane potential (at least, at the local scale that is considered in the present model), which in turns affect D r G 0 k .The expressions for the D r G 0 k of each process of the model can be found in Table 3.

É :
The former is the set of species involved only in the internal reactions, while the latter includes the controlled species and the species involved also in the external reactions.Thus, the rate equations for internal and exchanged species can be respectively written as where V Y ref is the volume of the compartment to which species Y belongs, and I Y is the exchange current either accounting for the external reactions (named flux control 42 ) or modeling additional processes that are responsible for the homeostasis of the controlled species (named concentration control 42 ).On the one hand ATP c and ADP c are involved in the external reactions SERCA and HYD and hence

Second law for mitochondrial metabolism
In general, the second law of thermodynamics for open CRNs can be written as 38,41 T 32) where G is the (semi-grand) Gibbs free energy of the system, while _ w nc and _ w driv , respectively referred to as the nonconservative work rate and the driving work rate, are related to the energetic cost of maintaining CRNs out of equilibrium via the exchange of species fYg.Since G is a state function, its time derivative vanishes at steady state as well as when averaged over one period in the case of an oscillatory regime.
In the following, we use a topological analysis 38,41 to derive the explicit expressions of the nonconservative work rate and the driving work rate for mitochondrial metabolism.
Remark.The rate Equations 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, and 23 are coarse-grained, namely, each reactive process represents a sequence of out-of-equilibrium elementary reactions involving intermediate species whose dynamical behavior is not described.Each of these elementary reactions might affect the energetics of the whole system.Nevertheless, under the assumption of the existence of a time scale separation between the evolution of the species accounted by the dynamical model and the coarse-grained intermediate species, our thermodynamic analysis characterizes the correct energetics of the whole system, as proven elsewhere. 40,43

Conservation laws and emergent cycles
For our model, the stoichiometric matrix S encoding the net stoichiometric coefficients of internal species X and exchanged species Y in the internal reactions k i reads The 13 (linearly-independent) left-null eigenvectors of S, encoded as rows of the matrix which therefore satisfies LS = 0, define the conservation laws.Indeed, for the every row l of L (labeled using chemical symbols in Equation 34for reasons that will be explained in the following section), the quantity L l = P ; are named emergent cycles.They define sequences of reactions that overall leave the abundances of the internal species unchanged, since S X c r1 = 0 and S X c r2 = 0 by definition, while interconverting the exchanged species as they undergo the effective reactions

Potential and force species
The exchanged species fYg can be classified as either potential species fY p g or as force species fY f g.The potential species are the largest subset of exchanged species such that the submatrix L b Yp of L for the broken conservation laws and the potential species can be inverted, i.e., ðL b Yp Þ À 1 exists.The force species are the remaining species fY f g = fYg\fY p g.As discussed in other references, 38,41  In Equation 33, the horizontal lines split S into (Equation 37) In Equation 34the horizontal and vertical lines split L into (Equation 38) with u and b for unbroken and broken conservation laws, respectively.Note that the existence of ðL b Yp Þ À 1 defines a representation of the broken conservation laws given by ðL b Yp Þ À 1 L b where every broken conservation law involves only one potential species.This physically means that each quantity L l corresponding to a broken conservation law stops being conserved once a specific potential species is exchanged.Furthermore, no new quantities L l stop being conserved when the force species are exchanged.For this reason, each conservation law l in Equation 34 can always be labeled using the chemical symbol of the potential species that, once exchanged, would make L l a nonconserved quantity.

Nonconservative work
The general expression of the nonconservative work rate is given by 39) where 40) is the vector of the nonconservative forces, I Y f is the vector collecting the exchange currents of the force species, and m Y f (resp.m Yp ) is the vector of the chemical/electrochemical potential of the force (resp.potential) species.In our model, there are only two force species, namely, ATP c and H + c .The two corresponding nonconservative forces are given by  42) while the exchange currents are 43) and 44)  45) and 46) We numerically compute the nonconservative forces using the latter expressions.The expression of the nonconservative work rate for our model then becomes 47)

Driving work
The driving work rate is in general given by the sum of two contributions, _ w driv = _ w ch driv + _ w in driv ; (Equation 48) namely, the chemical driving work rate 49) and the interaction driving work rate _ w in driv = V ½e G in ð½Z; ½eÞ$d t ½e; (Equation 50) Here, ½Z = ð½X; ½YÞ is a vector collecting the concentrations of both internal and exchanged species, ½e is a vector collecting the concentrations of other interacting species (i.e., Na + m , Ca 2+ m ) that are not interconverted by the internal reactions and that do not appear in S, V ½e is the gradient with respect to ½e and G in ð½Z; ½eÞ is the interaction Gibbs free energy, whose exact expression depends on the model used to describe interactions. 41n our model, we compute the driving work rate over one period t p in the oscillatory regime (as it vanishes at steady state).The specific expression of the average chemical driving work rate is given by In conclusion, the chemical driving work rate over one period vanishes: w ch driv = 0 : (Equation 52) We cannot determine the explicit expression of the interacting driving work rate since our model does not provide the interaction Gibbs free energy G in ð½Z;½eÞ.Thus, we compute the driving work over one period by calculating the total entropy production rate as the sum of the individual EPR of each internal reaction k i = {ANT, F1, OX, CS, ACO, IDH, KGDH, SL, SDH, FH, MDH}, 53) from which we subtract the nonconservative work over one period: 54)

Thermodynamic efficiency
In subsection 2, we determined the specific expressions of work rates in the second law of thermodynamics (Equation 32), accounting for the free energy exchanges between mitochondria and their surroundings that balance dissipation and maintain the mitochondrial metabolism out of equilibrium.The thermodynamic efficiency of mitochondrial metabolism is defined by identifying which of these terms play the role of free energy input and output.To do so, we further split the nonconservative work rate (Equation 47) by recognizing that the effective reaction (Equation 35) is given by the sum of 2 (mass balanced) reactions (Equation 56) The effective reaction 55 corresponds to the net production (output) of free energy (ATP c ) by mitochondrial metabolism, while the effective reactions 56 and 36 represent a free energy input generated by the TCA cycle, the electron transport chain and the F1F0-ATPase, respectively in its ATP synthase or hydrolyzing mode.The latter mode optimizes the proton driving force while the the first maximizes ATP synthesis.This allows us to write the nonconservative work rate (Equation 47) as (Equation 57) The thermodynamic efficiency at steady state is thus given by (Equation 58) In the oscillatory regime we also have to account for the free energy provided by the driving work and, therefore, the thermodynamic efficiency averaged over a period is (Equation 59) In both cases, the thermodynamic efficiency quantify the amount of energy released by the synthesis of ATP normalized by the amount of energy injected in the mitochondria.In the main text, we use the notation h to refer to the average efficiency at steady-state or in the oscillatory regime.

Sensitivity analysis
We investigated the robustness of the rescuing effect of Ca 2+ on the thermodynamic efficiency of mitochondrial metabolism.First, we added perturbations to all the leading constants (V max , k f , r) of the model.Random perturbations in the range of [-10, 10]% and selected from a uniform distribution are applied to each parameter of the model (Table 6).In Figures S6A and S6B, we show the resulting bifurcation diagrams of mitochondrial efficiency in dependence on ½IP 3 (for ½AcCoA = 10 mM) and ½AcCoA (for ½IP 3 = 0:20 mM).For each realization, the thermodynamic consistency criterion is satisfied, i.e. all Gibbs free energies of reactions are negative.In all simulations, an increase in efficiency at the onset of Ca 2+ oscillations is visible, which supports the robustness of the results against cellular variability.We also performed a complementary sensitivity analysis with a larger perturbation range ([-75,75]%) to assess the impact of variability of one leading constant at a time on steady-state regime (½IP 3 = 0:20 mM) and oscillatory regimes (½IP 3 = 0:50 mM) with ½AcCoA = 10 mM.The fold-change to reference conditions (Table 6) in key variables and the absolute change period are shown in Figures S6C and S6D, respectively.
Cellular respiration becomes endergonic (DG Ox > 0) upon moderate to large variations in k Hyd , r F1 (decrease by 20% or more) and r res (increase by 20% or more) or upon large variations in V SERCA max , V ANT max (decrease by 50% or more) and V IDH max (increase by 50% or more).The reaction catalyzed by fumarate hydratase also becomes endergonic (DG FH > 0 -data not shown) for moderate to large perturbations of r res (decrease by 20% or more) and large perturbations of V CS max (decrease by 75%) or r F1 (increase by 50% or more).Those perturbations indeed lead to an accumulation of malate, which affects the spontaneity of the reaction.We exclude these thermodynamically unfavorable conditions from the remaining of the discussion.The efficiency is mainly impacted by the two parameters related to OXPHOS (r res and r F1 ), SERCA-independent ATPc hydrolysis (k Hyd ), by Ca 2+ exchanges (V LEAK , V IP3R max , V UNI max , V NCX max and, to a lesser extent, V SERCA max ), nucleotide exchanges (V ANT max ) and some enzymes of the TCA cycle (V CS max and V IDH max ).Generally, Ca 2+ transfer to mitochondria enhances the efficiency while Ca 2+ release in the cytosol decreases the steady-state efficiency.Large perturbations in V IP3R max (+75%) or in V IDH max (-75%) induce oscillations, respectively via stimulation of GPCRs and by limitation of SERCA pumps by ATP c deficiency, respectively.Ca 2+ oscillations have a beneficial impact on the efficiency in the first case because the increased production of NADH by activated mitochondrial dehydrogenases results in an increased production of ATP c .In the second case, the global rate of the TCA cycle is limited by the rate of isocitrate dehydrogenase, which leads to reduced DJ, efficiency and ATP c level.
Similar observations hold in the oscillatory regime.DG Ox is even more susceptible to become positive upon changes in k Hyd , r F1 , r res , V ANT max , and V IDH max , and variations in V NCX max (decrease by 50% or more) and V CS max (increase by 50% or more) can also make cellular respiration endergonic, probably because increasing DJ affects the spontaneity of the reaction.Note that decreasing V IP3R max (by 75%) leads to the disappearance of oscillations and the associated efficiency is lower than for the oscillatory regimes generated by larger values of that parameter.The kinetics of mitochondrial exchange processes (V NCX max and V UNI max ) have a more pronounced impact on the efficiency of oscillatory regimes as compared to steady-state efficiency, with a beneficial impact of mitochondrial Ca 2+ accumulation (low V NCX max and high V UNI max ) on the efficiency.Overall, the robustness of the rescuing effect of Ca 2+ oscillations on the thermodynamic efficiency of mitochondria is supported by the sensitivity analysis.

Figure 2 .
Figure 2. Kinetic behavior of the system (A and B) Ca 2+ c and ATP c concentrations over time for ½AcCoA = 10 mM and (A) ½IP 3 = 0:30 mM or (B) ½IP 3 = 0:50 mM.(C) Effect of ½IP 3 and ½AcCoA on the oscillation period.(D and E) Average concentration of Ca 2+ c and ATP c as a function of (D) ½IP 3 for ½AcCoA = 10 mM or as a function of (E) ½AcCoA for ½IP 3 = 0:20 mM.Empty and filled dots represent steady-state and oscillatory regimes, respectively, and the boundaries of the shaded areas correspond to the minimum and maximum concentrations.Parameter values are given in Table4.FigureS3Billustrates the behavior of ½ATP c for an extended range of ½AcCoA and ½IP 3 .

Figure 3 .
Figure 3. Stimulating the Ca 2+ signaling machinery impacts the dissipation and efficiency of mitochondrial metabolism via the Ca 2+ -metabolism cross-talk (A) Nonconservative work contributions, driving work and dissipation for different ½IP 3 .The driving work represents less than 0.017% of the EPR.At high stimulation, oscillations disappear in the favor of a nonequilibrium steady-state regime.(B) Efficiency and ATP c concentration as a function of ½IP 3 .(C) Efficiency as a function of the total dissipation.(D and E) Plots corresponding to (B-C) for V SERCA max = 0:096 mM s À 1 .Empty and filled dots correspond to steady-state or period-averaged quantities, respectively.Unless specified otherwise, parameter values are the same as in Figure 2D.An analogous thermodynamic behavior is observed upon stimulation by AcCoA (Figure S4).

Figure 4 .
Figure 4. Efficiency of mitochondrial metabolism regulated by Ca 2+ signaling, Ca 2+ sequestration fluxes and Ca 2+ -dependent TCA cycle reaction fluxes, without and with the coupling of SERCA pumps to ATP c hydrolysis (A) Stimulating Ca 2+ release by IP 3 Rs (that is, increasing ½IP 3 ) monotonically increases the efficiency in the uncoupled case, which strongly contrasts with the nonmonotonic dependency on ½IP 3 in the coupled case.While varying over different ranges, the period evolves according to the same trends in both cases.(B) Both systems display similar responses in their efficiency upon variations in ½AcCoA, but the oscillation period of the uncoupled system is decreasing as ½AcCoA increases.Empty and filled dots correspond respectively to steady-state and oscillatory regimes -note the use of linear colorbar schemes for the period.(C) Phase portraits of Ca 2+ -dependent TCA cycle currents (purple) and mitochondrial Ca 2+ uptake (green), namely J IDH , J KGDH and J UNI , vs. ER Ca 2+ uptake, namely J SERCA .Note that J IDH and J KGDH are indistinguishable.Symbols denote steady-state values.Triangles (½IP 3 = 0:24 mM) and dotted curves (½IP 3 = 0:32 mM) correspond to the uncoupled case, while circles (½IP 3 = 0:18 mM) and solid curves (½IP 3 = 0:26 mM) represent the coupled case.(A) ½AcCoA = 10 mM, (B) ½IP 3 = 0:30 mM for the uncoupled case and ½IP 3 = 0:20 mM for the coupled case, (C) ½AcCoA = 1 mM.k Hyd = 1:3310 À 1 mM s À1 for the uncoupled case.The other parameter values are the same as in Figures2 and 3.
(A) and 4.10 (B-C) of Moein, 36 respectively.The evolution of the ATP:ADP ratio in dependence of ½IP 3 and ½AcCoA (Figures S1 -before overstimulation -and S2) can be compared to Figures 4.

YY
½Y would be a conserved quantity if mitochondria were closed systems, namely, if I Y = 0 c Y. When I Y s0, only 3 out of the 13 conservation laws corresponding to the last three rows of L in Equation 34 involve exclusively internal species and their corresponding quantities L l = P X L l X ½X are still conserved.These conservation laws are said to be unbroken.The other conservation laws correspond to quantities L l = P ½Y that are not conserved anymore and are, therefore, named broken conservation laws.The 2 (linearly-independent) right-null eigenvectors of S X (i.e., the (sub)stoichiometric matrix for the internal species), ANT FI OX CS ACO IDH KGDH SL SDH FH MDH

Table 2 .
Fluxes of the system

iScience Article Table 4 .
Reference parameter values 33Forward rate constant of ACO 12.5 s À1 Cortassa et al.33(Continued on next page)

Table 4 .
Continued (Continued on next page) iScience Article

Table 5 .
Experimental values of concentrations, concentration ratios, mitochondrial membrane potential and oscillation period used for the calibration of the model liver tissue (rat) Wilhelm and Hirrlinger 89 , Williamson et al.

TABLE
d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d METHOD DETAILS B Kinetic model B Concepts of biothermodynamics B Nonequilibrium thermodynamic analysis

Table 2 )
is associated to each chemical reaction k and is normalized with respect to the volume of the corresponding cell compartment (cytosol, ER or mitochondria).The concentration of the controlled species {Pi c , Pi m , Na + c , Na + m , IP 3 , H + c , H + m , O 2 , H 2 O c , H 2 O m , AcCoA, CoA, CO 2 , CoQH 2 , CoQ} is fixed in time (see also Figure Mitochondria can be considered as engines converting ADP c into ATP c via 11 so-called internal reactions Y ˛ÈATP c ; ADP c ; Pi m ; H + c ; H + m ; O 2 ; H 2 O m ; AcCoA; CoA; CO 2 ; CoQ; CoQH 2 J HYD Á ; on the other hand, the other exchanged species (i.e., Pi m , H + m , H + m , O 2 , H 2 O m , AcCoA, CoA, CO 2 , CoQH 2 , CoQ) are controlled species and hence I Y = À P whose stoichiometry is defined by S Y c r1 and S Y c r2 , respectively, with S Y the (sub)stoichiometric matrix for the exchanged species.
˛ÈADP c ; Pi m ; H + m ; O 2 ; H 2 O m ; AcCoA; CoA; CO 2 ; CoQ; CoQH 2 this partitioning is not unique, but different choices do not change the conclusions of the thermodynamic analysis.
Equation 43 is obtained by writing the rate equation for ATP c (Equation9) according to Equation 31b: 1 d d t ½ATP c = J ANT + I ATPc .Equation 44 is obtained by recognising that H + c is a controlled species, i.e., d t ½H + c = 0, and using again Equation 31b, we can write 1 d d t ½H + c = 0 = 10 J OX À 3 J F1 + I H + c .Notice that the nonconservative forces in Equations 41 and 42 correspond exactly to D r G 0 r1 and D r G 0 r2 of the effective reactions in Equations 35 and 36.They can therefore be rewritten as 33 À d t m CoQH 2 L CoQH 2 33 À d t m CoQ L Notice that the terms corresponding to uncharged controlled species (i.e., Pi m , O 2 , H 2 O m , AcCoA, CoA, CO 2 , CoQH 2 , CoQ) vanish since their chemical potential is constant over time.For the charged controlled species H + m , the quantity L H + m is still conserved since ½H + m and ½H + c are constant.Hence, 1 ADP c + Pi m #