Jason N. Bazil

Modeling Mitochondrial Bioenergetics with Integrated Volume Dynamics

Mathematical models of mitochondrial bioenergetics provide powerful analytical tools to help interpret experimental data and facilitate experimental design for elucidating the supporting biochemical and physical processes. As a next step towards constructing a complete physiologically faithful mitochondrial bioenergetics model, a mathematical model was developed targeting the cardiac mitochondrial bioenergetic based upon previous efforts, and corroborated using both transient and steady state data. The model consists of several modified rate functions of mitochondrial bioenergetics, integrated calcium dynamics and a detailed description of the K+-cycle and its effect on mitochondrial bioenergetics and matrix volume regulation. Model simulations were used to fit 42 adjustable parameters to four independent experimental data sets consisting of 32 data curves. During the model development, a certain network topology had to be in place and some assumptions about uncertain or unobserved experimental factors and conditions were explicitly constrained in order to faithfully reproduce all the data sets. These realizations are discussed, and their necessity helps contribute to the collective understanding of the mitochondrial bioenergetics.

A Global Parallel Model Based Design of Experiments Method to Minimize Model Output Uncertainty

Model-based experiment design specifies the data to be collected that will most effectively characterize the biological system under study. Existing model-based design of experiment algorithms have primarily relied on Fisher Information Matrix-based methods to choose the best experiment in a sequential manner. However, these are largely local methods that require an initial estimate of the parameter values, which are often highly uncertain, particularly when data is limited. In this paper, we provide an approach to specify an informative sequence of multiple design points (parallel design) that will constrain the dynamical uncertainty of the biological system responses to within experimentally detectable limits as specified by the estimated experimental noise. The method is based upon computationally efficient sparse grids and requires only a bounded uncertain parameter space; it does not rely upon initial parameter estimates. The design sequence emerges through the use of scenario trees with experimental design points chosen to minimize the uncertainty in the predicted dynamics of the measurable responses of the system. The algorithm was illustrated herein using a T cell activation model for three problems that ranged in dimension from 2D to 19D. The results demonstrate that it is possible to extract useful information from a mathematical model where traditional model-based design of experiments approaches most certainly fail. The experiments designed via this method fully constrain the model output dynamics to within experimentally resolvable limits. The method is effective for highly uncertain biological systems characterized by deterministic mathematical models with limited data sets. Also, it is highly modular and can be modified to include a variety of methodologies such as input design and model discrimination.

A Minimal Model for the Mitochondrial Rapid Mode of Ca2+ Uptake Mechanism

Mitochondria possess a remarkable ability to rapidly accumulate and sequester Ca2+. One of the mechanisms responsible for this ability is believed to be the rapid mode (RaM) of Ca2+ uptake. Despite the existence of many models of mitochondrial Ca2+ dynamics, very few consider RaM as a potential mechanism that regulates mitochondrial Ca2+ dynamics. To fill this gap, a novel mathematical model of the RaM mechanism is developed herein. The model is able to simulate the available experimental data of rapid Ca2+ uptake in isolated mitochondria from both chicken heart and rat liver tissues with good fidelity. The mechanism is based on Ca2+ binding to an external trigger site(s) and initiating a brief transient of high Ca2+ conductivity. It then quickly switches to an inhibited, zero-conductive state until the external Ca2+ level is dropped below a critical value (∼100–150 nM). RaMs Ca2+- and time-dependent properties make it a unique Ca2+ transporter that may be an important means by which mitochondria take up Ca2+ in situ and help enable mitochondria to decode cytosolic Ca2+ signals. Integrating the developed RaM model into existing models of mitochondrial Ca2+ dynamics will help elucidate the physiological role that this unique mechanism plays in mitochondrial Ca2+-homeostasis and bioenergetics.

Dynamic buffering of mitochondrial Ca2+ during Ca2+ uptake and Na+-induced Ca2+ release

In cardiac mitochondria, matrix free Ca2+ ([Ca2+]m) is primarily regulated by Ca2+ uptake and release via the Ca2+ uniporter (CU) and Na+/Ca2+ exchanger (NCE) as well as by Ca2+ buffering. Although experimental and computational studies on the CU and NCE dynamics exist, it is not well understood how matrix Ca2+ buffering affects these dynamics under various Ca2+ uptake and release conditions, and whether this influences the stoichiometry of the NCE. To elucidate the role of matrix Ca2+ buffering on the uptake and release of Ca2+, we monitored Ca2+ dynamics in isolated mitochondria by measuring both the extra-matrix free [Ca2+] ([Ca2+]e) and [Ca2+]m. A detailed protocol was developed and freshly isolated mitochondria from guinea pig hearts were exposed to five different [CaCl2] followed by ruthenium red and six different [NaCl]. By using the fluorescent probe indo-1, [Ca2+]e and [Ca2+]m were spectrofluorometrically quantified, and the stoichiometry of the NCE was determined. In addition, we measured NADH, membrane potential, matrix volume and matrix pH to monitor Ca2+-induced changes in mitochondrial bioenergetics. Our [Ca2+]e and [Ca2+]m measurements demonstrate that Ca2+ uptake and release do not show reciprocal Ca2+ dynamics in the extra-matrix and matrix compartments. This salient finding is likely caused by a dynamic Ca2+ buffering system in the matrix compartment. The Na+- induced Ca2+ release demonstrates an electrogenic exchange via the NCE by excluding an electroneutral exchange. Mitochondrial bioenergetics were only transiently affected by Ca2+ uptake in the presence of large amounts of CaCl2, but not by Na+- induced Ca2+ release.

