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


Introduction
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 (Ca 2+ ) cycle.Ca 2+ controls the contraction of myocytes as well as mitochondrial energy production [1].In particular, a dysfunction in mitochondrial Ca 2+ handling may very likely be associated to cardiac rhythm pathologies [2].The exact role of the mitochondria in regulating the intracellular Ca 2+ 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, Ca 2+ regulation through the calcium uniporter and the Na + /Ca 2+ 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.

Derivation of a new thermodynamical model with simplified flux functions
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 + /Ca 2+ exchange and Ca 2+ 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.
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: (2.1) We also kept the conservation of the total concentration of nicotinamide adenine dinucleotide (NAD) and adenine nucletotide in the mitochondrial matrix : The mitochondrial transmembrane potential difference ∆ψ and the calcium concentration [Ca 2+ ] m are also variables of the model, so that the dynamics of the mitochondrial activity can be modeled using the following ODE system: where C m and f m 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 C m , 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 C m and f m are given in Table 1.We detail in the next section the expressions of each flux J i .

Derivation of flux expressions
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.2.1.1.Metabolization of substrates 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 (PDH a ) 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 PDH a on mitochondrial calcium, as it is believed to stimulate the phosphatase reaction that produces PDH a .They proposed an expression of the fraction of PDH a : where u 1 is the dependence of PDH-P dephosphorylation on Mg 2+ , u 2 is the ratio between the maximal rate of kinase and the maximal rate of phosphatase, and K Ca is the Ca 2+ -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 J PDH and has the following expression: which parameters J max , u 1 , u 2 , K Ca , and K PDHnad are given in appendix in Table 4.The value of the parameter J max depends on the external substrate: glutamate, malate or both.

Oxygen consumption and proton pump
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 J H,resp and oxygen consumption J O : where J cycle 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 J cycle and its parameters are stated in section B.1 and in Table 5.

ADP phosphorylation -ATP hydrolysis
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.2were repeated to obtain the expressions of the rate of ATP hydrolysis and its associated proton pumping: where the coefficient 3 corresponds to the stoichiometry between ATP hydrolysis and its associated proton pump.The expression of J cycle is given in appendix B.2, and its parameters, as well as ρ F1 , in Table 6.

Adenine nucleotide translocator (ANT)
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: 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 handling
Calcium is exchanged between mitochondria and the cytoplasm through membrane channels and pores such as the Ca 2+ uniporter, the Na + /Ca 2+ 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 Ca 2+ uniporter and the Na + /Ca 2+ exchanger as calcium regulators.
We assume that the Na + /Ca 2+ exchanger is electro-neutral (two Na + for one Ca 2+ ) [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): Several studies proposed similar expressions for the Na + /Ca 2+ 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 Ca 2+ uniporter J uni depends on the transmembrane voltage ∆ψ and on the cytoplasmic calcium concentration [Ca 2+ ] c , as follows: The parameters of the calcium uniporter flux are given in Table 9.

External proton leakage
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.

Simplification of the model with thermodynamical variables
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.

Change of variables and new considerations
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: where K resp = 1.35 × 10 18 is the redox equilibrium constant of the reaction (2.1) [4].
Using the conservation relation (2.2), the NADH redox potential reads: (2.13) Gibbs free energy.It describes the energy balance between mitochondrial ATP and ADP.It is expressed in units of mV as follows: where K F1 = 1.71 × 10 6 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: We kept the variables for the transmembrane potential ∆ψ and the mitochondrial calcium concentration [Ca 2+ ] m since they are of main interest for us.
In addition to these four state variables, we also considered the cytoplasmic calcium concentration ([Ca 2+ ] c ), and the cytoplasmic ADP concentration ([ADP] c ) or equivalently the cytoplasmic Gibbs free energy change due to adenosine variations: 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: (2.17)
New ODE system.The new ODE system that models the cardiac mitochondrial activity was then written as follows: where and 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 J ADP,ext and [Ca 2+ ] ext are source terms for external ADP and Ca 2+ 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 mgprot −1 to mM.The parameter γ is simply 10

