Integrated computational model of the bioenergetics of isolated lung mitochondria

Integrated computational modeling provides a mechanistic and quantitative framework for describing lung mitochondrial bioenergetics. Thus, the objective of this study was to develop and validate a thermodynamically-constrained integrated computational model of the bioenergetics of isolated lung mitochondria. The model incorporates the major biochemical reactions and transport processes in lung mitochondria. A general framework was developed to model those biochemical reactions and transport processes. Intrinsic model parameters such as binding constants were estimated using previously published isolated enzymes and transporters kinetic data. Extrinsic model parameters such as maximal reaction and transport velocities were estimated by fitting the integrated bioenergetics model to published and new tricarboxylic acid cycle and respirometry data measured in isolated rat lung mitochondria. The integrated model was then validated by assessing its ability to predict experimental data not used for the estimation of the extrinsic model parameters. For example, the model was able to predict reasonably well the substrate and temperature dependency of mitochondrial oxygen consumption, kinetics of NADH redox status, and the kinetics of mitochondrial accumulation of the cationic dye rhodamine 123, driven by mitochondrial membrane potential, under different respiratory states. The latter required the coupling of the integrated bioenergetics model to a pharmacokinetic model for the mitochondrial uptake of rhodamine 123 from buffer. The integrated bioenergetics model provides a mechanistic and quantitative framework for 1) integrating experimental data from isolated lung mitochondria under diverse experimental conditions, and 2) assessing the impact of a change in one or more mitochondrial processes on overall lung mitochondrial bioenergetics. In addition, the model provides important insights into the bioenergetics and respiration of lung mitochondria and how they differ from those of mitochondria from other organs. To the best of our knowledge, this model is the first for the bioenergetics of isolated lung mitochondria.

Integrated computational modeling provides a mechanistic and quantitative framework for describing lung mitochondrial bioenergetics. Thus, the objective of this study was to develop and validate a thermodynamically-constrained integrated computational model of the bioenergetics of isolated lung mitochondria. The model incorporates the major biochemical reactions and transport processes in lung mitochondria. A general framework was developed to model those biochemical reactions and transport processes. Intrinsic model parameters such as binding constants were estimated using previously published isolated enzymes and transporters kinetic data. Extrinsic model parameters such as maximal reaction and transport velocities were estimated by fitting the integrated bioenergetics model to published and new tricarboxylic acid cycle and respirometry data measured in isolated rat lung mitochondria. The integrated model was then validated by assessing its ability to predict experimental data not used for the estimation of the extrinsic model parameters. For example, the model was able to predict reasonably well the substrate and temperature dependency of mitochondrial oxygen consumption, kinetics of NADH redox status, and the kinetics of mitochondrial accumulation of the cationic dye rhodamine 123, driven by mitochondrial membrane potential, under different respiratory states. The latter required the coupling of the integrated bioenergetics model to a pharmacokinetic model for the mitochondrial uptake of rhodamine 123 from buffer. The integrated bioenergetics model provides a mechanistic and quantitative framework for 1) integrating experimental data from isolated lung mitochondria under diverse experimental conditions, and 2) assessing the impact of a change in one or more mitochondrial processes on overall lung mitochondrial bioenergetics. In addition, the model provides important insights into the bioenergetics and respiration of lung mitochondria and how they differ from those of mitochondria from other organs. To the best of our knowledge, this model is the first for the bioenergetics of isolated lung mitochondria. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Introduction
In addition to the vital function of ATP production, mitochondria are involved in other important cellular functions, including apoptosis, calcium homeostasis, oxygen sensing, and nitric oxide signaling [1][2][3][4][5]. Mitochondria are also an important source of reactive oxygen species (ROS) under physiological and pathophysiological conditions, and are a primary target of oxidative stress [1,2,6]. Under normal conditions, mitochondria-derived ROS are involved in important signaling pathways, but when produced in excess, they contribute to oxidative stress, which is a key factor in the pathogenesis of lung diseases [2,7].
Mitochondrial respiration accounts for 80-85% of total lung ATP, and glucose is the most prevalent oxidizable substrate in lungs under physiological conditions [8]. There is ample evidence that mitochondrial dysfunction plays a key role in the pathogenesis of lung diseases, including acute lung injury (ALI), which is one of the most common causes of admission to intensive care units [2]. Using animal models of ALI and its most severe form, acute respiratory distress syndrome (ARDS), previous studies have reported significant changes in different mitochondrial processes, including decreases in the activity of complexes I and II, decreases in TCA cycle enzyme activities, dissipation of membrane potential, and impairment of ATP generation [9][10][11][12][13][14]. In addition, alteration in lung mitochondrial bioenergetics has been linked to pulmonary edema formation, which is a cardinal feature of ALI/ARDS [8]. Bongard et al. demonstrated that the integrity of the pulmonary endothelial barrier is dependent on normal lung mitochondrial bioenergetics [10].
There is a wealth of information regarding the impact of lung injury on mitochondrial enzymes, transporters and macromolecules [10,15]. However, because of the interdependence of mitochondrial processes, measurement of a change in one or more mitochondrial processes in response to injury is not sufficient to predict the functional impact of this change on overall mitochondrial bioenergetics. Integrated computational modeling provides a mechanistic and quantitative framework for describing mitochondrial bioenergetics, for quantifying the impact of a change in one or more processes on overall mitochondrial bioenergetics, and for identifying potential mitochondrial targets for diagnostic, prognostic, and therapeutic purposes. This type of modeling has been used to describe mitochondrial bioenergetics of other organs, including heart and skeletal muscle under physiological and pathophysiological conditions [16][17][18]. However, to the best of our knowledge, no such model exists for bioenergetics of isolated lung mitochondria.
Therefore, the objective of this study was to use new and published experimental data to develop and validate a thermodynamically-constrained integrated computational model of the bioenergetics of isolated lung mitochondria. Key results of this study are 1) development and validation of the first integrated computational model of the bioenergetics of isolated lung mitochondria, 2) coupling of this integrated model to a pharmacokinetics model of the mitochondrial uptake and release of the cationic membrane potential probe rhodamine 123 during different respiratory states for quantifying mitochondrial membrane potential, and 3) model simulations that provide important insights into the bioenergetics and respiration of isolated lung mitochondria and how they differ from those of mitochondria isolated from other organs.
Care and Use Committees of the Zablocki Veterans Affairs Medical Center, the Medical College of Wisconsin, and Marquette University.