Catalytic Coupling of Oxidative Phosphorylation, ATP Demand, and Reactive Oxygen Species Generation

Competing models of mitochondrial energy metabolism in the heart are highly disputed. In addition, the mechanisms of reactive oxygen species (ROS) production and scavenging are not well understood. To deepen our understanding of these processes, a computer model was developed to integrate the biophysical processes of oxidative phosphorylation and ROS generation. The model was calibrated with experimental data obtained from isolated rat heart mitochondria subjected to physiological conditions and workloads. Model simulations show that changes in the quinone pool redox state are responsible for the apparent inorganic phosphate activation of complex III. Model simulations predict that complex III is responsible for more ROS production during physiological working conditions relative to complex I. However, this relationship is reversed under pathological conditions. Finally, model analysis reveals how a highly reduced quinone pool caused by elevated levels of succinate is likely responsible for the burst of ROS seen during reperfusion after ischemia.

A biophysically based mathematical model for the catalytic mechanism of glutathione reductase.

Glutathione reductase (GR) catalyzes the reduction of oxidized glutathione (GSSG) to reduced glutathione (GSH) using NADPH as the reducing cofactor, and thereby maintains a constant GSH level in the system. GSH scavenges superoxide (O2(*-)) and hydroxyl radicals (OH) nonenzymatically or by serving as an electron donor to several enzymes involved in reactive oxygen species (ROS) detoxification. In either case, GSH oxidizes to GSSG and is subsequently regenerated by the catalytic action of GR. Although the GR kinetic mechanism has been extensively studied under various experimental conditions with variable substrates and products, the catalytic mechanism has not been studied in terms of a mechanistic model that accounts for the effects of the substrates and products on the reaction kinetics. The aim of this study is therefore to develop a comprehensive mathematical model for the catalytic mechanism of GR. We use available experimental data on GR kinetics from various species/sources to develop the mathematical model and estimate the associated model parameters. The model simulations are consistent with the experimental observation that GR operates via both ping-pong and sequential branching mechanisms based on relevant concentrations of its reaction substrate GSSG. Furthermore, we show the observed pH-dependent substrate inhibition of GR activity by GSSG and bimodal behavior of GR activity with pH. The model presents a unique opportunity to understand the effects of products on the kinetics of GR. The model simulations show that under physiological conditions, where both substrates and products are present, the flux distribution depends on the concentrations of both GSSG and NADP(+), with ping-pong flux operating at low levels and sequential flux dominating at higher levels. The kinetic model of GR may serve as a key module for the development of integrated models for ROS-scavenging systems to understand protection of cells under normal and oxidative stress conditions.

A parallel algorithm for reverse engineering of biological networks

Dynamic biological systems, such as gene regulatory networks (GRNs) and protein signaling networks, are often represented as systems of ordinary differential equations. Such equations can be utilized in reverse engineering these biological networks, specifically since identifying these networks is challenging due to the cost of the necessary experiments growing with at least the square of the size of the system. Moreover, the number of possible models, proportional to the number of directed graphs connecting nodes representing the variables in the system, suffers from combinatorial explosion as the size of the system grows. Therefore, exhaustive searches for systems of nontrivial complexity are not feasible. Here we describe a practical and scalable algorithm for determining candidate network interactions based on decomposing an N-dimensional system into N one-dimensional problems. The algorithm was tested on in silico networks based on known biological GRNs. The computational complexity of the network identification is shown to increase as N(2) while a parallel implementation achieves essentially linear speedup with the increasing number of processing cores. For each in silico network tested, the algorithm successfully predicts a candidate network that reproduces the network dynamics. This approach dramatically reduces the computational demand required for reverse engineering GRNs and produces a wealth of exploitable information in the process. Moreover, the candidate network topologies returned by the algorithm can be used to design future experiments aimed at gathering informative data capable of further resolving the true network topology.