Surface fitting to simpler expressions
Expressions J i 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 J O defined in Eq (2.7) has a sigmoidal shape when plotted as a function of its variables (∆E resp , ∆ψ).The proposed simplified expression for this rate is: where is a sigmoid function transitioning from s 1 to s 2 around the threshold x 0 with slope k.We restricted the fitting interval of ∆E resp to ∆E resp (ε), ∆E resp (N tot − ε) , with ε = 10 −2 , 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 l 2 error of the fit, which is given by , was 0.014 (see Figure 2(a)).
ADP phosphorylation / ATP hydrolysis rate.Similarly, the ADP phosphorylation rate J F1F0 from Eq (2.8) has a sigmoidal profile as function of each of its variables ∆G p,m and ∆ψ.We proposed the following simplified expression: JF1F0 (∆G p,m , ∆ψ) = σ ∆G p,m , p 10 , 0, p 11 , p 12 σ (∆ψ, 1, p 15 , p 13 , p 14 ) . (2.23) The fit was performed on the same interval for ∆ψ, and ∆G p,m (ε), ∆G p,m (A tot,m − ε) for the Gibbs free energy.The relative l 2 error for this fit was 0.056 (see Figure 2(b)).Using the change of variable (2.13), and plotting the resulting expression as function of its variables ([Ca 2+ ] m , ∆E resp ), we observed that it could be fitted by a simplified expression with 4 parameters: The mitochondrial calcium concentration range of variability was set to [0, 4 × 10 −4 ] nmol mgprot −1 , and the l 2 relative difference due to the fit was 0.006 (see Figure 3(a)).
Calcium uniporter.Plotting the flux through the calcium uniporter as a function of ∆ψ and [Ca 2+ ] 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: Even though the simplified expression is not differentiable everywhere in its range of state variables, the original values were properly fitted with a relative l 2 error of 0.037 (see Figure 3(b)).The range of variation of the cytoplasmic calcium concentration was set to [ε, 5] µM.
Other expressions.The remaining flux expressions J NaCa , J ANT and J leak were kept unchanged, with adaption to the new variables and grouping of their parameters: The parameter p 20 and p 21 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 l 2 error and allowed for easier parameter calibration.

Solving the ODE system
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 an initial guess is computed via the Euler method, then this guess is improved using trapezoidal rule: where h is the time step and was chosen equal to 10 −5 min.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 ∆E resp and ∆G p,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, N tot ) and (0, A tot,m ) to ∆E resp and ∆G p,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.

Experiments
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 O 2 to H 2 O on a platinum disk (at a potential of -0.8 V vs the internal Ag/AgCl reference electrode).O 2 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 O 2 maximum concentration).Figure 4 shows the experimental setup used to acquire 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.