Isolation of lung mitochondria
Adult male Sprague-Dawley rats (320-360g) were anesthetized (pentobarbital sodium 50-100 mg/kg, ip) and the lungs were rapidly exposed and cleared of residual blood with 50 ml cold perfusion solution (physiologic saline buffered with 10 mM HEPES, pH 7.4, and containing 5.5 mM glucose) via the right ventricle. The lungs were then removed from the chest, and the trachea, large airways and large vessels were removed, after which the peripheral lung was placed in an ice-cold homogenization buffer (pH 7.4) containing 10 mM HEPES, 200 mM mannitol, 70 mM sucrose, 1 mM EGTA, 2% fatty acid-free BSA, and protease inhibitor cocktail Set III (50 μl/g lung tissue; Calbiochem, San Diego, CA) and minced over ice. Lung tissue was then homogenized using a Teflon pestle. The resulting homogenate was then centrifuged (Sorvall Superspeed RC-5B, Norwalk, CT) at 2,000 × g and 4˚C for 15 minutes. The supernatant was transferred to a clean tube and centrifuged at 17,800 × g at 4˚C for 15 minutes. The resulting supernatant was discarded and the remaining pellet was resuspended in 5 ml ice-cold homogenization solution and centrifuged at 17,800 × g at 4˚C for 15 minutes. The supernatant was discarded and the final pellet was resuspended in 0.3-4 ml ice-cold buffer (same as the homogenization buffer without BSA or the protease inhibitor cocktail) and stored on ice to be used for the respirometry and membrane potential studies described below. Mitochondrial protein was determined using the Pierce BCA protein assay with bovine serum albumin as the standard. Mitochondrial yield was~2 mg/rat lung.

