Research article Special Issues

A simple model of cardiac mitochondrial respiration with experimental validation

  • Cardiac mitochondria are intracellular organelles that play an important role in energy metabolism and cellular calcium regulation. In particular, they influence the excitation-contraction cycle of the heart cell. A large number of mathematical models have been proposed to better understand the mitochondrial dynamics, but they generally show a high level of complexity, and their parameters are very hard to fit to experimental data.

    We derived a model based on historical free energy-transduction principles, and results from the literature. We proposed simple expressions that allow to reduce the number of parameters to a minimum with respect to the mitochondrial behavior of interest for us. The resulting model has thirty-two parameters, which are reduced to twenty-three after a global sensitivity analysis of its expressions based on Sobol indices.

    We calibrated our model to experimental data that consists of measurements of mitochondrial respiration rates controlled by external ADP additions. A sensitivity analysis of the respiration rates showed that only seven parameters can be identified using these observations. We calibrated them using a genetic algorithm, with five experimental data sets. At last, we used the calibration results to verify the ability of the model to accurately predict the values of a sixth dataset.

    Results show that our model is able to reproduce both respiration rates of mitochondria and transitions between those states, with very low variability of the parameters between each experiment. The same methodology may apply to recover all the parameters of the model, if corresponding experimental data were available.

    Citation: Bachar Tarraf, Emmanuel Suraniti, Camille Colin, Stéphane Arbault, Philippe Diolez, Michael Leguèbe, Yves Coudière. A simple model of cardiac mitochondrial respiration with experimental validation[J]. Mathematical Biosciences and Engineering, 2021, 18(5): 5758-5789. doi: 10.3934/mbe.2021291

    Related Papers:

    [1] Sebastian Builes, Jhoana P. Romero-Leiton, Leon A. Valencia . Deterministic, stochastic and fractional mathematical approaches applied to AMR. Mathematical Biosciences and Engineering, 2025, 22(2): 389-414. doi: 10.3934/mbe.2025015
    [2] Grazia Policastro, Vincenzo Luongo, Luigi Frunzo, Nick Cogan, Massimiliano Fabbricino . A mechanistic mathematical model for photo fermentative hydrogen and polyhydroxybutyrate production. Mathematical Biosciences and Engineering, 2023, 20(4): 7407-7428. doi: 10.3934/mbe.2023321
    [3] Floriane Lignet, Vincent Calvez, Emmanuel Grenier, Benjamin Ribba . A structural model of the VEGF signalling pathway: Emergence of robustness and redundancy properties. Mathematical Biosciences and Engineering, 2013, 10(1): 167-184. doi: 10.3934/mbe.2013.10.167
    [4] Daniele Bernardo Panaro, Andrea Trucchia, Vincenzo Luongo, Maria Rosaria Mattei, Luigi Frunzo . Global sensitivity analysis and uncertainty quantification for a mathematical model of dry anaerobic digestion in plug-flow reactors. Mathematical Biosciences and Engineering, 2024, 21(9): 7139-7164. doi: 10.3934/mbe.2024316
    [5] Shuai Miao, Lijun Wang, Siyu Guan, Tianshu Gu, Hualing Wang, Wenfeng Shangguan, Weiding Wang, Yu Liu, Xue Liang . Integrated whole transcriptome analysis for the crucial regulators and functional pathways related to cardiac fibrosis in rats. Mathematical Biosciences and Engineering, 2023, 20(3): 5413-5429. doi: 10.3934/mbe.2023250
    [6] Martina Bukač, Sunčica Čanić . Longitudinal displacement in viscoelastic arteries:A novel fluid-structure interaction computational model, and experimental validation. Mathematical Biosciences and Engineering, 2013, 10(2): 295-318. doi: 10.3934/mbe.2013.10.295
    [7] Maria Pia Saccomani, Karl Thomaseth . Calculating all multiple parameter solutions of ODE models to avoid biological misinterpretations. Mathematical Biosciences and Engineering, 2019, 16(6): 6438-6453. doi: 10.3934/mbe.2019322
    [8] Yu Ichida, Yukihiko Nakata . Global dynamics of a simple model for wild and sterile mosquitoes. Mathematical Biosciences and Engineering, 2024, 21(9): 7016-7039. doi: 10.3934/mbe.2024308
    [9] Krzysztof Fujarewicz, Marek Kimmel, Andrzej Swierniak . On Fitting Of Mathematical Models Of Cell Signaling Pathways Using Adjoint Systems. Mathematical Biosciences and Engineering, 2005, 2(3): 527-534. doi: 10.3934/mbe.2005.2.527
    [10] Abdon Atangana, Jyoti Mishra . Analysis of nonlinear ordinary differential equations with the generalized Mittag-Leffler kernel. Mathematical Biosciences and Engineering, 2023, 20(11): 19763-19780. doi: 10.3934/mbe.2023875
  • Cardiac mitochondria are intracellular organelles that play an important role in energy metabolism and cellular calcium regulation. In particular, they influence the excitation-contraction cycle of the heart cell. A large number of mathematical models have been proposed to better understand the mitochondrial dynamics, but they generally show a high level of complexity, and their parameters are very hard to fit to experimental data.

    We derived a model based on historical free energy-transduction principles, and results from the literature. We proposed simple expressions that allow to reduce the number of parameters to a minimum with respect to the mitochondrial behavior of interest for us. The resulting model has thirty-two parameters, which are reduced to twenty-three after a global sensitivity analysis of its expressions based on Sobol indices.

    We calibrated our model to experimental data that consists of measurements of mitochondrial respiration rates controlled by external ADP additions. A sensitivity analysis of the respiration rates showed that only seven parameters can be identified using these observations. We calibrated them using a genetic algorithm, with five experimental data sets. At last, we used the calibration results to verify the ability of the model to accurately predict the values of a sixth dataset.

    Results show that our model is able to reproduce both respiration rates of mitochondria and transitions between those states, with very low variability of the parameters between each experiment. The same methodology may apply to recover all the parameters of the model, if corresponding experimental data were available.



    Mitochondria are intracellular organelles that have several roles, including production of ATP which provides energy to cell mechanisms. They also participate in the intracellular calcium (Ca2+) cycle. Ca2+ controls the contraction of myocytes as well as mitochondrial energy production [1]. In particular, a dysfunction in mitochondrial Ca2+ handling may very likely be associated to cardiac rhythm pathologies [2]. The exact role of the mitochondria in regulating the intracellular Ca2+ concentration is not fully understood.

    The different mitochondrial functions are coupled one to another, which makes their activity complex. It is characterized by many ionic flux, which can be either exchanges of components between the mitochondria and its environment, or reactions inside the mitochondrial matrix.

    To better understand this activity, a collection of mathematical models have been proposed, which was recently reviewed by Takeuchi and Matsuoka [3]. Most of these models are based on the pioneering work of Magnus and Keizer [4] and generally show a high level of complexity. These models are generally developed following a cumulative approach, adding new mechanisms and biological reactions to previous models, which leads to dozens of state variables and hundreds of parameters. This large number of equations and parameters is not compatible with calibration techniques required for validation against experimental data. Indeed, using too large sets of parameters, which all have uncertainties, reduces the identifiability of the parameters, and often results in cases of overfitting. Several groups, including Bertram et al. [5], proposed a different modeling approach by simplifying the mathematical expressions and reducing their number of parameters, still capturing the mitochondrial dynamics.

    In this paper, we propose a new model that has a number of parameters small enough to be fitted with the available experimental data. We followed the approach initiated by Bertram et al. and applied it to cardiac mitochondria. We also describe the methodology used to actually identify the relevant parameters of the model, which we successfully applied to mitochondrial respiration data.

    The derivation of a simple model is motivated by the deciphering of the interplay between intracellular calcium dynamics and mitochondria and cardiac electrophysiology. To this aim, we need a model simple to be coupled to an existing cardiac ionic model, such as the Mitchell-Schaeffer [6] or Ten Tuscher-Panvilov models [7], which do not account explicitely for the mitochondrial activity. In a cardiac simulation tool at the tissue scale, information on the ionic state must be stored at each degree of freedom, i.e., several hundred thousands of times. Using a complex mitochondrial model with dozens of state variables would considerably increase the computational burden of cardiac simulations. Therefore, there is interest in developing models that make the balance between simplicity and fidelity to the biological phenomena.

    In order to derive our model, we revisited thoroughly each flux expression starting from the original work of Hill [8], which has been the base of the work of Magnus & Keizer, and most of the subsequent mitochondrial models. We simplified the final expressions by surface fitting to simpler analytical expressions (sec. 2.1 and 2.2), as done by Bertram et al. [5]. Additionally, we adopted an original modeling approach by writing the model in terms of thermodynamics variables instead of concentrations.

    The final simplified model consists of a system of 6 ordinary differential equations (ODE) for 6 state variables: intra- and extra-mitochondrial calcium and ADP-based phosphate energy potential, transmembrane potential difference and the NADH redox potential. These variables are coupled through the right hand sides of the ODEs that are functions representing the internal and external ionic flux. These functions describe seven main features of mitochondrial oxidative phosphorylation: substrate metabolization, oxygen consumption by the respiratory chain, ATP synthesis trough phosphorylation, ATP/ADP translocation, Ca2+ regulation through the calcium uniporter and the Na+/Ca2+ exchanger and finally proton leakage.

    Even though our model is relatively simple, it involves 32 parameters overall. Hence, it remains very challenging to fit all these parameters to experimental data. Consequently, we proposed a calibration methodology that starts with a global sensitivity analysis to eliminate parameters with small influence on the flux functions (sec. 4.2). In a second step, we performed the sensitivity analysis on the output of the model that can be directly linked to experimental data (sec. 4.3), which consisted in our case of oxygen consumption rates by isolated mitochondria (sec. 3.1). Finally (sec. 5), we used the results of this sensitivity analysis to calibrate the remaining parameters to experimental respiratory rates, using a genetic optimization algorithm.

    We obtained parameters that are able to mimic steady states of the mitochondrial respiration. This set of parameters does not vary by more than 2% across the five different experimental datasets that were available, attesting the relevance of our model.

    The model proposed in this article, and more importantly the modeling and parameter identification approaches are a crucial intermediate step towards a model that describe the regulation of intracellular calcium by cardiac mitochondria, with calibrated parameters.

    We based our model on the work of Magnus and Keizer [4] (MK model) that discussed initially the mitochondrial activity in the pancreatic β-cell. We chose to keep the seven main mechanisms that were considered by the MK model: substrate metabolization, respiration or oxygen consumption, ADP phosphorylation, adenine nucleotide translocation, Na+/Ca2+ exchange and Ca2+ uniporter and proton leakage (see Figure 1). These mechanisms are modeled by ionic flux expressions that are functions of the state variables and parameters of the model. The expressions of the respiration and oxidative phosphorylation flux were adapted by Magnus and Keizer from a previously published proton pump model of Pietrobon and Caplan [9], who used the Hill-diagram method [8] to derive their equations.

    Figure 1.  The seven reactions and flux described by our model, indicated by boxes. State variables are indicated by bold font.

    The MK model was adopted later on by several authors, including Nguyen et al. in [10,11], Bertram et al. in [5] and Cortassa et al. in [12,13,14], while adding more mechanics and biological pathways to model the mitochondria in a cardiac cell context.

    In those models, the input of the respiratory chain flux are written in terms of concentration of mitochondrial NADH ([NADH]m), and are driven by the oxydation-reduction equation:

    NADH+H++12O2NAD++H2O. (2.1)

    ADP phosphorylation input flux are written in terms of mitochondrial ADP concentration ([ADP]m).

    We also kept the conservation of the total concentration of nicotinamide adenine dinucleotide (NAD) and adenine nucletotide in the mitochondrial matrix :

    Ntot= [NADH]m+[NAD+]m, (2.2)
    Atot,m= [ATP]m+[ADP]m. (2.3)

    The mitochondrial transmembrane potential difference Δψ and the calcium concentration [Ca2+]m are also variables of the model, so that the dynamics of the mitochondrial activity can be modeled using the following ODE system:

    {ddt[NADH]m= (JPDHJO),ddt[ADP]m= (JANTJF1F0),ddtΔψ= (JH,respJH,ATPJANTJleak2Juni)/Cm,ddt[Ca2+]m= fm(JuniJNaCa), (2.4)

    where Cm and fm are the mitochondrial membrane capacitance and the fraction of matrix free calcium respectively. The third equation of the system 2.4 comes from applying Kirchhoff's current law to mitochondrial membranes. The capacitance Cm, which is usually in units of Farad, is expressed in nmol mV−1 mgprot−1 to match the unit of the ionic fluxes through the membrane. The values of Cm and fm are given in Table 1. We detail in the next section the expressions of each flux Ji.

    Table 1.  Parameters associated to the ODE system (2.18) and Eqs (2.19), (2.20), (2.21).
    Parameter Description Value Unit Reference
    fm Fraction of free matrix calcium concentration 3 × 104 dimensionless [4]
    Cm Membrane capacitance 1.450 × 103 nmol mV−1 mgprot−1 [4]
    fc Fraction of free cytosolic calcium concentration 1 × 102 dimensionless [15]
    γ Volume ratio × flux unit conversion 0.0105 dimensionless Calibrated (sec. 5)
    Ntot Total NADH/NAD+ concentration 8 nmol mgprot−1 [4]
    Atot,m Total adenine nucleotide conc. (mitochondrial) 12 nmol mgprot−1 [4]
    Atot,c Total adenine nucleotide conc. (cytosol) 2 mM [15]

     | Show Table
    DownLoad: CSV

    In this section, we explain the biological role of each of the seven considered mechanisms, and we detail the derivation of their associated flux functions. Expressions from the current literature were derived by successive contributions, starting from first mathematical statements of the various biochemical processes involved. Consequently, these processes and their current expressions are difficult to relate. Hence, we started our derivation process from the original papers.

    In the cytoplasm, glucose is metabolized to pyruvate, which enters the mitochondria through a carrier. Then, pyruvate enters the Krebs cycle through the help of the pyruvate dehydrogenase (PDH) enzyme, which leads subsequently to NADH production. PDH has an active form (PDHa) that becomes completely inactivated when phosphorylated (PDH-P). The conversion between these two forms is regulated by the kinase and phosphatase enzymes. Magnus and Keizer [15] discussed the dependence of the fraction of PDHa on mitochondrial calcium, as it is believed to stimulate the phosphatase reaction that produces PDHa. They proposed an expression of the fraction of PDHa:

    fPDHa=11+u2[1+u1(1+[Ca2+]mKCa)2], (2.5)

    where u1 is the dependence of PDH-P dephosphorylation on Mg2+, u2 is the ratio between the maximal rate of kinase and the maximal rate of phosphatase, and KCa is the Ca2+-dependant affinity of the PDH-P phosphatase. Bertram et al. [5] then used this expression, with the addition of an explicit dependence on the ratio [NADH]m/[NAD+]m, to model the mitochondrial NADH production considering glycolysis, pyruvate dehydrogenase and its process through the Krebs cycle as a single compartment. This rate of production is denoted by JPDH and has the following expression:

    JPDH=Jmax×11+u2[1+u1(1+[Ca2+]mKCa)2]×1[NADH]m[NAD+]m+KPDHnad, (2.6)

    which parameters Jmax, u1, u2, KCa, and KPDHnad are given in appendix in Table 4. The value of the parameter Jmax depends on the external substrate: glutamate, malate or both.

    Table 2.  Parameters associated to the rates and flux of our model (section 2.2.2) computed from fitting to the original expressions (section 2.1).
    Parameter Value Unit
    p0 3.222×102 nmol mgprot−1 min−1
    p1 1.093×10-1 mV−1
    p2 1.330×102 mV−1
    p3 1.507×10-2 mV−1
    p4 1.119×103 mV
    p5 1.082×103 mV
    p6 3.530×101 nmol mgprot−1 min−1
    p7 4.038×10-1 dimensionless
    p8 2.300×104 mgprot nmol−1
    p9 1 dimensionless
    p10 9.791×102 nmol mgprot−1 min−1
    p11 1.171×10-2 mV
    p12 2.588×102 mV
    p13 -5.596×10-2 mV−1
    p14 1.761×102 mV
    p15 -3.613×10-5 nmol mgprot−1 min−1
    p16 1×103 nmol mgprot−1 min−1
    p17 8×10-1 dimensionless
    p18 6.156×105 dimensionless
    p19 5×10-1 dimensionless
    p20 2.740×10-3 nmol mgprot−1 min−1
    p21 5.245×10-2 mV−1
    p22 2.276×101 min−1
    p23 3×10-3 nmol mgprot−1
    p24 1.006 min−1 mV−1
    p25 8.4×101 mV
    p26 6.25×10-1 nmol mgprot−1

     | Show Table
    DownLoad: CSV
    Table 3.  Parameters of the model provided by the genetic algorithm, minimizing the cost function (3.2) on 5 different sets of experimental data.
    Parameter Unit Average Standard deviation
    p2 mV 1.314×102 1.15
    p5 mV 1.035×103 7.32
    p11 mV 1.194×10-2 1.76×10-4
    p12 mV 2.823×102 1.52
    p14 mV 1.7×102 2.58
    p21 mV−1 5.502×10-2 1.38×10-3
    γ dimensionless 1.057×10-2 1.007×10-3

     | Show Table
    DownLoad: CSV

    During the oxydation of the electron carrier NADH, two electrons are transported through the electron transport chain (ETC). This transport is coupled with proton pumping across the inner mitochondrial membrane, which generates a proton gradient. The gradient is then used to power the oxidative phosphorylation through the ATP synthase complex.

    We considered a model of a respiration driven proton pump that is adapted from the six-state proton pump model of Pietrobon & Caplan [9]. This model is constructed using the more general theory of coupled flows and membrane transport, developped by Hill [8], that is extended from the Altman and King diagrams method to be later on known as Altman-King-Hill diagram method.

    With this method, we obtained the following expression for the rates of proton ejection JH,resp and oxygen consumption JO:

    JH,resp=12JO=6ρrespJcycle, (2.7)

    where Jcycle is a highly complex expression with 13 parameters. The coefficients 12 and 6 account for the stoichiometry of reaction (2.1) between electrons, protons and oxygen atoms. The parameter ρresp=0.4 nmol mgprot−1 is the concentration of proton pumps introduced by Magnus & Keizer [4]. The derivation of Jcycle and its parameters are stated in section B.1 and in Table 5.

    The proton electrochemical gradient created by the respiratory chain is used by the F1F0-ATPase complex to synthesize ATP by phosphorylation of ADP. Using the fact that ATP hydrolysis is coupled to proton ejection, Pietrobon and Caplan [9] modeled the mitochondrial F1F0-ATPase by a six-states proton pump. We built the ADP phosphorylation flux from this model, adopting the same modeling approach as for the respiratory chain (see appendix B.2).

    Similar calculations as in section 2.1.2 were repeated to obtain the expressions of the rate of ATP hydrolysis and its associated proton pumping:

    JH,ATP=3JF1F0=3ρF1Jcycle, (2.8)

    where the coefficient 3 corresponds to the stoichiometry between ATP hydrolysis and its associated proton pump. The expression of Jcycle is given in appendix B.2, and its parameters, as well as ρF1, in Table 6.

    The ANT transports the mitochondrial ATP from the matrix to the cytoplasm, and the cytoplasmic ADP to the matrix. Magnus & Keizer wrote down an expression for the rate of this translocator:

    JANT=Jmax,ANT10.05[ATP]c×0.45×0.8[ADP]m0.45[ADP]c×0.05[ATP]mexp(FRTΔψ)(1+0.05[ATP]c0.45[ADP]cexp(fFRTΔψ))(1+0.45×0.8[ADP]m0.05[ATP]m). (2.9)

    This rate depends primarily on the membrane potential Δψ, as well as on the concentration of mitochondrial and cytoplasmic ATP and ADP respectively. In the expression (2.9), the coefficients leading the adenine nucleotide concentrations are the fractions of unbound ATP and ADP, existing in ionized form, in the matrix and cytoplasm compartments. Other parameters are listed in Table 7.

    Calcium is exchanged between mitochondria and the cytoplasm through membrane channels and pores such as the Ca2+ uniporter, the Na+/Ca2+ exchanger and the mitochondrial permeability transition pore (mPTP). The exact role of the mPTP remains currently not fully understood. In our model we considered only the Ca2+ uniporter and the Na+/Ca2+ exchanger as calcium regulators.

    We assume that the Na+/Ca2+ exchanger is electro-neutral (two Na+ for one Ca2+) [16,17,18]. Thus the flux does not depend on the membrane potential of the mitochondria, and we modeled the exchanger by a two-substrates Michaelis-Menten kinetic model (see Table 8 for parameters):

    JNaCa=Jmax,NaCa1(1+(KNa[Na]c)2)(1+KCa[Ca2+]m) (2.10)

    Several studies proposed similar expressions for the Na+/Ca2+ exchanger using different modeling approaches [4,19].

    For the calcium uniporter, Magnus and Keizer modeled the ionic flux using the Nernst-Planck equation, and then added the calcium dependence to the equation through an allosteric model [4]. The rate of the calcium internalized by the Ca2+ uniporter Juni depends on the transmembrane voltage Δψ and on the cytoplasmic calcium concentration [Ca2+]c, as follows:

    Juni=Jmax,uni[Ca2+]cKtrans(1+[Ca2+]cKtrans)3(1+[Ca2+]cKtrans)4+L(1+[Ca2+]cKact)na.2FRT(ΔψΔψ)1exp(2FRT(ΔψΔψ)). (2.11)

    The parameters of the calcium uniporter flux are given in Table 9.

    Several mitochondrial models [4,12] consider the external leakage of protons across the mitochondrial membrane to be a linear function of the transmembrane potential or the proton-motive force. However it has been shown [20,21] that this leakage increases exponentially with the proton-motive force.We chose such an exponential function to model the protons leakage, as described in section 2.2.2.

    The model discussed in section 2.1 is rather complex in terms of number of parameters and flux expressions, which makes the fitting to experimental data mathematically impossible.

    In this section we propose a model with simplified flux functions, which involve a reduced number of parameters, but still capture the same dynamics as the complete model. Indeed, each flux expression is a function of 1 or 2 model variables, except the ANT flux, function of 3 variables. So, it was possible derive simpler expression for each of these flux by surface fitting, using as few parameters as possible.

    We also chose to describe the mitochondrial behavior in terms of thermodynamic variables (actually related to an electrochemical force induced by chemical gradients of concentrations), instead of the concentrations, since they are more informative to biologists than absolute concentrations.

    Our model uses the following state variables (in bold face within this section).

    NADH redox potential. It is a logarithmic function of the redox couple NADH and NAD+, and is expressed in units of mV as follows:

    ΔEresp=RTFln(Kresp[NADH]m[NAD+]m), (2.12)

    where Kresp=1.35×1018 is the redox equilibrium constant of the reaction (2.1) [4].

    Using the conservation relation (2.2), the NADH redox potential reads:

    ΔEresp([NADH]m)=RTFln(Kresp[NADH]mNtot[NADH]m). (2.13)

    Gibbs free energy. It describes the energy balance between mitochondrial ATP and ADP. It is expressed in units of mV as follows:

    ΔGp,m=RTFln(KF1[ATP]m[ADP]m[Pi]m), (2.14)

    where KF1=1.71×106 mM is the equilibrium constant for the ATP hydrolysis reaction [4], and the inorganic phosphate concentration [Pi]m is considered as a fixed parameter.

    The conservation relation (2.3) gives:

    ΔGp,m([ADP]m)=RTFln(KF1Atot,m[ADP]m[ADP]m[Pi]m). (2.15)

    We kept the variables for the transmembrane potential Δψ and the mitochondrial calcium concentration [Ca2+]m since they are of main interest for us.

    In addition to these four state variables, we also considered the cytoplasmic calcium concentration ([Ca2+]c), and the cytoplasmic ADP concentration ([ADP]c) or equivalently the cytoplasmic Gibbs free energy change due to adenosine variations:

    ΔGp,c=RTFln(KF1,c[ATP]c[ADP]c[Pi]c). (2.16)

    These two variables were added because they are the controlable inputs of the experiments. Since we add the variable [ADP]c to the model, we assumed the following conservation relation:

    Atot,c=[ATP]c+[ADP]c. (2.17)

    New ODE system. The new ODE system that models the cardiac mitochondrial activity was then written as follows:

    {ddtΔEresp= φ1(ΔEresp)(JsubJO),ddtΔGp,m= φ2(ΔGp,m)(JF1F0JANT),ddtΔGp,c= φ3(ΔGp,c)(γJANTddtJADP,ext(t)),CmddtΔψ= 12JO3JF1F0JANTJleak2Juni,ddt[Ca2+]m= fm(JuniJNaCa),ddt[Ca2+]c= fc˜γ(JNaCaJuni)+ddt[Ca2+]ext(t) (2.18)

    where

    φ1(ΔEresp)=12RTF(K2resp+exp(2FRTΔEresp))2K2respNtotexp(2FRTΔEresp), (2.19)
    φ2(ΔGp,m)=RTF(1+[Pi]mK1F1exp(FRTΔGp,m))2Atot,m[Pi]mK1F1exp(FRTΔGp,m), (2.20)

    and

    φ3(ΔGp,c)=RTF(1+[Pi]cK1F1,cexp(FRTΔGp,c))2Atot,c[Pi]cK1F1,cexp(FRTΔGp,c), (2.21)

    are functions that come from the derivation chain rule when mapping concentrations to thermodynamical variables. These functions have units of mV nmol−1 mgprot. The functions JADP,ext and [Ca2+]ext are source terms for external ADP and Ca2+ respectively, allowing us to mimic various experimental conditions.

    The scalar parameter γ regroups two factors: the mitochondria over cytoplasm volume ratio, and a unit conversion factor from nmol to mM. The parameter ˜γ is simply 103×γ to account for the unit of [Ca2+]c in µM. The values of the parameters used in Eqs (2.18), (2.19), (2.20) and (2.21) are listed in Table 1.

    Expressions Ji derived previously in section 2.1 show a high level of complexity in terms of number of parameters and non-linearity. To simplify some of these expressions, we first plotted them as function of their state variables. The shapes of these graphs lead us to propose simpler expressions, with fewer parameters, that can reproduce the same dependencies. We then calibrated the parameters of the simplified expressions so as they fit the original ones. Bertram et al. had a similar approach [5]. However, we used surface fitting instead of several unidimensional fits, and performed the identification using a least-squares method on the whole domain spanned by the state variables. The simplified expressions are stated below and the calibrated set of parameters is listed in Table 2.

    Oxygen consumption rate. The oxygen consumption rate JO defined in Eq (2.7) has a sigmoidal shape when plotted as a function of its variables (ΔEresp,Δψ). The proposed simplified expression for this rate is:

    ˜JO(ΔEresp,Δψ)=σ(ΔEresp,0,σ(Δψ,p0,0,p1,p2),p3,σ(Δψ,p4,p5,p1,p2)), (2.22)

    where

    σ(x,s1,s2,k,x0)=s1+s2s12(1+tanh(k(xx0)))

    is a sigmoid function transitioning from s1 to s2 around the threshold x0 with slope k.

    We restricted the fitting interval of ΔEresp to [ΔEresp(ε),ΔEresp(Ntotε)], with ε=102, as this variable goes to ± for very low concentrations of NADH or NAD+. The range of Δψ was set to [100,200] mV. In this state space, the relative l2 error of the fit, which is given by Joriginal˜Jfittedl2Joriginall2, was 0.014 (see Figure 2(a)).

    Figure 2.  Fit of the flux JO (2(a)) and JF1F0 (2(b)) defined in Eqs (2.22) and (2.23), respectively. The surface height quantifies the value of the flux and the fit error is computed as the absolute difference between the fitted and the original expressions.

    ADP phosphorylation / ATP hydrolysis rate. Similarly, the ADP phosphorylation rate JF1F0 from Eq (2.8) has a sigmoidal profile as function of each of its variables ΔGp,m and Δψ. We proposed the following simplified expression:

    ˜JF1F0(ΔGp,m,Δψ)=σ(ΔGp,m,p10,0,p11,p12)σ(Δψ,1,p15,p13,p14). (2.23)

    The fit was performed on the same interval for Δψ, and [ΔGp,m(ε),ΔGp,m(Atot,mε)] for the Gibbs free energy. The relative l2 error for this fit was 0.056 (see Figure 2(b)).

    Substrate metabolization. The substrate metabolization flux expression (2.6) has 5 parameters. Using the change of variable (2.13), and plotting the resulting expression as function of its variables ([Ca2+]m,ΔEresp), we observed that it could be fitted by a simplified expression with 4 parameters:

    ˜Jsub=p6(1p7exp(p8[Ca2+]m))1exp(FRTΔErespln(Kresp))+p9 (2.24)

    The mitochondrial calcium concentration range of variability was set to [0,4×104] nmol mgprot−1, and the l2 relative difference due to the fit was 0.006 (see Figure 3(a)).

    Figure 3.  Fit of the flux Jsub (3(a)) and Juni (3(b)) defined in Eqs (2.24) and (2.25), respectively.

    Calcium uniporter. Plotting the flux through the calcium uniporter as a function of Δψ and [Ca2+]c suggested a product of two linear functions as a simplification. A limitation is that the calcium uniporter operates generally in one direction (from the cytoplasm to the matrix), and thus takes only positive values. For this reason, we chose to take only the positive part of that product:

    ˜Juni([Ca2+]c,Δψ)=max(0,p24(Δψp25)([Ca2+]cp26)) (2.25)

    Even though the simplified expression is not differentiable everywhere in its range of state variables, the original values were properly fitted with a relative l2 error of 0.037 (see Figure 3b). The range of variation of the cytoplasmic calcium concentration was set to [ε,5] µM.

    Other expressions. The remaining flux expressions JNaCa,JANT and Jleak were kept unchanged, with adaption to the new variables and grouping of their parameters:

    ˜JNaCa([Ca2+]m)=p22[Ca2+]mp23+[Ca2+]m, (2.26)
    ˜JANT(ΔGp,m,Δψ,ΔGp,c)=p16(1p17exp(FRT(ΔGp,cΔGp,mΔψ)))(1+p17p18exp(FRT(ΔGp,cp19Δψ)))(1+p18exp(FRTΔGp,m)), (2.27)
    ˜Jleak(Δψ)=p20exp(p21Δψ). (2.28)

    The parameter p20 and p21 associated the the proton leakage flux were computed by fitting the expressions (2.28) to the data from [22].

    By deriving the ODE system (2.18) with the flux expressions described in the previous section, we rewrote a classic model of mitochondrial respiration and calcium handling with thermodynamical variables instead of concentrations. The complex functions from biochemical modeling describe all the molecular mechanisms involved in all flux. They were replaced by simpler phenomenological functions, in which the molecular mechanisms cannot be identified anymore. Yet these functions compare to the original ones up to a 5.6% difference in relative l2 error and allowed for easier parameter calibration.

    In order to solve numerically the ODE system (2.18) with the simplified flux expressions, we chose a robust, yet simple, numerical predictor-corrector scheme. The numerical method is described as follows: consider a system of first order ordinary differential equations

    y=f(t,y),y(t0)=y0,

    an initial guess is computed via the Euler method, then this guess is improved using trapezoidal rule:

    ˜yi+1=yi+hf(ti,yi),yi+1=yi+12h(f(ti,yi)+f(ti+1,˜yi+1)),

    where h is the time step and was chosen equal to 105min. This solver was implemented from scratch in Python (available in the supplementary materials), and produced an output in 100 seconds. Special care was taken for the variables ΔEresp and ΔGp,m, since the change of variable introduced in section 2.2.1 maps the variables [NADH]m and [ADP]m in the bounded intervals (0,Ntot) and (0,Atot,m) to ΔEresp and ΔGp,m respectively, which both belong to the unbounded interval (,). Hence it could be numerically problematic if not treated with caution, specially for very low values of NADH and ADP.

    The experimental data were obtained as detailed in [23]. Cardiac mitochondria were extracted from Wistar Male rats (from Janvier Labs, France). Populations of isolated mitochondria were analyzed within the day following their preparation and subsequent determination of protein concentration by the Bratford method. All effectors of the respiratory chain used herein were first prepared as concentrated mother solutions (500 mM for the substrates of the respiratory chain, glutamate and malate, and 100 mM for ADP; all compounds were purchased from Sigma France) in the specific respiration buffer for mammalian cardiac mitochondria. All solutions were kept on ice during the experiments.

    The mitochondrial oxygen consumption rates were collected by chronoamperometry with a Clark electrode [24], which is the gold standard in bioenergetics. This electrochemical system determines the dissolved oxygen concentration of a liquid medium by measuring the current resulting from the 4-electron reduction of O2 to H2O on a platinum disk (at a potential of -0.8 V vs the internal Ag/AgCl reference electrode). O2 diffused towards the platinum through a Teflon membrane protecting it from the solution. The Clark electrode was inserted in an oxygraphy chamber filled with 6 mL of respiration buffer under constant stirring, and thermostated at 28 degree C with a water gasket. Substrates of the respiratory chain were injected at 5 mM final concentration, then mitochondria were injected to reach a concentration of 0.4 mgprot mL−1 in the respiration buffer. ADP was further introduced at various concentrations (0.166 to 1 mM). The current at the Clark electrode was continuously measured, and the oxygen concentration determined by a simple proportionality factor (210 µM of O2 maximum concentration). Figure 4 shows the experimental setup used to acquire data.

    Figure 4.  Oxygraphy chamber used to collect respiration data.

    Among the available data, we selected those with substrates containing glutamate and malate. We could not select the data with succinate because it activates another complex of the respiratory chain that was not included in our model. We calibrated our model to the data of 5 representative experiments with different time instants and concentrations of ADP additions. A sixth dataset was used to challenge the predictive properties of the model. Each dataset gathers observations on the complete population of mitochondria in the oxygraphy chamber.

    As explained in section 3.1, the measured experimental data reflect the evolution of the concentration of O2 due to mitochondrial respiration, following successive ADP addditions. Since the cytoplasmic calcium is not taken into account in this set of experiments, the calcium source term [Ca2+]ext(t) is set to zero in the ODE system (2.18).

    Each cytoplasmic ADP addition is modeled through the source term JADP,ext(t), where

    JADP,ext=[ADP]c(t=0)+i[ADP]+c,i(1exp(tt+ADP,iτ+ADP))1t>t+ADP,i(t), (3.1)

    with parameters [ADP]+c,i and t+ADP,i representing the ADP concentration and the time instant of the i-th ADP addition into the system, respectively. The parameter τ+ADP accounts for the diffusion of the ADP from the tip of the syringe to the mitochondria. This time is short (<5 seconds) compared to the duration of the experiment, as there is an agitator at the bottom of the chamber. Although the parameters [ADP]+c,i and t+ADP,i are given by the experimental protocol, τ+ADP is an uncertain parameter since it depends on the way the addition was executed, and can vary from one experiment to another. Nevertheless, in our simulations we set τ+ADP=0.05 min, a reasonable constant value, since this parameter appeared not to be influential on the output of our model (see Figure 8).

    In order to complete the link between the experimental data and the observables of our model, we first computed the solution of the system (2.18) over the time interval [t0,tf], where t0 is the time where isolated mitochondria were added and tf the end of the experiment. Then, substituting the state variables (ΔEresp(t),Δψ(t)) of the solution in (2.7), we computed the single oxygen atom consumption rate over the time interval.

    In general, the results we obtained from simulations were close to the experimental data up to a multiplicative factor (usually between 30 and 40), accounting for various elements of the experiments that are not included in the model (e.g., the electrochemical process at the O2 electrode). Moreover, our model does not account for the variability shown by mitochondria populations between experiments (different solutions, differentlife times, etc.). For these reasons, our analysis focused on reproducing the ratios of respiration rates between different mitochondrial states rather than the absolute values of these rates.

    The parameter analysis and the calibration algorithm that will be described in sections 4 and 5 both rely on a cost function, which in our case will measure the discrepency between the output of our model and experimental data. The choice of a cost function for these methods determines the properties of the solution. In particular, we are interested in reproducing the mitochondrial states that can be clearly identified on the experimental data (see Figure 5). These states are characterized by activity rates that are nearly constant and by short transitions that occur in between [25]. We then used a cost function designed specifically to reproduce those states.

    Figure 5.  Targeted oxygen consumption rates JO,k (orange slopes) are plotted on top of experimental data (blue). Intervals Ik on which linear regressions were performed were picked manually for each experiment. On I1 and I3, after ADP additions, mitochondria are in the so called state 4 [25]. They are in state 3 on I2, after ADP from the first addition is depleted. On I0, mitochondria are in a state which is close to state 3 and which we call "substrate state", as there is no external ADP up to that point.

    Let us note Ik the intervals on which mitochondrial states can be identified and JO,k the corresponding respiration rates, which are obtained by linear fitting of the experimental data over the interval Ik. As explained in the previous section, we cannot compare directly the simulated and experimental absolute values of JO,k, therefore each rate for k>0 is normalized by the value for k=0.

    Our goal is to obtain simulated respiration rates JO,k, k>0, that are nearly constant on each interval Ik, with an accurate value. Then, the cost function can be formally written as:

    D(P)=w1k>0(¯JO,k(P)¯JO,0(P)JO,kJO,0)2+w2k>01|Ik|Ik1ddtJO,k(P)>ε(s)ds. (3.2)

    The first term of D is the difference between the experimental rates fitted over each interval Ik, and the simulated rates (¯JO,k) fitted over the same intervals. The second term measures the relative duration for which the derivative of JO,k is larger than a tolerance ε. The smaller this term gets, the more linear JO,k(P,t) is. The weights w1 and w2 balance the two terms. They were set to w1=0.04 and w2=1.

    Even though the proposed model in section 2.2.2 is simple and has only 32 parameters, some of these parameters may have small influence on the output of the flux or on the respiratory rates computed through the model. Setting these parameters to a fixed value somewhere in their range of variability reduces the number of uncertain parameters and hence is helpful for the calibration procedure. Amongst the wide variety of tools to perform sensitivity analysis, the Sobol method is effective for quantifying the influence of a certain parameter on the output of a model. Using the concept of variance decomposition of the output with respect to the parameters [26], it allows for global exploration of the parameters space, while taking into account their statistical distribution.

    In practice, our model is an equation that can be written as Y=f(p1,p2,,pk) where Y is the output and pi,i{1,,k} are the parameters. The total Sobol index for a parameter pi is STi=E[V(Y/pi)]/V(Y), where pi denotes all the parameters except pi. This index quantifies the influence of the parameter pi on the output of the model Y, relatively to the other parameters: the larger the value of the index, the larger the influence on the output. If STi0, then pi can be fixed to any value in its distribution, without significantly affecting V(Y).

    Calculating the Sobol index associated to a parameter pi is computationally expensive, because a large sampling of parameters is required, leading to a large number of model evaluations. To overcome this, we used the Saltelli sampling method [26] which reduces the total number of evaluations of f to (k+2)Ns, where Ns is the number of parameter samples. All the Sobol indices were computed using a sample size Ns=2000, thanks to the SALib Python library [27].

    In the next two sections, we first report the Sobol analysis on our simplified flux and rates expressions, then on the whole system (2.18). The first step gave us an indication on the relative influence of each parameter on its associated rate expression. The second step highlighted which parameters had an influence on the respiration rates from the model, that were compared to the experimental ones thanks to the cost function D, given in section 3.3.

    In this section, we consider each flux from section 2.2.2 as a separate model that depends on state variables and some parameters. For each flux and each parameter, Sobol indices were computed using the same ranges of state variables as for the surface fitting procedure (see Figure 6). We then computed the average and maximum of the total Sobol indices along a batch of 50 trajectories that were generated with random sets of parameters, each parameter being picked in a range of ±10% around the value determined previously from the fits. The resulting indices are presented in Figure 7. From this bar graph, we set a minimum Sobol index of 0.1, under which parameters were considered constant without changing the outcome of the model. We then excluded the parameters p0, p1, p3, p8, p10, p13, p15, p20 and p24 from our analysis.

    Figure 6.  Total Sobol index for the parameter p11 of ˜JF1F0 defined in Eq (2.23) plotted in the phase space (ΔGp,m,Δψ). White lines are trajectories computed with random samples of the whole set of parameters, on which we measured the Sobol index.
    Figure 7.  Average and maximum values of total Sobol indices for all parameters appearing in a flux function (Eqs (2.22) to (2.27)). Each parameter contributes to a single flux function. Sobol indices were evaluated on 50 trajectories generated with random parameters.

    The final objective of our sensitivity analysis was to determine which were the influential parameters of the model, in order to calibrated them versus the experimental measurements of mitochondrial oxygen consumption rates (sec. 3.1). Having fixed the nine parameters that were identified in the previous section, we used the cost function D from Eq 3.2 as output of the model.

    With a reduced total of 23 parameters (32-9), and using Ns=2000 samples, the computation of the Sobol indices required 50000 model evaluations. Even though this computational cost is high, each evaluation is independant of the others, so the task can be embarrassingly parallelized.

    The results of the Sobol analysis are presented in Figure 8. Six parameters only (p2,p5,p11,p12,p14, and p21) appeared to have important influence on the output D.

    Figure 8.  Total Sobol index representing the sensitivity of the cost function D defined in Eq (3.2) with respect to its parameters. The conditions used to compute these indices are of two consecutive additions of ADP at concentrations of 0.66 mM and 1 mM.

    The calibration of the 6 previously identified parameters p2,p5,p11,p12,p14, and p21 was required to be able to reproduce the experimental measurements as accurately as possible.

    In addition to the six above parameters, we calibrated the parameter γ that accounts for the volume ratio between the mitochondria and the cytoplasm, and for the flux unit conversions.

    The dimensionality of the problem and the non-smoothness of our model make the use of gradient-based calibration method not trivial. For this reason, we choose a genetic optimization algorithm to calibrate the parameters. We implemented the algorithm ourselves in Python (available in the supplementary materials).

    Each generation of the algorithm consisted in a set of 1800 samples of parameters. As previously done for the sensitivity analysis, parameters were picked in a range of ±10 % around their base value, computed initially from surface fitting (see section 2.2.2). The function to be minimized is the cost function (3.2) that was discussed in section 3.3. The selection of sets of parameters that occurs at each iteration of the algorithm was based on truncation, where only the better half of the population (i.e. with lower values of the cost function) will eventually be passed to the reproduction phase. During this phase, crossover and mutation of parent parameters are applied with a probabilty equal 0.7 and 0.05 respectively. In the crossover operation, the weighting coefficient between the two selected sets of parameters is picked randomly in (0, 1). The minimization was performed independently for the 5 experiments described in section 3.1, to account for the variability of the experiments (mitochondria, measurements, etc.). The algorithm converged in 21 to 25 iterations. Convergence was considered attained when |¯D(Pj+1)¯D(Pj)|105, where Pj is the population of the j-th generation of the parameters {p2,p5,p11,p12,p14,p21,γ}. The functional D decreases from 5.33 to 2.52 on average. We run 2-3 repeats of the algorithm for each experiment, to reduce the probability of obtaining a local minimum.

    After having calibrated the parameters, we computed their average values and standard deviation across experiments. We also ran our model with these individual sets of parameters in order to make a direct comparison with the experimental data.

    The values of the calibrated parameters are reported in Table 3, averaged over the five available experiments. With the exception of the volume ratio γ, the variability of the parameters across experiments is low (2%). This suggests that our model is robust for the modeling experiments on populations of isolated mitochondria.

    Figure 9 shows the comparison between our model run with calibrated parameters and respiration data. Our model reproduces the quick transitions between respiratory states of the mitochondria.

    Figure 9.  Calibrated oxygen consumptions (dashed red lines) compared to data from 5 experiments (blue solid lines). Arrows locate additions of mitochondria then ADP at two time instants.

    However, when the small quantities of ADP are added (Figure 9(d) and 9(e)), it is more difficult to identify the slope corresponding to mitochondrial state 3, and our model hardly reproduces the experimental data. This is confirmed when plotting the simulated respiration rates on top of the values fitted from data, as shown in Figure 10. It is also likely that the fitted respiration rate could not be one of a theoretical population of mitochondria which are all in the same state at all times. Nevertheless, the simulations correctly identified the periods of time when mitochondria are in state 3 or 4, even for these difficult conditions.

    Figure 10.  Calibrated respiration rates (red) compared to the rates fitted on experimental data (black), as described in section 3.3. The target respiration rates are only plotted on the intervals used to fit experimental data, and also used in the computation of the cost function.

    We recall that the cost function used for calibration is derived to recover the ratios between respiration rates, instead of their absolute values. Thus, the simulated results plotted in Figures 9 and 10 were scaled by a constant value determined during the initial mitochondrial inactivated state (before any addition).

    Additionally, the very sharp transitions of the simulations shown in Figure 10, cannot be seen on the experiments. This is expected as well, since the experimental measurement system displays a response-time (of about 1 second) and because a population effect tends to smoothen the total respiration of all mitochondria.

    In order to assess the predictive properties of our model, we ran a simulation with the averaged parameters given in Table 3. In Figure 11, we compare the output of the model with a sixth separated experimental dataset, which was not used for calibration. There is good agreement between the two time series. It is confirmed by the fact that the value of the cost function is 2.44, which is lower than the average 2.52 obtained at the end of the calibration.

    Figure 11.  Comparison between experimental data and a simulation using average parameters calibrated from five other experiments. 11(a): oxygen concentration, 11(b): respiration rates.

    We proposed a simple model of cardiac mitochondrial activity, including ATP production through oxidative phosphorylation and calcium handling. The model has simple mathematical expressions that were fitted to revised versions of the more complex MK equations.

    The fact that the parameters of the MK model were set for pancreatic β cells can be seen as a limitation of our model. However, all those parameters are now grouped under our new flux functions parameters, which were eventually calibrated to respiration data of cardiac mitochondria from rats. Therefore, the origin of the parameters is not a major concern as we showed that they can be fitted to other types of cells. Nevertheless, parameters associated with the calcium handling were not calibrated in this work, because of lack of associated data. The ability of the model for accurate predictions for calcium associated mechanisms is yet to be confirmed.

    Another possible limitation of our simplified model is that its parameters do not have a direct biological meaning, like the parameters in models based on physics or physiology. However, this is offset by the robustness of our model, which cannot be overfitted. This was confirmed by the sensitivity analysis that was performed on the flux expressions, where only 9 parameters appeared to be non-influential (at least for the experimental setups that we considered). This indicates that the simplified expressions balance correctly the complexity of the observations.

    We further made use of the available experimental data to perform a global sensitivity analysis, based on Sobol indices, on the respiration rates. In this step, a lot more parameters were excluded, which is explained by the sparsity of the experimental data at hand. In fact, this second sensitivity analysis step was an unformal indicator of the non-identifiability of most of the parameters, with respect to respiration rates data. The choice for the Sobol method was motivated by its ability to hightlight the interactions between the parameters. We also wanted to use a method that was able to analyze the entirety of our set of parameters in one go, and to account for the fact that they are statistically distributed. Our model with few parameters allowed for the use of such a computationally expensive method. With a very descriptive but highly complex model, we may have been able to conduct an analysis on a subset of parameters only, as done for example by Wu et al. [28].

    In our work, the seven influential parameters that were indentified from the sensitivity analysis were calibrated seperately, using an optimization algorithm, to five data sets of respiration rates controlled by different concentrations of ADP addition. Then the prediction capability of our model was assessed on a sixth data set that was not used in the calibration. Results show that our model is capable of reproducing respiratory states of cardiac mitochondria, and the transition between them, as described by biologists.

    The methodology that we developed and tested with experimental data on mitochondrial respiration activities can be applied to other types of datasets. For instance, other parameters can be determined with data on external calcium or ATP/ADP concentrations. Ultimately, all the parameters of the model may be determined, with the certainty that the model will not be overfitted.

    We acknowledge Audrey Sémont and Emma Abell from IHU Liryc, who prepared the mitochondria that were used in the experiments. We also thank Mark Potse for his help all along this work.

    This study received financial support from the French Government as part of the "Investments of the Future" program managed by the Agence Nationale de la Recherche grant ANR-10-IAHU-04 and also from the project MITOCARD (ANR-17-CE11-0041). Sobol analyses and calibrations were carried out using the PlaFRIM experimental computing testbed, supported by Inria, CNRS (LaBRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr/).

    All authors declare no conflicts of interest in this paper.

    Table 4.  Parameters of the rate of production of NADH by pyruvate dehydrogenase (2.6).
    Parameter Description Value Unit Reference
    u1 Parameter for the dependence of PDH-P dephosphorylation on Mg2+ 1.5 dimensionless [15]
    u2 Ratio between the maximal rate of kinase and the maximal rate of phosphatase 1 dimensionless [15]
    KCa Affinity for the PDH-P phosphatase 0.5×10-1 µM [15]
    Jmax Maximum rate for PDH flux 7.157×101 nmol mgprot−1 min−1 [5]
    kPDHnad Parameter for the dependence of JPDH on NAD 1 dimensionless [5]

     | Show Table
    DownLoad: CSV
    Table 5.  Parameters of the respiration driven proton pump flux and the oxygen consumption rate, involved in expressions (B.1) to (B.6).
    Parameter Description Value Unit Reference
    α12 Unimolecular rate constant for the transition 12 400 s−1 [9]
    α21 Unimolecular rate constant for the transition 21 5 s−1 [9]
    α23 Bimolecular rate constant for the transition 23 3.405×104 nmol−1/2 mgprot1/2 s−1 [4]
    α32 Unimolecular rate constant for the transition 32 8×104 s−1 [9]
    α34 Unimolecular rate constant for the transition 34 400 s−1 [9]
    α43 Bimolecular rate constant for the transition 43 1.59×102 nmol−1/2 mgprot1/2 s−1 [4]
    α45 Unimolecular rate constant for the transition 45 40 s−1 [9]
    α54 Unimolecular rate constant for the transition 54 0.4 s−1 [9]
    α56 Unimolecular rate constant for the transition 56 100 s−1 [9]
    α65 Unimolecular rate constant for the transition 65 1×105 s−1 [9]
    α16,0 Rate constant for the transition 16 at Δψ=0 130 s−1 [9]
    α61,0 Rate constant for the transition 61 at Δψ=0 1×1012 s−1 [9]
    ρresp Proton pump concentration 0.4 nmol mgprot−1 [4]
    ΔψB Boundary potential parameter 50 mV [9]

     | Show Table
    DownLoad: CSV
    Table 6.  Parameters of the oxidative phosphorylation rate and the proton uptake flux, involved in expressions (2.8) and (B.7).
    Parameter Description Value Unit Reference
    β12 Unimolecular rate constant for the transition 12 400 s−1 [9]
    β21 Unimolecular rate constant for the transition 21 40 s−1 [9]
    β23 Bimolecular rate constant for the transition 23 5×103 mgprot nmol−1 s−1 [9]
    β32 Unimolecular rate constant for the transition 32 5×103 s−1 [9]
    β34 Unimolecular rate constant for the transition 34 100 s−1 [9]
    β43 Bimolecular rate constant for the transition 43 5×104 M−1 mgprot nmol−1 s−1 [4]
    β45 Unimolecular rate constant for the transition 45 100 s−1 [9]
    β54 Unimolecular rate constant for the transition 54 100 s−1 [9]
    β56 Unimolecular rate constant for the transition 56 1000 s−1 [9]
    β65 Unimolecular rate constant for the transition 65 1000 s−1 [9]
    β16,0 Rate constant for the transition 16 at Δψ=0 4.98×107 s−1 [9]
    β61,0 Rate constant for the transition 61 at Δψ=0 100 s−1 [9]
    ρF1 Proton pump concentration 0.7 nmol [4]

     | Show Table
    DownLoad: CSV
    Table 7.  Adenine nucleotide translocator parameters of Eq (2.9).
    Parameter Description Value Unit Reference
    Jmax,ANT Maximum exchange rate 1000 nmol mgprot−1 min−1 [4]
    f Fraction Δψ 0.5 dimensionless [4]

     | Show Table
    DownLoad: CSV
    Table 8.  Na+/Ca2+ exchanger parameters of Eq (2.10).
    Parameter Description Value Unit Reference
    Jmax,NaCa Maximum rate for the Na/Ca exchanger 25 nmol mgprot−1 min−1 [4]
    KNa Michaelis-Menten constant 9.4 mM [4]
    KCa Michaelis-Menten constant 0.003 nmol [4]
    [Na]c Cytoplasmic sodium concentration 30 mM [4]

     | Show Table
    DownLoad: CSV
    Table 9.  Calcium uniporter parameters of Eq (2.11).
    Parameter Description Value Unit Reference
    Jmax,uni Maximum rate for the Ca2+ uniporter 300 nmol mgprot−1 min−1 [4]
    Ktrans Dissociation constant for uniporter translocated calcium 19 µM [4]
    Kact Dissociation constant for uniporter activating calcium 0.38 µM [4]
    na Activation cooperativity parameter 2.8 dimensionless [4]
    L Equilibrium constant for uniporter conformations 110 dimensionless [4]
    Δψ Δψ offset 91 mV [4]

     | Show Table
    DownLoad: CSV

    Let us recall the oxydation-reduction equation (2.1) of the respiratory chain that we considered in our model:

    1/2NADH+1/2H++1/4O21/2NAD++1/2H2O.

    This reaction is in fact a sum of three other reactions (occurring at complexes I, III, IV) where protons are pumped, since the oxydation of NADH involves the previous three complexes. Hence, our proton pump model takes into consideration all the protons pumped at three different energy conserving sites. We considered that for each electron transferred through the ETC, the total number of protons pumped in the three energy conserving sites is 6. This stoichiometry is coherent with the previously discussed model of Magnus and Keizer.

    The whole reaction can be represented schematically by the diagram in Figure 12. On this diagram:

    Figure 12.  The six state Altman-King-Hill diagram used to illustrate the respiration driven proton pump.

    ● The nodes represent possible states of the enzymatic complex.

    ● The edges represent possible transitions between these possible states.

    ● States on the left side are in the matrix and states on the right side are on the exterior of the mitochondria.

    ● The interval between the two sides can be seen as the mitochondrial membrane thickness.

    ● To each edge we can associate two rate constants αij and αji that represent the two possible transitions between states i and j. These rates can be either constants or functions of a variable of the ODE system.

    ● The arcs associated with some edges represent the binding or unbinding of molecules to the enzymatic complex.

    ● The completion of a single cycle in the counter clockwise direction lead to the ejection of 6H+ from inside the matrix to the outside.

    The dashed edge associating states 2 and 5 illustrates the possibility of reaction slip, also known as intrinsic uncoupling. Pietrobon and Caplan as well as Magnus and Keizer considered those slips in their models. However, we made the assumption that the reaction is fully coupled to the outward proton flow, and all possible proton leaks are taken into account in the inward leak added in section 2.2.2.

    The expression of the cycle rate (the average number of occurrences of the complete cycle per second) was discussed by Hill [8], and is given by:

    Jcycle=60NΠ+ΠΣ (B.1)

    where N is the total number of respiratory enzymatic complexes in the inner membrane and the coefficient 60 converts $ s−1 to min−1 . \Pi_\pm $ are the product of the rate constants of the complete cycle, in the counter clockwise and clockwise directions respectively:

    Π+=α12α23α34α45α56α61,Π=α16α65α54α43α32α21. (B.2)

    In order to define Σ, we need to introduce the Hill definition of what are called directional diagrams. For each state i of the cycle (in our case 6 states), we can associate a collection of directional diagrams for this state i, which illustrates all the possible ways leading to the state i taking into consideration 2 rules:

    ● The path towards state i must not be a complete cycle.

    ● The maximum number of edges must be used (5 in our case).

    As example, Figure 13 lists the six possible directional diagrams leading to state 1. States 2 to 5 also have six similar diagrams, built following the same principles. To each diagram corresponds a coefficient that is the product of the rates associated to the edges that constitute the path to state i. Then, we define Σi as the sum of these coefficients. For example:

    Figure 13.  The six possible directional diagrams that lead to state 1.
    Σ1=α23α34α45α56α61+α65α54α43α32α21+α34α45α56α61α21+α61α54α43α32α21+α56α61α43α32α21+α45α56α61α32α21. (B.3)

    The denominator Σ in (B.1) is then obtained as the sum of coefficients for all states:

    Σ=Σ1+Σ2+Σ3+Σ4+Σ5+Σ6 (B.4)

    Basically, Σ is the sum of 6×6=36 elements, and each element is the product of 5 rate constants. Even though that is a very large number of terms, it is still a simplification as we did not consider the intrinsic uncoupling in our model (i.e the slip transitions between the states 2 and 5), contrary to the original models of Pietrobon and Caplan or Magnus and Keizer.

    The reaction rate is hence a function of all rates of the cycle. Those can be either of first order (constants) or pseudo-first order (function of the ODE variables). Typically, the binding of a molecule to the enzymatic complex is associated to a pseudo first order rate constant. In our case, and since the respiratory variable is NADH, that would be the case for transitions (2)(3) and (4)(3). We can then write

    α23=α23NADH,α43=α43NAD+.

    The membrane potential plays a role in the (1)(6) transitions since it is where the conformational change in the enzymatic complexe happens i.e the reorientation of the binding sites. Following thermodynamics considerations, Pietrobon and Caplan were able to give expressions of the rate constants α16 and α61 as function of the membrane potential Δψ and the boundary potential parameter ΔψB:

    α16=α16,0e(nF/2RT)(ΔψΔψB) = α16(0)e(3FRT)(ΔψΔψB),α61=α61,0e(nF/2RT)(ΔψΔψB) = α61(0)e(3FRT)(ΔψΔψB), (B.5)

    where F=96480 Cmol1 is the Faraday constant, R=8.315 VCmol1, K1 is the perfect gaz constant and T=310.1 K is the temperature.

    The parameters α16,0 and α61,0 are the values of the rate contants α16 and α61, when Δψ=ΔψB, respectively. All other rates are considered constant as in the MK model.

    A summary of parameters involved in the expression of Jcycle is given in Table 5. In the final expression of the respiration flux, the parameter N is replaced by the concentration of proton pumps ρresp=0.4 nmol mgprot−1 introduced by Magnus and Keizer.

    Finally, as 6 protons are ejected and half a dioxygen is consumed during one cycle, the expressions of the rates used in our model are:

    JH,resp=12×JO=6ρrespJcycle. (B.6)

    We used the same method to derive the proton flux occuring during ATP synthesis as for the respiratory chain. The Hill-diagram used for this reaction correspond to a model of a six-state proton pump proposed by Pietrobon and Caplan [9] (see Figure 14). We made the same assumption as for the respiratory chain, by neglectic the leakage and slippage of the proton pump.

    Figure 14.  The six-state Altman-King-Hill diagram for an ATPase-driven proton pump. The proton stoichiometry is taken from the work of Pietrobon and Caplan [9].

    We note βij the rate constants associated to this diagram. As for the respiration driven proton pump model, all rates are first order, except those related to our ODE system variables. The latter ones are pseudo-first order rates, and they can be written as function of the state variables (Δψ,[ADP]m) as follows:

    β23= β23(Atot,m[ADP]m),β43= β43[ADP]m[Pi]m,β16= β16(0)exp(1.5FRT(ΔψΔψB))β61= β61(0)exp(1.5FRT(ΔψΔψB)).

    The intermediate quantities Π+,Π,Σ, and Jcycle are calculated the same way as in section B.1, only replacing the rate constants αij by the βij associated to the diagram 14.

    Since we express the rate of the oxidative phosphorylation and the proton uptake by mitochondria, a minus sign is added to the rates of ATP hydrolysis and its associated proton pump:

    JH,ATP=3JF1F0=3ρF1Jcycle, (B.7)

    where ρF1 is the concentration of the F1F0-ATPase-driven proton pump, in replacement of N in the definition (B.1) of Jcycle.



    [1] P. Diolez, V. Deschodt-Arsac, G. Raffard, C. Simon, P. D. Santos, E. Thiaudière, et al., Modular regulation analysis of heart contraction: application to in situ demonstration of a direct mitochondrial activation by calcium in beating heart, Am. J. Physiol.-Reg., I., 293 (2007), R13–R19.
    [2] M. Luo, M. E. Anderson, Mechanisms of altered Ca2+ handling in heart failure, Circ. Res., 113 (2013), 690–708.
    [3] A. Takeuchi, S. Matsuoka, Integration of mitochondrial energetics in heart with mathematical modelling, J. Physiol.-London, 598 (2020), 1443–1457.
    [4] G. Magnus, J. Keizer, Minimal model of β-cell mitochondrial Ca2+ handling, Am. J. Physiol.-Cell Ph., 273 (1997), C717–C733. doi: 10.1152/ajpcell.1997.273.2.C717
    [5] R. Bertram, M. G. Pedersen, D. S. Luciani, A. Sherman, A simplified model for mitochondrial ATP production, J. Theor. Biol., 243 (2006), 575–586. doi: 10.1016/j.jtbi.2006.07.019
    [6] C. C. Mitchell, D. G. Schaeffer, A two-current model for the dynamics of cardiac membrane, B. Math. Biol., 65 (2003), 767–793.
    [7] K. H. ten Tusscher, A. V. Panfilov, Alternans and spiral breakup in a human ventricular tissue model, Am. J. Physiol. Heart C., 291 (2006), H1088–1100.
    [8] T. L. Hill, Free energy transduction in biology: the steady-state kinetic and thermodynamic formalism, Academic Press, 1977.
    [9] D. Pietrobon, S. R. Caplan, Flow-force relationships for a six-state proton pump model: intrinsic uncoupling, kinetic equivalence of input and output forces, and domain of approximate linearity, Biochemistry-US, 24 (1985), 5764–5776. doi: 10.1021/bi00342a012
    [10] M.-H. T. Nguyen, M. S. Jafri, Mitochondrial calcium signaling and energy metabolism, Ann. N.Y. Acad. Sci., 1047 (2005), 127–137.
    [11] M.-H. T. Nguyen, S. J. Dudycha, M. S. Jafri, Effect of Ca2+ on cardiac mitochondrial energy production is modulated by Na+ and H+ dynamics, Am. J. Physiol.-Cell Ph., 292 (2007), C2004–C2020.
    [12] S. Cortassa, M. A. Aon, E. Marbân, R. L. Winslow, B. O'Rourke, An integrated model of cardiac mitochondrial energy metabolism and calcium dynamics, Biophys. J., 84 (2003), 2734–2755. doi: 10.1016/S0006-3495(03)75079-6
    [13] S. Cortassa, B. O'Rourke, R. L. Winslow, M. A. Aon, Control and regulation of mitochondrial energetics in an integrated model of cardiomyocyte function, Biophys. J., 96 (2009), 2466–2478.
    [14] L. D. Gauthier, J. L. Greenstein, S. Cortassa, B. O'Rourke, R. L. Winslow, A computational model of reactive oxygen species and redox balance in cardiac mitochondria, Biophys. J., 105 (2013), 1045–1056. doi: 10.1016/j.bpj.2013.07.006
    [15] G. Magnus, J. Keizer, Model of β-cell mitochondrial calcium handling and electrical activity. I. Cytoplasmic variables, Am. J. Physiol.-Cell Ph., 274 (1998), C1158–C1173.
    [16] M. D. Brand, The stoichiometry of the exchange catalysed by the mitochondrial calcium/sodium antiporter, Biochem. J., 229 (1985), 161–166.
    [17] T. Matsuda, K. Takuma, A. Baba, Na+-Ca2+ exchanger: physiology and pharmacology, Jpn. J. Pharmacol., 74 (1997), 1–20.
    [18] T. E. Gunter, D. R. Pfeiffer, Mechanisms by which mitochondria transport calcium, Am. J. Physiol.-Cell Ph., 258 (1990), C755–C786.
    [19] D. E. Wingrove, T. Gunter, Kinetics of mitochondrial calcium transport. I. Characteristics of the sodium-independent calcium efflux mechanism of liver mitochondria, J. Biol. Chem., 261 (1986), 15159–15165. doi: 10.1016/S0021-9258(18)66846-2
    [20] G. C. Brown, The leaks and slips of bioenergetic membranes, FASEB J., 6 (1992), 2961–2965.
    [21] M. D. Brand, L.-F. Chien, P. Diolez, Experimental discrimination between proton leak and redox slip during mitochondrial electron transport, Biochem. J., 297 (Pt 1) (1994), 27–29.
    [22] G. Krishnamoorthy, P. C. Hinkle, Non-ohmic proton conductance of mitochondria and liposomes, Biochemistry-US, 23 (1984), 1640–1645.
    [23] C. Colin, Développement de nouvelles méthodes de microscopie et d'électrochimie pour des analyses multi-paramétriques de mitochondries individuelles cardiaques (in French), Theses, Université de Bordeaux, 2020, Available from: https://tel.archives-ouvertes.fr/tel-03116368.
    [24] Carlos M. Palmeira, António J. Moreno (eds.), Mitochondrial Bioenergetics, vol. 1782 of Methods in Molecular Biology, Humana Press, 2008.
    [25] B. Chance, W. G.R., Respiratory enzymes in oxidative phosphorylation: III. The steady state, J. Biol. Chem., 217 (1955), 409–427. doi: 10.1016/S0021-9258(19)57191-5
    [26] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, et al., Global sensitivity analysis: the primer, John Wiley & Sons, 2008.
    [27] J. Herman, W. Usher, SALib: An open-source Python library for Sensitivity Analysis, J. Open Source Softw., 2 (2017), 97. doi: 10.21105/joss.00097
    [28] F. Wu, F. Yang, K. Vinnakota, D. Beard, Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology, J. Biol. Chem., 282 (2007), 24525–24537. doi: 10.1074/jbc.M701024200
  • mbe-18-05-291-s001.zip
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3385) PDF downloads(139) Cited by(0)

Figures and Tables

Figures(14)  /  Tables(9)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog