
The (2+1)-dimensional Chaffee-Infante equation (CIE) is a significant model of the ion-acoustic waves in plasma. The primary objective of this paper was to establish and examine closed-form soliton solutions to the CIE using the modified extended direct algebraic method (m-EDAM), a mathematical technique. By using a variable transformation to convert CIE into a nonlinear ordinary differential equation (NODE), which was then reduced to a system of nonlinear algebraic equations with the assumption of a closed-form solution, the strategic m-EDAM was implemented. When the resulting problem was solved using the Maple tool, many soliton solutions in the shapes of rational, exponential, trigonometric, and hyperbolic functions were produced. By using illustrated 3D and density plots to evaluate several soliton solutions for the provided definite values of the parameters, it was possible to determine if the soliton solutions produced for CIE are cuspon or kink solitons. Additionally, it has been shown that the m-EDAM is a robust, useful, and user-friendly instrument that provides extra generic wave solutions for nonlinear models in mathematical physics and engineering.
Citation: Naveed Iqbal, Muhammad Bilal Riaz, Meshari Alesemi, Taher S. Hassan, Ali M. Mahnashi, Ahmad Shafee. Reliable analysis for obtaining exact soliton solutions of (2+1)-dimensional Chaffee-Infante equation[J]. AIMS Mathematics, 2024, 9(6): 16666-16686. doi: 10.3934/math.2024808
[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 |
The (2+1)-dimensional Chaffee-Infante equation (CIE) is a significant model of the ion-acoustic waves in plasma. The primary objective of this paper was to establish and examine closed-form soliton solutions to the CIE using the modified extended direct algebraic method (m-EDAM), a mathematical technique. By using a variable transformation to convert CIE into a nonlinear ordinary differential equation (NODE), which was then reduced to a system of nonlinear algebraic equations with the assumption of a closed-form solution, the strategic m-EDAM was implemented. When the resulting problem was solved using the Maple tool, many soliton solutions in the shapes of rational, exponential, trigonometric, and hyperbolic functions were produced. By using illustrated 3D and density plots to evaluate several soliton solutions for the provided definite values of the parameters, it was possible to determine if the soliton solutions produced for CIE are cuspon or kink solitons. Additionally, it has been shown that the m-EDAM is a robust, useful, and user-friendly instrument that provides extra generic wave solutions for nonlinear models in mathematical physics and engineering.
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.
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++1╱2O2⟶NAD++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= (JPDH−JO),ddt[ADP]m= (JANT−JF1F0),ddtΔψ= (JH,resp−JH,ATP−JANT−Jleak−2Juni)/Cm,ddt[Ca2+]m= fm(Juni−JNaCa), | (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.
Parameter | Description | Value | Unit | Reference |
fm | Fraction of free matrix calcium concentration | 3 × 10−4 | dimensionless | [4] |
Cm | Membrane capacitance | 1.450 × 10−3 | nmol mV−1 mgprot−1 | [4] |
fc | Fraction of free cytosolic calcium concentration | 1 × 10−2 | 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] |
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.
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 |
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 |
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,ANT1−0.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(Δψ−Δψ∗)1−exp(−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)(Jsub−JO),ddtΔGp,m= φ2(ΔGp,m)(JF1F0−JANT),ddtΔGp,c= φ3(ΔGp,c)(γJANT−ddtJADP,ext(t)),CmddtΔψ= 12JO−3JF1F0−JANT−Jleak−2Juni,ddt[Ca2+]m= fm(Juni−JNaCa),ddt[Ca2+]c= fc˜γ(JNaCa−Juni)+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]mK−1F1exp(FRTΔGp,m))2Atot,m[Pi]mK−1F1exp(FRTΔGp,m), | (2.20) |
and
φ3(ΔGp,c)=RTF(1+[Pi]cK−1F1,cexp(FRTΔGp,c))2Atot,c[Pi]cK−1F1,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+s2−s12(1+tanh(k(x−x0))) |
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 ε=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 l2 error of the fit, which is given by ∣∣Joriginal−˜Jfitted∣∣l2∣∣Joriginal∣∣l2, was 0.014 (see Figure 2(a)).
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(1−p7exp(−p8[Ca2+]m))1exp(FRTΔEresp−ln(Kresp))+p9 | (2.24) |
The mitochondrial calcium concentration range of variability was set to [0,4×10−4] nmol mgprot−1, and the l2 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 [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+]c−p26)) | (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(1−p17exp(FRT(ΔGp,c−ΔGp,m−Δψ)))(1+p17p18exp(FRT(ΔGp,c−p19Δψ)))(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 10−5min. 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.
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(1−exp(−t−t+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.
Let us note Ik 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 Ik. 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 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)=w1∑k>0(¯JO,k(P)¯JO,0(P)−J∗O,kJ∗O,0)2+w2∑k>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/p∼i)]/V(Y), where p∼i 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 STi≃0, 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.
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.
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)|≤10−5, 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.
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 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.
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.
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] |
Parameter | Description | Value | Unit | Reference |
α12 | Unimolecular rate constant for the transition 1→2 | 400 | s−1 | [9] |
α21 | Unimolecular rate constant for the transition 2→1 | 5 | s−1 | [9] |
α∗23 | Bimolecular rate constant for the transition 2→3 | 3.405×104 | nmol−1/2 mgprot1/2 s−1 | [4] |
α32 | Unimolecular rate constant for the transition 3→2 | 8×104 | s−1 | [9] |
α34 | Unimolecular rate constant for the transition 3→4 | 400 | s−1 | [9] |
α∗43 | Bimolecular rate constant for the transition 4→3 | 1.59×102 | nmol−1/2 mgprot1/2 s−1 | [4] |
α45 | Unimolecular rate constant for the transition 4→5 | 40 | s−1 | [9] |
α54 | Unimolecular rate constant for the transition 5→4 | 0.4 | s−1 | [9] |
α56 | Unimolecular rate constant for the transition 5→6 | 100 | s−1 | [9] |
α65 | Unimolecular rate constant for the transition 6→5 | 1×105 | s−1 | [9] |
α16,0 | Rate constant for the transition 1→6 at Δψ=0 | 130 | s−1 | [9] |
α61,0 | Rate constant for the transition 6→1 at Δψ=0 | 1×1012 | s−1 | [9] |
ρresp | Proton pump concentration | 0.4 | nmol mgprot−1 | [4] |
ΔψB | Boundary potential parameter | 50 | mV | [9] |
Parameter | Description | Value | Unit | Reference |
β12 | Unimolecular rate constant for the transition 1→2 | 400 | s−1 | [9] |
β21 | Unimolecular rate constant for the transition 2→1 | 40 | s−1 | [9] |
β∗23 | Bimolecular rate constant for the transition 2→3 | 5×103 | mgprot nmol−1 s−1 | [9] |
β32 | Unimolecular rate constant for the transition 3→2 | 5×103 | s−1 | [9] |
β34 | Unimolecular rate constant for the transition 3→4 | 100 | s−1 | [9] |
β∗43 | Bimolecular rate constant for the transition 4→3 | 5×104 | M−1 mgprot nmol−1 s−1 | [4] |
β45 | Unimolecular rate constant for the transition 4→5 | 100 | s−1 | [9] |
β54 | Unimolecular rate constant for the transition 5→4 | 100 | s−1 | [9] |
β56 | Unimolecular rate constant for the transition 5→6 | 1000 | s−1 | [9] |
β65 | Unimolecular rate constant for the transition 6→5 | 1000 | s−1 | [9] |
β16,0 | Rate constant for the transition 1→6 at Δψ=0 | 4.98×107 | s−1 | [9] |
β61,0 | Rate constant for the transition 6→1 at Δψ=0 | 100 | s−1 | [9] |
ρF1 | Proton pump concentration | 0.7 | nmol | [4] |
Parameter | Description | Value | Unit | Reference |
Jmax,ANT | Maximum exchange rate | 1000 | nmol mgprot−1 min−1 | [4] |
f | Fraction Δψ | 0.5 | dimensionless | [4] |
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] |
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] |
Let us recall the oxydation-reduction equation (2.1) of the respiratory chain that we considered in our model:
1/2NADH+1/2H++1/4O2⟶1/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:
● 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:
Σ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=α∗23√NADH,α43=α∗43√NAD+. |
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 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 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.
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] |
A. Seadawy, A. Sayed, Soliton solutions of cubic-quintic nonlinear Schrödinger and variant Boussinesq equations by the first integral method, Filomat, 31 (2017), 4199–4208. http://dx.doi.org/10.2298/FIL1713199S doi: 10.2298/FIL1713199S
![]() |
[2] |
H. Yasmin, A. Alshehry, A. Ganie, A. Mahnashi, R. Shah, Perturbed Gerdjikov-Ivanov equation: soliton solutions via Backlund transformation, Optik, 298 (2024), 171576. http://dx.doi.org/10.1016/j.ijleo.2023.171576 doi: 10.1016/j.ijleo.2023.171576
![]() |
[3] |
A. Jawad, M. Petkovi, A. Biswas, Modified simple equation method for nonlinear evolution equations, Appl. Math. Comput., 217 (2010), 869–877. http://dx.doi.org/10.1016/j.amc.2010.06.030 doi: 10.1016/j.amc.2010.06.030
![]() |
[4] |
H. Khan, Shoaib, D. Baleanu, P. Kumam, J. Al-Zaidy, Families of travelling waves solutions for fractional-order extended shallow water wave equations, using an innovative analytical method, IEEE Access, 7 (2019), 107523–107532. http://dx.doi.org/10.1109/ACCESS.2019.2933188 doi: 10.1109/ACCESS.2019.2933188
![]() |
[5] |
H. Khan, S. Barak, P. Kumam, M. Arif, Analytical solutions of fractional Klein-Gordon and gas dynamics equations, via the (G'/G)-expansion method, Symmetry, 11 (2019), 566. http://dx.doi.org/10.3390/sym11040566 doi: 10.3390/sym11040566
![]() |
[6] |
H. Khan, R. Shah, J. Gomez-Aguilar, Shoaib, D. Baleanu, P. Kumam, Travelling waves solution for fractional-order biological population model, Math. Model. Nat. Phenom., 16 (2021), 32. http://dx.doi.org/10.1051/mmnp/2021016 doi: 10.1051/mmnp/2021016
![]() |
[7] |
K. Nisar, O. Ilhan, S. Abdulazeez, J. Manafian, S. Mohammed, M. Osman, Novel multiple soliton solutions for some nonlinear PDEs via multiple Exp-function method, Results Phys., 21 (2021), 103769. http://dx.doi.org/10.1016/j.rinp.2020.103769 doi: 10.1016/j.rinp.2020.103769
![]() |
[8] |
E. Yusufolu, A. Bekir, Solitons and periodic solutions of coupled nonlinear evolution equations by using the sine-cosine method, Int. J. Comput. Math., 83 (2006), 915–924. http://dx.doi.org/10.1080/00207160601138756 doi: 10.1080/00207160601138756
![]() |
[9] |
E. Fan, H. Zhang, A note on the homogeneous balance method, Phys. Lett. A, 246 (1998), 403–406. http://dx.doi.org/10.1016/S0375-9601(98)00547-7 doi: 10.1016/S0375-9601(98)00547-7
![]() |
[10] |
Z. Rahman, M. Zulfikar Ali, H. Roshid, Closed-form soliton solutions of three nonlinear fractional models through proposed improved Kudryashov method, Chinese Phys. B, 30 (2021), 050202. http://dx.doi.org/10.1088/1674-1056/abd165 doi: 10.1088/1674-1056/abd165
![]() |
[11] |
M. Ozisik, A. Secer, M. Bayram, On solitary wave solutions for the extended nonlinear Schrödinger equation via the modified F-expansion method, Opt. Quant. Electron., 55 (2023), 215. http://dx.doi.org/10.1007/s11082-022-04476-z doi: 10.1007/s11082-022-04476-z
![]() |
[12] |
Z. Lan, S. Dong, B. Gao, Y. Shen, Bilinear form and soliton solutions for a higher order wave equation, Appl. Math. Lett., 134 (2022), 108340. http://dx.doi.org/10.1016/j.aml.2022.108340 doi: 10.1016/j.aml.2022.108340
![]() |
[13] |
H. Hussein, H. Ahmed, W. Alexan, Analytical soliton solutions for cubic-quartic perturbations of the Lakshmanan-Porsezian-Daniel equation using the modified extended tanh function method, Ain Shams Eng. J., 15 (2024), 102513. http://dx.doi.org/10.1016/j.asej.2023.102513 doi: 10.1016/j.asej.2023.102513
![]() |
[14] |
G. Akram, M. Sadaf, S. Arshed, F. Sameen, Bright, dark, kink, singular and periodic soliton solutions of Lakshmanan-Porsezian-Daniel model by generalized projective Riccati equations method, Optik, 241 (2021), 167051. http://dx.doi.org/10.1016/j.ijleo.2021.167051 doi: 10.1016/j.ijleo.2021.167051
![]() |
[15] |
L. Li, E. Li, M. Wang, The (G'/G, 1/G)-expansion method and its application to travelling wave solutions of the Zakharov equations, Appl. Math. J. Chin. Univ., 25 (2010), 454–462. http://dx.doi.org/10.1007/s11766-010-2128-x doi: 10.1007/s11766-010-2128-x
![]() |
[16] |
X. Yang, Z. Deng, Y. Wei, A Riccati-Bernoulli sub-ODE method for nonlinear partial differential equations and its application, Adv. Differ. Equ., 2015 (2015), 117. http://dx.doi.org/10.1186/s13662-015-0452-4 doi: 10.1186/s13662-015-0452-4
![]() |
[17] |
S. Dai, Poincare-Lighthill-Kuo method and symbolic computation, Appl. Math. Mech., 22 (2001), 261–269. http://dx.doi.org/10.1007/BF02437964 doi: 10.1007/BF02437964
![]() |
[18] |
H. Durur, A. Kurt, O. Tasbozan, New travelling wave solutions for KdV6 equation using sub equation method, Applied Mathematics and Nonlinear Sciences, 5 (2020), 455–460. http://dx.doi.org/10.2478/amns.2020.1.00043 doi: 10.2478/amns.2020.1.00043
![]() |
[19] |
S. Bibi, S. Mohyud-Din, U. Khan, N. Ahmed, Khater method for nonlinear Sharma Tasso-Olever (STO) equation of fractional order, Results Phys., 7 (2017), 4440–4450. http://dx.doi.org/10.1016/j.rinp.2017.11.008 doi: 10.1016/j.rinp.2017.11.008
![]() |
[20] |
H. Ur Rehman, R. Akber, A. Wazwaz, H. Alshehri, M. Osman, Analysis of Brownian motion in stochastic Schrödinger wave equation using Sardar sub-equation method, Optik, 289 (2023), 171305. http://dx.doi.org/10.1016/j.ijleo.2023.171305 doi: 10.1016/j.ijleo.2023.171305
![]() |
[21] |
M. Alqhtani, K. Saad, R. Shah, W. Hamanah, Discovering novel soliton solutions for (3+1)-modified fractional Zakharov-Kuznetsov equation in electrical engineering through an analytical approach, Opt. Quant. Electron., 55 (2023), 1149. http://dx.doi.org/10.1007/s11082-023-05407-2 doi: 10.1007/s11082-023-05407-2
![]() |
[22] |
M. Al-Sawalha, S. Mukhtar, R. Shah, A. Ganie, K. Moaddy, Solitary waves propagation analysis in nonlinear dynamical system of fractional coupled Boussinesq-Whitham-Broer-Kaup equation, Fractal Fract., 7 (2023), 889. http://dx.doi.org/10.3390/fractalfract7120889 doi: 10.3390/fractalfract7120889
![]() |
[23] |
H. Yasmin, A. Alshehry, A. Ganie, A. Shafee, R. Shah, Noise effect on soliton phenomena in fractional stochastic Kraenkel-Manna-Merle system arising in ferromagnetic materials, Sci. Rep., 14 (2024), 1810. http://dx.doi.org/10.1038/s41598-024-52211-3 doi: 10.1038/s41598-024-52211-3
![]() |
[24] |
H. Yasmin, N. Aljahdaly, A. Saeed, R. Shah, Investigating symmetric soliton solutions for the fractional coupled Konno-Onno system using improved versions of a novel analytical technique, Mathematics, 11 (2023), 2686. http://dx.doi.org/10.3390/math11122686 doi: 10.3390/math11122686
![]() |
[25] |
S. El-Tantawy, H. Alyousef, R. Matoog, R. Shah, On the optical soliton solutions to the fractional complex structured (1+1)-dimensional perturbed Gerdjikov-Ivanov equation, Phys. Scr., 99 (2024), 035249. http://dx.doi.org/10.1088/1402-4896/ad241b doi: 10.1088/1402-4896/ad241b
![]() |
[26] |
L. Li, C. Duan, F. Yu, An improved Hirota bilinear method and new application for a nonlocal integrable complex modified Korteweg-de Vries (MKdV) equation, Phys. Lett. A, 383 (2019), 1578–1582. http://dx.doi.org/10.1016/j.physleta.2019.02.031 doi: 10.1016/j.physleta.2019.02.031
![]() |
[27] |
T. Han, Y. Jiang, Bifurcation, chaotic pattern and traveling wave solutions for the fractional Bogoyavlenskii equation with multiplicative noise, Phys. Scr., 99 (2024), 035207. http://dx.doi.org/10.1088/1402-4896/ad21ca doi: 10.1088/1402-4896/ad21ca
![]() |
[28] |
T. Han, Y. Jiang, J. Lyu, Chaotic behavior and optical soliton for the concatenated model arising in optical communication, Results Phys., 58 (2024), 107467. http://dx.doi.org/10.1016/j.rinp.2024.107467 doi: 10.1016/j.rinp.2024.107467
![]() |
[29] |
R. Ali, E. Tag-eldin, A comparative analysis of generalized and extended (G'/G)-Expansion methods for travelling wave solutions of fractional Maccari's system with complex structure, Alex. Eng. J., 79 (2023), 508–530. http://dx.doi.org/10.1016/j.aej.2023.08.007 doi: 10.1016/j.aej.2023.08.007
![]() |
[30] |
R. Ali, A. Hendy, M. Ali, A. Hassan, F. Awwad, E. Ismail, Exploring propagating soliton solutions for the fractional Kudryashov-Sinelshchikov equation in a mixture of liquid-gas bubbles under the consideration of heat transfer and viscosity, Fractal Fract., 7 (2023), 773. http://dx.doi.org/10.3390/fractalfract7110773 doi: 10.3390/fractalfract7110773
![]() |
[31] |
Y. Shi, C. Song, Y. Chen, H. Rao, T. Yang, Complex standard eigenvalue problem derivative computation for Laminar-Turbulent transition prediction, AIAA J., 61 (2023), 3404–3418. http://dx.doi.org/10.2514/1.J062212 doi: 10.2514/1.J062212
![]() |
[32] |
X. Cai, R. Tang, H. Zhou, Q. Li, S. Ma, D. Wang, et al., Dynamically controlling terahertz wavefronts with cascaded metasurfaces, Advanced Photonics, 3 (2021), 036003. http://dx.doi.org/10.1117/1.AP.3.3.036003 doi: 10.1117/1.AP.3.3.036003
![]() |
[33] |
A. She, L. Wang, Y. Peng, J. Li, Structural reliability analysis based on improved wolf pack algorithm AK-SS, Structures, 57 (2023), 105289. http://dx.doi.org/10.1016/j.istruc.2023.105289 doi: 10.1016/j.istruc.2023.105289
![]() |
[34] |
T. Ali, Z. Xiao, H. Jiang, B. Li, A class of digital integrators based on trigonometric quadrature rules, IEEE T. Ind. Electron., 71 (2024), 6128–6138. http://dx.doi.org/10.1109/TIE.2023.3290247 doi: 10.1109/TIE.2023.3290247
![]() |
[35] |
J. Dong, J. Hu, Y. Zhao, Y. Peng, Opinion formation analysis for expressed and private opinions (EPOs) models: reasoning private opinions from behaviors in group decision-making systems, Expert Syst. Appl., 236 (2024), 121292. http://dx.doi.org/10.1016/j.eswa.2023.121292 doi: 10.1016/j.eswa.2023.121292
![]() |
[36] |
C. Zhu, M. Al-Dossari, S. Rezapour, S. Shateyi, On the exact soliton solutions and different wave structures to the modified Schrödinger's equation, Results Phys., 54 (2023), 107037. http://dx.doi.org/10.1016/j.rinp.2023.107037 doi: 10.1016/j.rinp.2023.107037
![]() |
[37] |
C. Zhu, M. Al-Dossari, N. El-Gawaad, S. Alsallami, S. Shateyi, Uncovering diverse soliton solutions in the modified Schrödinger's equation via innovative approaches, Results Phys., 54 (2023), 107100. http://dx.doi.org/10.1016/j.rinp.2023.107100 doi: 10.1016/j.rinp.2023.107100
![]() |
[38] |
C. Zhu, S. Abdallah, S. Rezapour, S. Shateyi, On new diverse variety analytical optical soliton solutions to the perturbed nonlinear Schrödinger equation, Results Phys., 54 (2023), 107046. http://dx.doi.org/10.1016/j.rinp.2023.107046 doi: 10.1016/j.rinp.2023.107046
![]() |
[39] |
C. Zhu, S. Idris, M. Abdalla, S. Rezapour, S. Shateyi, B. Gunay, Analytical study of nonlinear models using a modified Schrödinger's equation and logarithmic transformation, Results Phys., 55 (2023), 107183. http://dx.doi.org/10.1016/j.rinp.2023.107183 doi: 10.1016/j.rinp.2023.107183
![]() |
[40] |
S. Lin, J. Zhang, C. Qiu, Asymptotic analysis for one-stage stochastic linear complementarity problems and applications, Mathematics, 11 (2023), 482. http://dx.doi.org/10.3390/math11020482 doi: 10.3390/math11020482
![]() |
[41] |
Y. Kai, Z. Yin, On the Gaussian traveling wave solution to a special kind of Schrödinger equation with logarithmic nonlinearity, Mod. Phys. Lett. B, 36 (2022), 2150543. http://dx.doi.org/10.1142/S0217984921505436 doi: 10.1142/S0217984921505436
![]() |
[42] |
Y. Kai, J. Ji, Z. Yin, Study of the generalization of regularized long-wave equation, Nonlinear Dyn., 107 (2022), 2745–2752. http://dx.doi.org/10.1007/s11071-021-07115-6 doi: 10.1007/s11071-021-07115-6
![]() |
[43] |
Y. Kai, Z. Yin, Linear structure and soliton molecules of Sharma-Tasso-Olver-Burgers equation, Phys. Lett. A, 452 (2022), 128430. http://dx.doi.org/10.1016/j.physleta.2022.128430 doi: 10.1016/j.physleta.2022.128430
![]() |
[44] |
S. Noor, A. Alshehry, A. Khan, I. Khan, Analysis of soliton phenomena in (2+1)-dimensional Nizhnik-Novikov-Veselov model via a modified analytical technique, AIMS Mathematics, 8 (2023), 28120–28142. http://dx.doi.org/10.3934/math.20231439 doi: 10.3934/math.20231439
![]() |
[45] |
I. Ullah, K. Shah, T. Abdeljawad, Study of traveling soliton and fronts phenomena in fractional Kolmogorov-Petrovskii-Piskunov equation, Phys. Scr., 99 (2024), 055259. http://dx.doi.org/10.1088/1402-4896/ad3c7e doi: 10.1088/1402-4896/ad3c7e
![]() |
[46] |
M. Bilal, J. Iqbal, R. Ali, F. Awwad, E. Ismail, Exploring families of solitary wave solutions for the fractional coupled higgs system using modified extended direct algebraic method, Fractal Fract., 7 (2023), 653. http://dx.doi.org/10.3390/fractalfract7090653 doi: 10.3390/fractalfract7090653
![]() |
[47] |
R. Ali, Z. Zhang, H. Ahmad, Exploring soliton solutions in nonlinear spatiotemporal fractional quantum mechanics equations: an analytical study, Opt. Quant. Electron., 56 (2024), 838. http://dx.doi.org/10.1007/s11082-024-06370-2 doi: 10.1007/s11082-024-06370-2
![]() |
[48] |
M. Akbar, N. Ali, J. Hussain, Optical soliton solutions to the (2+1)-dimensional Chaffee-Infante equation and the dimensionless form of the Zakharov equation, Adv. Differ. Equ., 2019 (2019), 446. http://dx.doi.org/10.1186/s13662-019-2377-9 doi: 10.1186/s13662-019-2377-9
![]() |
[49] |
R. Sakthivel, C. Chun, New soliton solutions of Chaffee-Infante equations using the exp-function method, Z. Naturforsch. A, 65 (2010), 197–202. http://dx.doi.org/10.1515/zna-2010-0307 doi: 10.1515/zna-2010-0307
![]() |
[50] |
T. Sulaiman, A. Yusuf, M. Alquran, Dynamics of lump solutions to the variable coefficients (2+1)-dimensional Burger's and Chaffee-infante equations, J. Geom. Phys., 168 (2021), 104315. http://dx.doi.org/10.1016/j.geomphys.2021.104315 doi: 10.1016/j.geomphys.2021.104315
![]() |
[51] |
S. Demiray, U. Bayrakcı, Construction of soliton solutions for Chaffee-Infante equation, AKU J. Sci. Eng., 21 (2021), 1046–1051. http://dx.doi.org/10.35414/akufemubid.946217 doi: 10.35414/akufemubid.946217
![]() |
[52] |
A. Seadawy, A. Ali, A. Altalbe, A. Bekir, Exact solutions of the (3+1)-generalized fractional nonlinear wave equation with gas bubbles, Sci. Rep., 14 (2024), 1862. http://dx.doi.org/10.1038/s41598-024-52249-3 doi: 10.1038/s41598-024-52249-3
![]() |
Parameter | Description | Value | Unit | Reference |
fm | Fraction of free matrix calcium concentration | 3 × 10−4 | dimensionless | [4] |
Cm | Membrane capacitance | 1.450 × 10−3 | nmol mV−1 mgprot−1 | [4] |
fc | Fraction of free cytosolic calcium concentration | 1 × 10−2 | 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] |
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 |
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 |
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] |
Parameter | Description | Value | Unit | Reference |
α12 | Unimolecular rate constant for the transition 1→2 | 400 | s−1 | [9] |
α21 | Unimolecular rate constant for the transition 2→1 | 5 | s−1 | [9] |
α∗23 | Bimolecular rate constant for the transition 2→3 | 3.405×104 | nmol−1/2 mgprot1/2 s−1 | [4] |
α32 | Unimolecular rate constant for the transition 3→2 | 8×104 | s−1 | [9] |
α34 | Unimolecular rate constant for the transition 3→4 | 400 | s−1 | [9] |
α∗43 | Bimolecular rate constant for the transition 4→3 | 1.59×102 | nmol−1/2 mgprot1/2 s−1 | [4] |
α45 | Unimolecular rate constant for the transition 4→5 | 40 | s−1 | [9] |
α54 | Unimolecular rate constant for the transition 5→4 | 0.4 | s−1 | [9] |
α56 | Unimolecular rate constant for the transition 5→6 | 100 | s−1 | [9] |
α65 | Unimolecular rate constant for the transition 6→5 | 1×105 | s−1 | [9] |
α16,0 | Rate constant for the transition 1→6 at Δψ=0 | 130 | s−1 | [9] |
α61,0 | Rate constant for the transition 6→1 at Δψ=0 | 1×1012 | s−1 | [9] |
ρresp | Proton pump concentration | 0.4 | nmol mgprot−1 | [4] |
ΔψB | Boundary potential parameter | 50 | mV | [9] |
Parameter | Description | Value | Unit | Reference |
β12 | Unimolecular rate constant for the transition 1→2 | 400 | s−1 | [9] |
β21 | Unimolecular rate constant for the transition 2→1 | 40 | s−1 | [9] |
β∗23 | Bimolecular rate constant for the transition 2→3 | 5×103 | mgprot nmol−1 s−1 | [9] |
β32 | Unimolecular rate constant for the transition 3→2 | 5×103 | s−1 | [9] |
β34 | Unimolecular rate constant for the transition 3→4 | 100 | s−1 | [9] |
β∗43 | Bimolecular rate constant for the transition 4→3 | 5×104 | M−1 mgprot nmol−1 s−1 | [4] |
β45 | Unimolecular rate constant for the transition 4→5 | 100 | s−1 | [9] |
β54 | Unimolecular rate constant for the transition 5→4 | 100 | s−1 | [9] |
β56 | Unimolecular rate constant for the transition 5→6 | 1000 | s−1 | [9] |
β65 | Unimolecular rate constant for the transition 6→5 | 1000 | s−1 | [9] |
β16,0 | Rate constant for the transition 1→6 at Δψ=0 | 4.98×107 | s−1 | [9] |
β61,0 | Rate constant for the transition 6→1 at Δψ=0 | 100 | s−1 | [9] |
ρF1 | Proton pump concentration | 0.7 | nmol | [4] |
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] |
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] |
Parameter | Description | Value | Unit | Reference |
fm | Fraction of free matrix calcium concentration | 3 × 10−4 | dimensionless | [4] |
Cm | Membrane capacitance | 1.450 × 10−3 | nmol mV−1 mgprot−1 | [4] |
fc | Fraction of free cytosolic calcium concentration | 1 × 10−2 | 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] |
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 |
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 |
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] |
Parameter | Description | Value | Unit | Reference |
α12 | Unimolecular rate constant for the transition 1→2 | 400 | s−1 | [9] |
α21 | Unimolecular rate constant for the transition 2→1 | 5 | s−1 | [9] |
α∗23 | Bimolecular rate constant for the transition 2→3 | 3.405×104 | nmol−1/2 mgprot1/2 s−1 | [4] |
α32 | Unimolecular rate constant for the transition 3→2 | 8×104 | s−1 | [9] |
α34 | Unimolecular rate constant for the transition 3→4 | 400 | s−1 | [9] |
α∗43 | Bimolecular rate constant for the transition 4→3 | 1.59×102 | nmol−1/2 mgprot1/2 s−1 | [4] |
α45 | Unimolecular rate constant for the transition 4→5 | 40 | s−1 | [9] |
α54 | Unimolecular rate constant for the transition 5→4 | 0.4 | s−1 | [9] |
α56 | Unimolecular rate constant for the transition 5→6 | 100 | s−1 | [9] |
α65 | Unimolecular rate constant for the transition 6→5 | 1×105 | s−1 | [9] |
α16,0 | Rate constant for the transition 1→6 at Δψ=0 | 130 | s−1 | [9] |
α61,0 | Rate constant for the transition 6→1 at Δψ=0 | 1×1012 | s−1 | [9] |
ρresp | Proton pump concentration | 0.4 | nmol mgprot−1 | [4] |
ΔψB | Boundary potential parameter | 50 | mV | [9] |
Parameter | Description | Value | Unit | Reference |
β12 | Unimolecular rate constant for the transition 1→2 | 400 | s−1 | [9] |
β21 | Unimolecular rate constant for the transition 2→1 | 40 | s−1 | [9] |
β∗23 | Bimolecular rate constant for the transition 2→3 | 5×103 | mgprot nmol−1 s−1 | [9] |
β32 | Unimolecular rate constant for the transition 3→2 | 5×103 | s−1 | [9] |
β34 | Unimolecular rate constant for the transition 3→4 | 100 | s−1 | [9] |
β∗43 | Bimolecular rate constant for the transition 4→3 | 5×104 | M−1 mgprot nmol−1 s−1 | [4] |
β45 | Unimolecular rate constant for the transition 4→5 | 100 | s−1 | [9] |
β54 | Unimolecular rate constant for the transition 5→4 | 100 | s−1 | [9] |
β56 | Unimolecular rate constant for the transition 5→6 | 1000 | s−1 | [9] |
β65 | Unimolecular rate constant for the transition 6→5 | 1000 | s−1 | [9] |
β16,0 | Rate constant for the transition 1→6 at Δψ=0 | 4.98×107 | s−1 | [9] |
β61,0 | Rate constant for the transition 6→1 at Δψ=0 | 100 | s−1 | [9] |
ρF1 | Proton pump concentration | 0.7 | nmol | [4] |
Parameter | Description | Value | Unit | Reference |
Jmax,ANT | Maximum exchange rate | 1000 | nmol mgprot−1 min−1 | [4] |
f | Fraction Δψ | 0.5 | dimensionless | [4] |
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] |
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] |