Mitochondrial oxygen consumption
Mitochondrial oxygen consumption (respiration) was measured polarographically at 23˚C with a Strathkelvin 782 2-Channel Oxygen System (Strathkelvin Instrument). Briefly, 0.55 ml respiration buffer (130 mM KCL, 5 mM K 2 HPO 4 , 20 mM MOPS, 0.1 mM EDTA, 0.001 mM Na 4 P 2 O 7 and 0.1% BSA, pH 7.2) was transferred to the reaction chamber. After 2 minutes, 0.55 mg mitochondrial protein (1 mg/ml) was transferred to the 0.55 ml reaction chamber. Stock solution containing respiratory substrates (pyruvate (10 mM) + malate (5 mM), glutamate (10 mM) + malate (5 mM), or succinate (7 mM) at pH 7.2) was introduced into the reaction chamber using a Hamilton syringe. Once the oxygen consumption rate reached steady state (state 2), ADP (100 µM or 50 µM) was added to stimulate the rate of oxygen utilization (state 3). State 4 was measured as the rate of oxygen consumption following depletion of ADP in the reaction chamber. The respiratory control index (RCI) for pyruvate + malate, malate + glutamate, or succinate was then calculated as the ratio of state 3 and state 4 respiratory rates. The rate of uncoupled oxygen consumption was measured in the presence of the uncoupler carbonyl cyanide-4-(trifluoromethoxy) phenylhydrazone (FCCP, 2 µM).
For a subset of experiments, maximal complex IV activity was measured at 30˚C using the complex IV substrate tetramethyl-p-phenylenediamine (TMPD, 0.3 mM) as a direct artificial electron donor in the presence of antimycin A (1 µM, complex III inhibitor). Ascorbate (1.3 mM) was also added to the medium to maintain TMPD in its reduced form. Furthermore, the integrity of the mitochondrial outer-membrane was assessed by evaluating the ability of exogenous cytochrome c (10 µM) added to the medium to stimulate oxygen consumption.
Mitochondrial membrane potential. As previously described [19], mitochondrial membrane potential (ΔC m ) studies were performed at room temperature (23˚C) using a Photon Technology International (PTI) QuantaMaster spectrofluorometer (HORIBA Scientific, Edison, New Jersey) that monitored and recorded the R123 fluorescent signal (503/527 nm converted to R123 concentration by assuming a linear relationship between R123 medium concentration and fluorescent signal.

Model development.
In what follows, all figure and table numbers with prefix "A" are found in S1 Supporting Information. The integrated model of the bioenergetics of isolated lung mitochondria was developed based on the thermodynamic and kinetic properties of the mitochondrial metabolic reactions and transporters [17,18,21]. As shown in Fig 1, the model consists of three regions: the extra-mitochondrial region (buffer), the mitochondrial matrix region, and the inter-membrane space (IMS) region. Moreover, the model accounts for the dynamics of forty state variables, including thirty-seven metabolite concentrations in the matrix and buffer regions, two ion concentrations (matrix H + and buffer H + concentrations), and mitochondrial membrane potential (ΔC m ). Volumes of the three regions along with initial concentrations of major metabolites are given in Table 1.
Fifteen reaction fluxes are included in this integrated model, including ten reactions in the tricarboxylic acid cycle (TCA or Krebs cycle) and five reactions in the respiratory system (complex I-complex V). Unlike liver mitochondria and kidney mitochondria, which can produce oxaloacetate (OXA) from pyruvate, pyruvate carboxylase activity is extremely low in lung mitochondria [27]. Thus, pyruvate carboxylase is not included in this model. In isolated mitochondrial experiments, bovine serum albumin (BSA) was used to bind free fatty acids and their esters [28]. Therefore, fatty acid oxidation was not considered in this isolated mitochondrial model. All the reactions in the integrated bioenergetics model are listed in Table A2 in S1 Supporting Information. Additionally, the model includes ten transport processes to account for the exchange of key metabolic species between the mitochondrial matrix and IMS regions. The mitochondrial inner membrane is impermeable to most metabolites. Therefore, except for proton leak, the exchange of species between the mitochondrial matrix and IMS regions is assumed to be catalyzed by specific transporters, such as inorganic phosphate carrier (PIC) and adenine nucleotide translocase (ANT). The proton leak between the IMS and mitochondrial matrix regions is modeled using modified Goldman-Hodgkin-Katz equation [18]. All the transporters in the model are listed in Table A3 in S1 Supporting Information.
The mitochondrial outer membrane has porins that allow free transport of small molecules [18]. Therefore, we assumed that the concentration differences of most of the metabolic species between the IMS and extra-mitochondrial space to be negligible. However, cytochrome c,

Parameter Value Source
Mitochondria matrix volume (V m ) 1 µl/mg mitochondria protein [18,22,23] Mitochondria IMS volume (V ims ) Assume to be 10% of V m Total cytochrome c content 0.33 nmol/mg mitochondria [24] Total pyridine nucleotide content (NAD + NADH) 1.73 nmol/mg mitochondria [24] Total flavin adenine dinucleotide concentration (FAD + FADH2) 0.7 mM [25] Total adenine nucleotide content (ATP and ADP) 6.4 nmol/mg mitochondria [24] Total coenzyme A content (SCoA + ACoA + CoA) 0.93 nmol/mg mitochondria [24] Aspartate + glutamate concentration in mitochondria matrix 12 mM [26] Total ubiquinone concentration (UQ + UQH 2 ) 0.52 mM [23] which has a molecular weight of about 12 KD, and hence too large to pass through the mitochondrial outer membrane under normal conditions, was assumed to exist only in the IMS region. Reaction and transport flux expressions. In the proposed integrated bioenergetics model, all enzymatic reactions are assumed to follow a generalized random-ordered rapidequilibrium kinetic mechanism [29], which encompass all possible kinetic mechanisms as special cases. Thus, for each reaction, all pertinent substrates must bind to the enzyme together before the catalysis can take place. For a general multi-substrate and multi-product enzymatic reaction: where S i is i th substrate, P j is j th product, N s and N p are the number of substrates and products, respectively, and α i and β j are the corresponding stoichiometric coefficients. The general form of the reaction flux, J, accounting for the thermodynamic (Haldane) constraint is given by: where K Si and K Pj are binding constants corresponding to substrates and products, respectively; [S i ] and [P j ] are concentrations of substrate i and product j, respectively; V maxf is the maximum forward reaction rate; and K 0 eq is the apparent equilibrium constant for the reaction, which is the value of the equilibrium-state reaction quotient (i.e. ratio of product of product concentrations over product of substrate concentrations) at specified thermodynamic conditions (i.e. temperature, ionic strength and pH). Detailed derivations of this general reaction flux equation and its various special cases are included in S1 Supporting Information.
In the presence of cofactor pairs (e.g. NADH and NAD + , ATP and ADP, GTP and GDP, CoA and ACoA, or CoA and SCoA), the generalized reaction flux Eq (2) can be modified appropriately so as to not include any interactive cofactor product terms [29]. The assumption is that the substrate and product represented as a cofactor pair bind with a given enzyme at the same binding site, in which case the resulting reaction flux equation does not include the corresponding substrate and product multiplication terms in the denominator of Eq 2. Eq 2 can also be suitably modified to account for other reaction kinetic mechanisms (e.g. sequentialordered mechanisms, ping-pong mechanisms). The flux expressions of all the enzymatic reactions are included in the S1 Supporting Information.
In this integrated bioenergetics model, all the transport processes across the mitochondrial inner membrane, with the exception of the proton leak, are catalyzed by specific metabolite transporters involving two different species (co-transporters and anti-porters). For generality, a random-ordered rapid-equilibrium binding mechanism is assumed for all the transporters, i.e., the transporter can bind to the two species in an arbitrary order, similar to that described for a reversible two-substrate two-product enzymatic reaction. Thus, the general form of the transport flux expression is assumed to be the same as that for the reaction flux expression (i.e. Eq 2). The flux expressions for all the metabolite transporters are included in S1 Supporting Information.
Accounting for the effect of pH on equilibrium constants. The proposed integrated bioenergetics model accounts for the pH dependence of the apparent equilibrium constants for proton-releasing reactions (Eq 3) and for proton-consumption reactions (Eq 4) [17,23,30]: where K 00 eq is the reaction apparent equilibrium constant (K 0 eq ) at pH of 7, and Δ r G '0 , R and T are the standard Gibbs free energy of the reaction at pH = 7, gas constant, and temperature, respectively.
Accounting for the effect of temperature on enzymatic reaction and transport rates. To allow for model simulations at different temperatures, we account for the temperature effects on the maximal reaction and transport rates (V max 's and T max 's) using the Q 10 temperature coefficient [23].
where T is the temperature, R 2 /R 1 is the correction coefficient (i.e. the ratio of maximal reaction and transport velocities at temperature T 2 and T 1 ), and Q 10 is the temperature coefficient which was set to 2.5.
Governing ordinary differential equations for the integrated bioenergetics model. The governing ordinary differential equations describing the dynamic changes in the concentrations (C) of the forty biochemical species were derived based on the principle of mass balance. The change in the concentration of a given species within the mitochondria matrix region and extra-mitochondria (buffer) region are given by: where J mj is the j th reaction flux in the mitochondrial matrix and J m-e,j is the j th transport flux between mitochondria matrix region and extra-mitochondria buffer region, V m and V e are the volumes of the mitochondrial matrix region and the extra-mitochondria region, respectively. For the extra-mitochondria region which does not include chemical reactions, the right hand side of Eq 7 contains just transport fluxes. Detailed mass balance equations are included in S1 Supporting Information. The rate of change of mitochondrial membrane potential (ΔC m ) is given by Eq 8 [18]: where C imm is the capacitance of the inner mitochondrial membrane; J CI , J CIII , J CIV and J CV are the reaction fluxes of complex I, complex III, complex IV and complex V, respectively;J Hleak is the proton leak between the IMS and mitochondria matrix regions; and J ANT is the transport flux characterizing ATP/ADP exchange via the ANT. MATLAB function 'ode15s' (MathWorks Inc.) was used to solve the system of governing differential equations.

Methods used for estimating the values of the unknown model parameters
As described below, some model parameters were set at published values, while others were estimated from new and published data, including isolated enzyme and transporter kinetics (see S1 Supporting Information), TCA cycle kinetics [27], and respirometry [24,27].
Estimation of the intrinsic model parameters. The binding constants (K's) of various substrates and products for enzymes and transporters were either set to published values [18,31] or estimated using published kinetic data from isolated mitochondrial enzymes and transporters. Some of the isolated enzyme and transporter kinetic data were from heart, kidney or liver mitochondria due to the scarcity of such data from lungs. The assumption is that the binding constants are intrinsic model parameters, and hence their values should be organindependent. For a given enzyme or transporter, the MATLAB function 'fmincon', a nonlinear program solver, was used with a least-squares objective function to estimate the values of the K's from the published kinetic data for the enzyme or transporter by fitting a given enzyme or transporter model to pertinent kinetic data. Models of reaction and transport fluxes for various enzymatic reactions and metabolite transporters, parameter estimation results, as well as the kinetic data are included in S1 Supporting Information.
Estimation of the extrinsic model parameters. With the binding constants (K's) of the various enzymes and transporters known, a genetic algorithm (GA, MATLAB function 'ga') was then used to estimate the values of the extrinsic model parameters (i.e. V maxf s and T maxf s) of the integrated bioenergetics model that best fit new and published dynamic kinetic data obtained from isolated lung mitochondria, including TCA cycle kinetics data [27], and respirometry data at different respiratory states [24,27]. Thus, the resulting values of extrinsic model parameters are specific to isolated rat lung mitochondria. GA is a global parameter estimation algorithm that mimics the process of natural selection and can escape local minima [32]. For the integrated model extrinsic parameter estimation, the objective function ƒ used is: where x i,j and X i,j are the model solutions and the corresponding experimental data at the i th time point and j th data set, respectively. N is number of data points and M is the number of data sets used for the parameter estimation.

Pharmacokinetic model for the mitochondrial uptake and release of the cationic dye rhodamine 123 (R123)
To validate the integrated bioenergetics model, we assessed its ability to predict experimental data that were not used for estimating the values of the unknown model parameters. To that end, we measured mitochondrial membrane potential (ΔC m ) using the fluorescent cationic dye R123 in isolated rat lung mitochondria. The fluorescent intensity was converted to dye concentration and compared with model simulations under the same experimental conditions. To simulate R123 dye disposition in isolated mitochondria, the integrated bioenergetics model was coupled with a modified version of a pharmacokinetic model by Gan et al. for the cellular uptake and mitochondrial accumulation of R123 [33]. Two additional state variables were needed, namely R123 concentrations in the mitochondria matrix [R123] m and in the extra-mitochondrial buffer[R123] e . The process which accounted for the transport of R123 from the buffer to the matrix driven by an electrochemical gradient was modeled using the Goldman-Hodgkin-Katz equation [20].
where F is Faraday' s constant, z is R123 valence, R is the gas constant, T is temperature, and p is the permeability coefficient of the membrane which was set to 1.38 mol/(liter mitochondria)/s/M. The ordinary differential equations that describe the change in the concentrations of R123 in matrix and buffer are: where V APP m and V APP e are the apparent volumes of mitochondrial matrix region (V m ) and extra-mitochondria buffer region (V e ), respectively [33].

Estimated values of the unknown model parameters
Most of the intrinsic parameters of the integrated bioenergetics model, such as the binding constants (K's) of the enzymatic reactions and transport processes, were estimated based on previously published isolated enzyme or transporter kinetic data shown in Figs A5-A14 in S1 Supporting Information. For instance, the parameters in Table A4 in S1 Supporting Information were estimated by fitting Eq A16 in S1 Supporting Information to the data in Fig A5 in S1 Supporting Information using the MATLAB function "fmincon". The solid lines superimposed on the data in Figs A5-A14 in S1 Supporting Information are fits of corresponding enzyme/ transporter models to the data in S1 Supporting Information. Since the binding constants of the enzymatic reactions and transport processes are intrinsic model parameters, they are assumed to be tissue-independent [23]. It is worth noting that the extrinsic model parameter T maxf s and V maxf s values estimated from the isolated enzyme and transporter kinetic data were not used in the integrated bioenergetics model since those parameters are expected to be dependent on the tissue source of mitochondria.
With the K's of the various enzymes and transporters known, the next step was to estimate the extrinsic parameters (i.e. T maxf s and V maxf s) of the integrated bioenergetics model using kinetic data from isolated rat lung mitochondria. Multiple datasets form isolated rat lung mitochondria were used, including previously published TCA cycle metabolite uptake/release dynamic data (Fig 2) and respirometry data at 30˚C (Fig 3) by Evans et al. [27], and newly measured respirometry data at room temperature (23˚C, Fig 4). All of the extrinsic parameters (Table 2) of the integrated bioenergetics model were estimated by simultaneously fitting the integrated model solution to kinetic data from isolated rat lung mitochondria (Figs 2-4) using genetic algorithm (MATLAB function "ga"). The estimated values of the extrinsic parameters of the integrated bioenergetics model are shown in Table 2.
As shown in Fig 2, mitochondrial pyruvate uptake, malate uptake, and citrate production/ release were measured in isolated rat lung mitochondria by Evans et al. in the presence of 5 mM malate and varying pyruvate concentrations (top panel) or 5 mM pyruvate and varying malate concentrations (bottom panel) in buffer [27]. Results show that pyruvate uptake, malate uptake, and citrate release increased with increasing buffer pyruvate concentration at 5 mM of malate (top panel) and increasing malate concentration at 5 mM of pyruvate (bottom panel), and that those flux rates saturated when pyruvate and malate concentrations were greater than 5 mM.
In addition to the TCA cycle metabolite uptake/release data, the study by Evans et al. also provided mitochondrial oxygen consumption rates using different metabolic substrates [27]. As shown in Fig 3, mitochondria oxygen consumption rates under states 2 and 3 were measured in the presence of different complex I substrates (pyruvate + malate, α-ketoglutarate (AKG) + pyruvate, and citrate + pyruvate) and complex II (succinate) substrates, with the buffer temperature set at 30˚C [27]. The data show that states 2 and 3 oxygen consumption rates tended to be higher in the presence of complex II substrate as compared to those in the presence of complex I substrates. However, only oxygen consumption rates are provided in this study, whereas key information such as length of state 3 and oxygen dynamic changes under state 3 conditions are lacking. Thus for additional information on rat lung mitochondria electron transport chain (ETC) for the estimation of pertinent extrinsic model parameters, we measured oxygen consumption (Fig 4) at 23˚C using pyruvate + malate or glutamate + malate as complex I substrates or succinate as complex II substrate. The resulting average oxygen  Table 2 Fig 3 (Evans et al. data), the maximum rate of proton leak (T maxf, LEAK ) was adjusted to account for those differences. Therefore, for the data in Fig 4, T maxf, LEAK was set at 40% of the value used to fit the data in Fig 3 (Table 2). Computational modeling of lung mitochondrial bioenergetics However, for experimental data from the same study (Figs 2 and 3), T maxf, LEAK was set to the same value.

Measures of identifiability and estimatility of the extrinic parameters of the integrated model
To assess the identifiability and estimability of the extrinsic parameters of the integrated bioenergetics model (Table 2), we estimated the parameters' normalized sensitivity coefficients and a matrix of correlation coefficients between the model parameters. The normalized sensitivity coefficients provide information about the contribution of each of the extrinsic model parameters to the overall model solution, whereas the correlation coefficient matrix provides information about the degree of interdependence between the various model parameters. For a given extrinsic parameter, the normalized sensitivity coefficient was determined using Eq 12 [34]: where E is the sum of squared difference between experimental data (Figs 2-4) and the Computational modeling of lung mitochondrial bioenergetics integrated bioenergetics model fit, and y i is the estimated value of i th intrinsic parameter. @E @y i was approximated using the central difference method with 0.1% change in y i . The matrix of correlation coefficients between the model parameters was evaluated at the values in Table 2 that best fit the integrated model to the data (Figs 2-4). The correlation coefficient (CC ij ) between the ith parameter and jth parameter was determined using Eq 13 [35]: where np is the number of model parameters, HH is the inverse of the product of the transpose of the Jacobian matrix and the Jacobian matrix evaluated at the values of the model parameters in Table 2 that best fit the integrated bioenergetics model to the data in

Experimental and computational results and model validation
Oxygen consumption. Respirometry data measured at 23˚C are shown in Fig 4. Averaged time-course respirometry data (mean ± SE) are shown in the presence of complex I substrates pyruvate + malate (Fig 4A, n = 4), complex I substrates malate + glutamate (Fig 4B, n = 1), and complex II substrate succinate (Fig 4C, n = 4). A summary of the oxygen consumption rates (OCR, nmol/min/mg) show higher states 2, 3 and 4 respiratory rates with complex II substrate ( Fig 4F) as compare to those with complex I substrates (Fig 4D and 4E). Furthermore, the results show smaller respiratory rates, but higher respiratory control index (RCI, defined as the ratio of state 3 OCR to state 4 OCR) than those in Fig 3 (at 30˚C) for pyruvate + malate and succinate. Fig 4 also shows fits of the integrated bioenergetics model to the measured oxygen consumption kinetic data (solid lines, upper panels) and oxygen consumption rates (red bars, lower panels) generated using the values of the model parameters in Table 2, except for the value of the T maxf parameter descriptive of the proton leak process (T maxf, LEAK ), which was set at 40% of the value in Table 2. This is consistent with a less "leaky" and more coupled isolated rat lung mitochondria for the data in  Computational modeling of lung mitochondrial bioenergetics Complex IV activity and integrity of mitochondrial outer membrane. Maximal complex IV activity was measured by adding TMPD (artificial electron donor for complex IV) along with ascorbate to the medium in the presence of antimycin A. The measured oxygen consumption rates (OCR) are shown in Table 3, which also shows OCR in the presence of exogenous cytochrome C and succinate. The results in Table 3 show no significant difference between the oxygen consumption rates measured under state 3 conditions in the presence of  Values are mean ± SE. OCR, oxygen consumption rate; CytoC, cytrochrome c; TMPD, tetramethyl-pphenylenediamine; NA, not available. https://doi.org/10.1371/journal.pone.0197921.t003 Computational modeling of lung mitochondrial bioenergetics succinate and those measured in the presence of TMPD + ascorbate or TMPD + ascorbate + cytochrome C. These results are consistent with those reported by Fisher et al. [24,27]. Mitochondrial membrane potential. We quantified ADP-stimulated depolarization of ΔC m in mitochondria isolated from rat lungs using R123 in the presence of complex I (pyruvate + malate) or compex II (succinate) substrates at 23˚C. Fig 8A shows that in the presence of mitochondria and complex I substrates, the addition of ADP (state 3) stimulated a transient and reversible efflux of R123 from mitochondria, consistent with transient and reversible partial depolarization of ΔC m . The addition of 0.1 mM ADP resulted in a greater and longer membrane potential depolarization than the addition of 0.05 mM ADP. The addition of the mitochondrial uncoupler FCCP at the end of the protocol maximally depolarized ΔC m resulting in the maximal efflux of R123 from mitochondria [20]. Higher membrane potential (thus lower R123 concentration in buffer) was obtained with succinate as a complex II substrate using the same protocol. Fig 8A also shows the simulated (model predicted) R123 dye concentration in buffer using the integrated bionergetics model and the pharmacokinetic model of rhodamine 123 mitocondrial uptake with the bioenergetics model parameter values set to those in Table 2  Model prediction of ACoA and CoA steady-state data and dynamic NADH redox status data. In addition to membrane potential data (Fig 8), we evaluated the ability of the integrated model to predict steady-state ACoA and CoA fractions at different malate concentrations [27], and kinetic NADH redox status data under different respiratory states [24]. Fig 9A  and 9B show the ability of the integrated bioenergetics model to predict quit well both sets of data using the values of the model parameters in Table 2 Table 2. https://doi.org/10.1371/journal.pone.0197921.g008 Computational modeling of lung mitochondrial bioenergetics

Discussion
We have developed and validated the first integrated mechanistic computational model of the bioenergetics and respiration of isolated lung mitochondria using new and previously published experimental data. The model provides a mechanistic and quantitative framework for 1) integrating experimental data from isolated lung mitochondria under diverse experimental conditions, and 2) assessing the impact of a change in one or more mitochondrial processes on overall lung mitochondrial bioenergetics. In addition, as described below, this integrated model provides important insights into the bioenergetics and respiration of isolated lung mitochondria and how they differ from those of mitochondria from other organs.
Similar computational models of mitochondrial bioenergetics and respiration have been developed for mitochondria isolated from other organs, including the heart and skeletal muscles [16-18, 21, 23]. Such models differ from the current integrated model in many aspects, in part because of differences in mitochondrial bioenergetics and energy demand of lung tissue as compared to those of the heart and skeletal muscle. It is well known that mitochondrial content, major metabolic substrates, and control mechanisms vary widely between different organs. For instance, β-oxidation of fatty acids has been identified as the primary source of energy for heart under physiological conditions accounting for 70% of ATP production [36][37][38], whereas glucose is the primary source of energy for lung tissue and accounts for 80-85% of the total lung tissue ATP production under physiological conditions [8]. Moreover, lung cells are mostly non-excitable, and thus are postulated to behave very differently from excitable cells such as cardiac and skeletal muscle cells [39]. For instance, the ATP demand for cardiac cells between resting state and maximal exercise state could increase by~4 fold, while that for skeletal muscle cells could increase by as much as~100 fold between resting state and maximal exercise state [40,41]. Thus, key enzymes in existing computational models are modeled to be  Table 2. modulated by metabolic controllers, such as inorganic phosphate (Pi), ADP/ATP ratio, NADH/NAD + ratio, etc. However, in lung mitochondria where the energy requirement is much more stable, such control mechanisms may not be as important. Furthermore, most existing computational models place emphasis on complex I substrate-dependent respiration [18,23]. These differences were one of the reasons for the development and validation of an integrated model specific for the bioenergetics and respiration of isolated lung mitochondria.
Another reason for developing a comprehensive integrated model for bioenergetics and respiration of isolated lung mitochondria is the relatively complex flux equations for enzymes and transporters used in the existing heart and skeletal muscle models [16][17][18]. Those equations are not feasible for an isolated lung mitochondrial model due to the scarcity of experimental data needed to estimate the values of the large number of unknown parameters in the individual flux expressions and in the integrated mitochondrial model [27]. In the present study, a simplified general flux equation that accounts for apparent binding constants was developed and utilized for all enzymatic reactions and metabolite transporters. As such, the number of unknown model parameters was minimized for the individual flux expressions as well as for the integrated bioenergetics model. Yet, the model incorporated major pathways associated with lung mitochondrial bioenergetics, including metabolite transports and oxidations, TCA cycle reactions, electron transport chain (ETC) reactions, and oxidative phosphorylation.

Estimation of unknown model parameters, and integrated model predictions
For intrinsic model parameters, such as binding constants (K's), for major enzyme reactions and transport processes, their values were estimated using experimental data measured in isolated enzymes or fixed to published values, mostly from organs other than lungs due to lack of such data from lungs. The assumption is that for a given mitochondrial enzyme or transporter, intrinsic properties, such as binding constants K's, are organ-independent [18].
With intrinsic model parameters known, a challenging aspect of the integrated bioenergetics model development was estimation of the values of the extrinsic parameters ( Table 2). To that end, we relied on previously published experimental data from multiple studies in isolated rat lung mitochondria as well as new data that were obtained in the present study for further validation and corroboration of the model. As shown in Figs 2-4 and 8-9, the model is capable of fitting (Figs 2-4) and predicting (Figs 8-9) quite well all the data from isolated rat lung mitochondria with the same set of values for the extrinsic model parameters (Table 2).
A genetic algorithm (GA) was used to estimate the values of the extrinsic model parameters that best fit the kinetic data in Figs 2-4. Even though GA is a global parameter estimation algorithm that does not require initial estimates for the model parameters, good initial estimates can significantly reduce the computational time needed to identify optimal values of the model parameters. Good initial estimates can be obtained from boundary conditions such as pyruvate uptake rate, citrate production rate, and malate production rate.
For the TCA cycle experimental data measured by Evans et al. [27] and used for estimating the values of some of the model parameters, citrate formation was reported even when no exogenous pyruvate was added to the buffer medium (Fig 2C and 2F). In our integrated bioenergetics model, both pyruvate and malate are required to generate citrate. This resulted in a difference between model fits and experimental data when the pyruvate buffer concentration is zero. One possible explanation for this apparent inconsistency is that citrate measured at zero pyruvate buffer concentration, in the experiments by Evans et al., was from endogenous substrates, such as ACoA [27]. However, reported endogenous ACoA of~1 nmol/mg mitochondria [27] is too small to account for the measured 40 nmol citrate. Since citrate formation was observed only after 4 min of incubation time (as shown in Fig 2C and 2F), it is more likely that the source of citrate was other metabolic pathways that are not included in this model such as amino acids from mitochondrial proteins [42]. Another potential source of medium pyruvate is the presence of some endogenous pyruvate in the isolated mitochondria.
Identifiability and estimability of the extrinsic parameters of the integrated model Fig 5 shows that proton leak extrinsic parameter T maxf, LEAK has the largest normalized sensitivity coefficient among the 25 extrinsic model parameters (Table 2) of the integrated model. This could be because mitochondrial membrane leakiness not only affects oxygen consumption rate, RCI, and membrane potential, but also substrate consumption rate at state 2 since proton leak is the only pathway in this model for consumption of the energy provided by substrate oxidation under state 2 conditions (in the absence of ADP).
The intrinsic model parameters for glutamate-hydrogen or glutamate-hydroxyl antiporter (GLUH), nucleoside diphosphokinase (NDK), fumarate hydratase (FH), and α-ketoglutarate (2-oxoglutarate) malate exchanger (OME) reactions have the smallest normalized sensitivities coefficients (Fig 5). GLUH exists in rat liver mitochondria and heart mitochondria [43], but does not exist in rat brain mitochondria [43]. The existence of GLUH in rat lung mitochondria has not been confirmed, even though it was proposed that GLUH might be an important process for the replenishment of TCA cycle metabolites in the rat lung [44]. In the present study, the estimated GLUH activity and sensitivity are very low, indicating that GLUH is not required to fit the data in Figs 2-4 or to predict the data in Figs 8 and 9. On the other hand, glutamate oxaloacetate (GOT) is the major pathway for glutamate + malate-driven respiration. The GOT enzyme activity reported by Evans et al. [27,44] (975 nmol/min/mg protein) is both necessary and sufficient to reproduce the glutamate + malate respirometry data in Fig 4B. Since the NDK and FH reactions are running at near equilibrium, their corresponding extrinsic parameters are not identifiable without additional data.
We also estimated the matrix of correlation coefficients between the extrinsic parameters (Fig 6) of the integrated bioenergetics model for the parameter values in Table 2. The extrinsic parameter pairs with the highest positive correlation are those for PDH/PYRH and DCC (SUC)/SDH reactions. The extrinsic parameter pair with the highest negative correlation is that for TCC/CITDH reactions. Extrinsic parameters for equilibrium reactions such as NDH and FH have virtually no correlation with the extrinsic parameters of other reactions.

Membrane integrity of isolated mitochondria
Isolation of mitochondria from rat lung tissue is challenging due to the high lipid content and low mitochondria content [45]. For instance, cardiac myocyte mitochondria account for~35% of cell volume, compared to 1-2% in pulmonary endothelial cells, which account for 50% of the cells in the lung [45]. Thus, a relatively high concentration of bovine serum albumin is required in the isolation buffer to bind free fatty acid and lipids. In addition, the integrity of lung mitochondrial membrane in the isolated mitochondria is highly dependent on the isolation protocol used. For the data from Evans et al. (Fig 3), the measured state 2 oxygen consumption rate was higher than that measured in the present study (Fig 4). This could be due to differences in the mitochondria isolation protocols used. Thus for the model fit to the data in Figs 2-4, the value of the leak parameter (T maxf, LEAK , which accounts for state 2 respiratory rate) was adjusted to account for differences in the measured state 2 respiratory rate. Fig 7 show the effect of temperature and proton leak activity (T maxf, LEAK ) on mitochondrial oxygen consumption. With pyruvate + malate as metabolic substrates, both state 2 and state 3 oxygen consumption rates increase (Fig 7A), although the % increase in state 2 appears to be larger than that for state 3. As a result, there is an inverse relationship between proton leak activity and RCI (Fig 7B). In contrast, with succinate as the metabolic substrate, the state 3 oxygen consumption rate decreases as proton leak activity increases ( Fig  7E). Therefore, membrane leakiness has a higher impact on RCI when using succinate as substrate.

Simulations in
By altering the proton leak activity (T maxf, LEAK ), the integrated model was capable of simulating a wide range of mitochondrial experimental data performed under different experimental conditions (temperature, mitochondria density, pH, etc.) with the same set of values for the other model parameters (adjusted for temperature and pH effects). Fig 7C and 7F show the effect of temperature on dynamic responses of membrane potential and NADH redox states at state 3. Fig 7C also shows that ADP is consumed faster at higher temperature. As a result, the length of state 3 decreases as temperature increases.

Apparent functionally incomplete TCA cycle in isolated lung mitochondria
Isolated rat lung mitochondria and heart mitochondria show different behavior under similar experimental conditions. For instance, isolated rat lung mitochondria with pyruvate + malate as respiratory substrates produce significantly higher citrate than isolated cardiac mitochondria [26]. For lung mitochondria with pyruvate + malate as metabolic substrates, citrate production rate accounts for~80% of the pyruvate and malate uptake rate. On the other hand, for heart mitochondria, the consumed pyruvate and malate rate is approximately equal to amounts of citrate, α-ketoglutarate, succinate, and fumarate production rates [26].  Table 2.
The large citrate production rate in isolated rat lung mitochondria is consistent with an apparent functionally incomplete TCA cycle since most of the consumed carbon from pyruvate and malate is converted to citrate and released into the extra-mitochondrial buffer medium. Therefore, not enough substrate is left for the reaction to proceed beyond citrate in the TCA cycle. This might also be due to high NADH significantly inhibiting the CITDH reaction under physiological conditions, based on isolated enzyme experiments, as shown in Fig  7A in S1 Supporting Information (K NADH = 4.7 μM). Integrated bioenergetics model predicted simulations (Fig 10) also show that the reaction fluxes of CITDH and AKGDH are very low compared with those for PYRDH and MALDH.
Additional evidence of the apparent functionally incomplete TCA cycle in isolated rat lung mitochondria is the relatively high pyruvate/oxygen ratio since pyruvate consumption rate (10 nmol/min/mg mitochondria protein, as shown in Fig 2) is comparable to oxygen consumption rate (12 nmol/min/mg mitochondria protein, as shown in Fig 3) under the same experimental conditions with pyruvate + malate as metabolic substrates. If all the TCA cycle reactions are active, one would expect 4 NADH to be generated from each pyruvate (PYRDH, CITDH, MALDH and AKGDH), and 8 electrons to transfer to complex IV and reduce two molecules of oxygen. However, the data from Evans et al. (Figs 2 and 3) show that only 1.2 oxygen molecules are consumed for each pyruvate molecule. The high pyruvate/oxygen ratio in experiments suggests that only around 2.4 molecules of NADH are produced for each pyruvate consumed, consistent with only two of the four dehydrogenases being functional.
The observed large citrate production and the release to extra-mitochondrial buffer in isolated lung mitochondria with complex I substrate (pyruvate + malate) could be a means to ensure the availability of enough lipid synthesis for lung surfactant since citrate is an important source of acetyl units for ATP citrate lyase, an important step in fatty acid biosynthesis. The lung is an active organ that produces fatty acid to maintain lung surfactant balance, which is composed of approximately 90% lipids [46].
It is well known that calcium ions (Ca 2+ ) have a stimulating effect on key dehydrogenase enzymes in mitochondria [41]. Considering that not all of the dehydrogenase enzymes are active in isolated mitochondrial experiments (Fig 10), the effect of Ca 2+ control on the TCA cycle and bioenergetics might be underestimated in isolated mitochondrial experiments. For example, Vinnakota et al. [41] showed that physiological and supra-physiological Ca 2+ levels have no or only modest stimulatory effect on cardiac mitochondrial oxidative phosphorylation, respectively, in isolated mitochondrial respiring on complex I substrates (pyruvate + malate), which is insufficient to explain the Ca 2+ stimulating effect on oxidative phosphorylation in vivo in the heart. The apparent functionally incomplete TCA cycle in isolated lung mitochondria due to excess citrate production and release to buffer (Figs 2 and 10) could provide an explanation for why the stimulating effect of Ca 2+ on mitochondrial bioenergetics and respiration with pyruvate + malate as respiratory substrate is less evident in isolated cardiac mitochondria than in vivo. In addition, significant NADH inhibition of the CITDH reaction at physiological conditions could be another mechanism behind the non-stimulating effects of Ca 2+ on mitochondrial bioenergetics and respiration in isolated cardiac mitochondria with pyruvate + malate as respiratory substrates.

Rate-limiting steps in rat lung mitochondrial electron transport chain (ETC)
In the present study, experimental data and model sensitivity analysis reveal key information regarding the limiting step(s) in the rat lung mitochondrial ETC. Consistent with previously reported data [24,27], states 2 and 3 oxygen consumption rates in the presence of succinate were higher than those measured in the presence of pyruvate + malate.
In isolated heart mitochondria, simultaneous addition of complex I and complex II substrates increased state 3 oxygen consumption rate by up to 2-fold relative to that following the addition of complex I or complex II substrates alone [47]. In contrast to heart and liver mitochondria, our data show that convergent electron flow from both complex I substrates and complex II substrates is not additive (e.g. oxygen consumption rate with pyruvate + malate + succinate did not increase relative to that with succinate alone), indicating that in isolated rat lung mitochondria the rate limiting step(s) for the flow of electrons in the ETC is downstream from complexes I and II. Therefore, complexes III, IV, V and/or ANT could be the potential rate limiting step.
At the end of each isolated rat lung mitochondrial experiment in Fig 4, the uncoupler FCCP was added to the surrounding buffer to increase proton leak from the inter-membrane space to the mitochondria matrix. In agreement with a previous rat lung mitochondria study [24], our data (Fig 4 and Table 3) show that the oxygen consumption rate at uncoupled state 3 did not increase relative to that of ADP-stimulated state 3, suggesting that there is no excess ETC capacity over and above that for oxidative phosphorylation. Thus, electron flow under state 3 respiration is not limited by mitochondrial complex V or ANT. It is more likely that complex III and/or complex IV are the rate limiting steps in rat lung mitochondria.
Integrated model sensitivity analysis can reveal important information regarding rate limiting steps in mitochondrial ETC. A high sensitivity coefficient for a given parameter indicates that the model output is sensitive to the enzyme whose activity is described by that parameter. Therefore, an enzyme with a relatively high sensitivity coefficient is more likely to be the rate limiting step than an enzyme with a low sensitivity coefficient. Sensitivity analysis results (Fig 5) show that mitochondrial complex I and complex IV have much higher normalized sensitivity coefficients than complex III. One possible explanation is that with complex I substrates, there is excess electron transport capacity in ETC, and that mitochondrial electron flow is limited by complex I activity. However, when complex II substrate is present, electron flow from complexes I and II exceeds the maximum capacity of complex IV, and thus complex IV becomes the rate limiting step.
To prove that mitochondrial complex IV is the rate limiting step, maximum complex IV activity was measured using the complex IV substrate TMPD + ascorbate as direct electron donors. Results are compared with oxygen consumption rate with succinate as the metabolic substrate. As shown in Table 3, ADP-stimulated maximum respiration rate (state 3) in the presence of succinate was not different from that in the presence of TMPD + ascorbate. Similar results were reported by Fisher et al. [24,27] (Table 3). These results suggest that maximum complex IV activity is reached when succinate is used as the metabolic substrate, and that additional oxygen consumption cannot be stimulated by the addition of an uncoupler or complex I substrates to the surrounding medium.
Release of cytochrome c due to the relatively leaky mitochondrial outer-membrane could be a reason for the apparent low complex IV activity. If so, the addition of exogenous cytochrome c to the medium should markedly stimulate oxygen consumption. As shown in Table 3, adding exogenous cytochrome c had a small (12.5%) and insignificant effect (paired t-test, p = 0.11) on the measured rate in the absence of exogenous cytochrome c in the medium. Considering that a 10% increase of the cytochrome c stimulated oxygen consumption can be regarded as a sign of intact outer-membrane in heart and muscle mitochondria [48], a 12.5% increase in oxygen consumption may suggest a slightly "leaky" mitochondrial outer-membrane.

Pharmacokinetic model for R123 mitochondrial uptake from medium
A novel feature of the present study is the coupling of a pharmacokinetic model for R123 mitochondrial uptake to the integrated bioenergetics model to predict the measured R123 medium concentration and quantify mitochondrial membrane potential (Fig 8). R123 fluorescence intensity increases linearly with dye concentration at low concentrations (<1 µM) in the buffer [20]. Thus, we converted the R123 fluorescence intensity to R123 concentration to compare with model simulated R123 concentration in the buffer. Integrated model simulations in Fig 8 show good agreement with experimental data. In addition, integrated model predicted lung mitochondrial membrane potential (~140 mV) is close to the value (133 ± 4 (SE) mV) measured by Gan et al. in cultued bovine pulmonary artery endothelial cells [33] and to the value (129 ± 4 (SE) mV) meaured by Bongard et al. in cultured rat pulmonary microvascular endothelial cells [9].

Limitations and future directions
Metal ions such as calcium (Ca 2+ ) and magnesium (Mg 2+ ) ions play an important role in rat lung mitochondrial bioenergetics. In this integrated bioenergetics model, all the metal ion concentrations are assumed constant. Even though it is known that dehydrogenases are activated by Ca 2+ at nM level (~100-300 nM), Fisher et al. reported that mitochondrial oxygen consumption rates are inhibited as Ca 2+ concentration increases to µM-mM ranges [24]. However, the detailed mechanism for this biphasic effect is not known. Therefore, additional studies and data from isolated rat lung mitochondria are required to incorporate metal ions into the current model.
Another limitation of this model is that it does not differentiate between the forty different cell types in the lung. Endothelial cells, which account for~50% of all lung cells, are not rich in mitochondrial [33]. However, other lung cell types including contractile (smooth muscle cells), phagocytic (alveolar macrophages) and epithelial cells are rich in mitochondria. Thus, the model and the estimated values of model parameters should be interpreted as a mechanistic description of the average mitochondrial bioenergetics of all lung cells. Although the question regarding the mitochondrial bioenergetics and respiration of specific lung cell types is important, understanding the bioenergetics/respiration of mitochondria isolated from lung tissue and alteration in those bioenergetics/respiration is important and has functional implications regardless of the lung cell types involved.
The proposed integrated bioenergetics model is an important step towards the development and validation of a comprehensive thermodynamically-constrained computational model of isolated perfused rat lung tissue bioenergetics that includes cytosolic processes. Such a model will be used to evaluate the impact of ALI-or ARDS-induced changes in one or more mitochondrial or cytosolic processes on overall lung tissue bioenergetics, and for assessing the impact of targeting specific processes for mitigating the impact of ALI and ARDS on lung tissue bioenergetics.
Supporting information S1 Supporting Information. Reactions and transport processes included in the integrated bioenergetics model along with model governing mass balance equations. This file consists of four parts. Part A lists the reaction and transport processes in the mitochondrial bioenergetics model. Part B provides a description of the derivation of the generalized metabolic reaction and transport flux equations. Part C lists the flux expressions for specific metabolic reactions and transport processes. Part D lists the governing mass balance equations for the mitochondrial bioenergetics model. (DOCX)