Link with our model
As explained in section 3.1, the measured experimental data reflect the evolution of the concentration of O 2 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 [Ca 2+ ] ext (t) is set to zero in the ODE system (2.18).
Each cytoplasmic ADP addition is modeled through the source term J ADP,ext (t), where 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 of the system (2.18) over the time interval [t 0 , t f ], where t 0 is the time where isolated mitochondria were added and t f the end of the experiment.Then, substituting the state variables (∆E resp (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 O 2 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.

Choice of cost function for sensitivity analysis and data fitting
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.
Let us note I k the intervals on which mitochondrial states can be identified and J * O,k the corresponding respiration rates, which are obtained by linear fitting of the experimental data over the interval I k .As explained in the previous section, we cannot compare directly the simulated and experimental absolute values of J * O,k , therefore each rate for k > 0 is normalized by the value for k = 0. Our goal is to obtain simulated respiration rates J O,k , k > 0, that are nearly constant on each interval I k , with an accurate value.Then, the cost function can be formally written as: The first term of D is the difference between the experimental rates fitted over each interval I k , and the simulated rates (J O,k ) fitted over the same intervals.The second term measures the relative duration for which the derivative of J O,k is larger than a tolerance ε.The smaller this term gets, the more linear J O,k (P, t) is.The weights w 1 and w 2 balance the two terms.They were set to w 1 = 0.04 and w 2 = 1.O,k (orange slopes) are plotted on top of experimental data (blue).Intervals I k on which linear regressions were performed were picked manually for each experiment.On I 1 and I 3 , after ADP additions, mitochondria are in the so called state 4 [25].They are in state 3 on I 2 , after ADP from the first addition is depleted.On I 0 , 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.

Sobol analysis principle
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 (p 1 , p 2 , . . ., p k ) where Y is the output and p i , i ∈ {1, . . ., k} are the parameters.The total Sobol index for a parameter p i is , where p ∼i denotes all the parameters except p i .This index quantifies the influence of the parameter p i 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 S T i ≃ 0, then p i can be fixed to any value in its distribution, without significantly affecting V(Y).
Calculating the Sobol index associated to a parameter p i 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)N s , where N s is the number of parameter samples.All the Sobol indices were computed using a sample size N s = 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.

Sobol analysis results on the flux functions
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 p 0 ,p 1 , p 3 , p 8 , p 10 , p 13 , p 15 , p 20 and p 24 from our analysis.

Sobol analysis results on the complete model
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 N s = 2 000 samples, the computation of the Sobol indices required 50 000 model evaluations.Even though this computational cost is high, each evaluation is independant of the others, so the task can be embarrassingly parallelized.

Genetic algorithm and a posteriori evaluation
The calibration of the 6 previously identified parameters p 2 , p 5 , p 11 , p 12 , p 14 , and p 21 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 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.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(P j+1 ) − D(P j ) ≤ 10 −5 , where P j is the population of the j-th generation of the parameters {p 2 , p 5 , p 11 , p 12 , p 14 , p 21 , γ}.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.

Results
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.
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.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 responsetime (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.

Conclusions
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.B. Derivation of flux expressions for the respiratory chain and the ATP-synthesis using Hill diagrams

B.1. Respiratory chain
Let us recall the oxydation-reduction equation (2.1) of the respiratory chain that we considered in our model: 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: • 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 α i j 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: where N is the total number of respiratory enzymatic complexes in the inner membrane and the coefficient 60 converts s −1 to min −1 .Π ± are the product of the rate constants of the complete cycle, in the counter clockwise and clockwise directions respectively: 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:

Mathematical Biosciences and Engineering
Volume 18, Issue 5, 5758-5789.The denominator Σ in (B.1) is then obtained as the sum of coefficients for all states: 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 where F = 96480 Cmol −1 is the Faraday constant, R = 8.315 VCmol −1 , K −1 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 J cycle 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: (B.6)

B.2. ATP synthesis
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.We note β i j 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: The intermediate quantities Π + , Π − , Σ, and J cycle are calculated the same way as in section B.1, only replacing the rate constants α i j by the β i j 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:

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

Figure 2 .
Figure 2. Fit of the flux J O (2(a)) and J F1F0 (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.

Figure 5 .
Figure 5. Targeted oxygen consumption rates J *O,k (orange slopes) are plotted on top of experimental data (blue).Intervals I k on which linear regressions were performed were picked manually for each experiment.On I 1 and I 3 , after ADP additions, mitochondria are in the so called state 4[25].They are in state 3 on I 2 , after ADP from the first addition is depleted.On I 0 , 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.

Figure 6 .
Figure 6.Total Sobol index for the parameter p 11 of JF1F0 defined in Eq (2.23) plotted in the phase space (∆G p,m , ∆ψ).White lines are trajectories computed with random samples of the whole set of parameters, on which we measured the Sobol index.

Figure 7 .
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.

Figure 8 .
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.

Figure 9 .
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.

Figure 10 .
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.

Figure 11 .
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.

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

Figure 13 .
Figure 13.The six possible directional diagrams that lead to state 1.

Figure 14 .
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].

J
H,ATP = 3J F1F0 = −3ρ F1 J cycle , (B.7) where ρ F1 is the concentration of the F1F0-ATPase-driven proton pump, in replacement of N in the definition (B.1) of J cycle .

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).

Table 3 .
Parameters of the model provided by the genetic algorithm, minimizing the cost function (3.2) on 5 different sets of experimental data.

Table 6 .
Parameters of the oxidative phosphorylation rate and the proton uptake flux, involved in expressions (2.8) and (B.7).