Modeling the calcium sequestration system in isolated guinea pig cardiac mitochondria.

Under high Ca2+ load conditions, Ca2+ concentrations in the extra-mitochondrial and mitochondrial compartments do not display reciprocal dynamics. This is due to a paradoxical increase in the mitochondrial Ca2+ buffering power as the Ca2+ load increases. Here we develop and characterize a mechanism of the mitochondrial Ca2+ sequestration system using an experimental data set from isolated guinea pig cardiac mitochondria. The proposed mechanism elucidates this phenomenon and others in a mathematical framework and is integrated into a previously corroborated model of oxidative phosphorylation including the Na+/Ca2+ cycle. The integrated model reproduces the Ca2+ dynamics observed in both compartments of the isolated mitochondria respiring on pyruvate after a bolus of CaCl2 followed by ruthenium red and a bolus of NaCl. The model reveals why changes in mitochondrial Ca2+ concentration of Ca2+ loaded mitochondria appear significantly mitigated relative to the corresponding extra-mitochondrial Ca2+ concentration changes after Ca2+ efflux is initiated. The integrated model was corroborated by simulating the set-point phenomenon. The computational results support the conclusion that the Ca2+ sequestration system is composed of at least two classes of Ca2+ buffers. The first class represents prototypical Ca2+ buffering, and the second class encompasses the complex binding events associated with the formation of amorphous calcium phosphate. With the Ca2+ sequestration system in mitochondria more precisely defined, computer simulations can aid in the development of innovative therapeutics aimed at addressing the myriad of complications that arise due to mitochondrial Ca2+ overload.

Analysis of the Kinetics and Bistability of Ubiquinol:Cytochrome c Oxidoreductase

Ubiquinol:cytochrome c oxidoreductase, bc1 complex, is the enzyme in the respiratory chain of mitochondria responsible for the transfer reducing potential from ubiquinol to cytochrome c coupled to the movement of charge against the electrostatic potential across the mitochondrial inner membrane. The complex is also implicated in the generation of reactive oxygen species under certain conditions and is thus a contributor to cellular oxidative stress. Here, a biophysically detailed, thermodynamically consistent model of the bc1 complex for mammalian mitochondria is developed. The model incorporates the major redox centers near the Qo- and Qi-site of the enzyme, includes the pH-dependent redox reactions, accounts for the effect of the proton-motive force of the reaction rate, and simulates superoxide production at the Qo-site. The model consists of six distinct states characterized by the mobile electron distribution in the enzyme. Within each state, substates that correspond to various electron localizations exist in a rapid equilibrium distribution. The steady-state equation for the six-state system is parameterized using five independent data sets and validated in comparison to additional experimental data. Model analysis suggests that the pH-dependence on turnover is primarily due to the pKa values of cytochrome bH and Rieske iron sulfur protein. A previously proposed kinetic scheme at the Qi-site where ubiquinone binds to only the reduced enzyme and ubiquinol binds to only the oxidized enzyme is shown to be thermodynamically infeasible. Moreover, the model is able to reproduce the bistability phenomenon where at a given overall flux through the enzyme, different rates of superoxide production are attained when the enzyme is differentially reduced.

A mechanistic mathematical model for the catalytic action of glutathione peroxidase.

Abstract Glutathione peroxidase (GPx) is a well-known seleno-enzyme that protects cells from oxidative stress (e.g., lipid peroxidation and oxidation of other cellular proteins and macromolecules), by catalyzing the reduction of harmful peroxides (e.g., hydrogen peroxide: H2O2) with reduced glutathione (GSH). However, the catalytic mechanism of GPx kinetics is not well characterized in terms of a mathematical model. We developed here a mechanistic mathematical model of GPx kinetics by considering a unified catalytic scheme and estimated the unknown model parameters based on different experimental data from the literature on the kinetics of the enzyme. The model predictions are consistent with the consensus that GPx operates via a ping-pong mechanism. The unified catalytic scheme proposed here for GPx kinetics clarifies various anomalies, such as what are the individual steps in the catalytic scheme by estimating their associated rate constant values and a plausible rationale for the contradicting experimental results. The developed model presents a unique opportunity to understand the effects of pH and product GSSG on the GPx activity under both physiological and pathophysiological conditions. Although model parameters related to the product GSSG were not identifiable due to lack of product-inhibition data, the preliminary model simulations with the assumed range of parameters show that the inhibition by the product GSSG is negligible, consistent with what is known in the literature. In addition, the model is able to simulate the bi-modal behavior of the GPx activity with respect to pH with the pH-range for maximal GPx activity decreasing significantly as the GSH levels decrease and H2O2 levels increase (characteristics of oxidative stress). The model provides a key component for an integrated model of H2O2 balance under normal and oxidative stress conditions.