A computational approach to determine the mechanisms of altered cardiac mitochondrial bioenergetics in a mouse model of genetic cardiomyopathy

Heart Failure (HF) is the leading cause of death in the U.S. and world wide. HF can arise from a variety of factors, some being environmental others being genetic. Chromosome 1p36 deletion is the most common human chromosomal deletion syndrome, and the associated loss of the PRDM16 gene often results in HF. PRDM16 deletion has been associated with left ventricular non-compaction and dilated cardiomyopathies. The importance of PRDM16 during heart development is being increasingly recognized, however the mechanisms by which PRDM16 loss leads to HF remain largely unknown. Numerous previous studies have found that altered cardiomyocyte mitochondrial function is a hallmark of many heart failure phenotypes. Therefore, a deeper understanding of the mechanistic relationship between cardiac mitochondrial bioenergetics and heart failure in the context of PRDM16 deletion may provide a target for therapies aimed at mitigating the degenerative changes associated with PRDM16 loss in 1p36 deletion syndrome. Exploring mechanistic changes of cardiac mitochondria is difficult to do experimentally so the aim of this project was to develop a computational model of cardiac mitochondria to use in combination with experimentation to elucidate the bioenergetic changes associated with the loss of PRDM16, in an effort to further understand what might be leading to heart failure in PRDM16 null hearts. Thermodynamically constrained enzymatic reaction and transport flux equations were derived based on principles of Michaelis Menten kinetics. These flux equations were used to derive governing ordinary differential equations, which were then used as the basis of the computational model. The parameters needed to fit these flux equations with the kinetics of specific enzymes and transporters were determined through a combination of literature values and experimental data. Once the model was parameterized, it was then used to quantify the mechanistic contribution of altered mitochondrial enzymes and transporters to cardiomyocyte bioenergetics in Prdm16-null hearts.

ABSTRACT Heart Failure (HF) is the leading cause of death in the U.S. and world wide. HF can arise from a variety of factors, some being environmental others being genetic. Chromosome 1p36 deletion is the most common human chromosomal deletion syndrome, and the associated loss of the PRDM16 gene often results in HF. PRDM16 deletion has been associated with left ventricular non-compaction and dilated cardiomyopathies. The importance of PRDM16 during heart development is being increasingly recognized, however the mechanisms by which PRDM16 loss leads to HF remain largely unknown. Numerous previous studies have found that altered cardiomyocyte mitochondrial function is a hallmark of many heart failure phenotypes. Therefore, a deeper understanding of the mechanistic relationship between cardiac mitochondrial bioenergetics and heart failure in the context of PRDM16 deletion may provide a target for therapies aimed at mitigating the degenerative changes associated with PRDM16 loss in 1p36 deletion syndrome. Exploring mechanistic changes of cardiac mitochondria is difficult to do experimentally so the aim of this project was to develop a computational model of cardiac mitochondria to use in combination with experimentation to elucidate the bioenergetic changes associated with the loss of PRDM16, in an e↵ort to further understand what might be leading to heart failure in PRDM16 null hearts. Thermodynamically constrained enzymatic reaction and transport flux equations were derived based on principles of Michaelis Menten kinetics. These flux equations were used to derive governing ordinary di↵erential equations, which were then used as the basis of the computational model. The parameters needed to fit these flux equations with the kinetics of specific enzymes and transporters were determined through a combination of literature values and experimental data. Once the model was parameterized, it was then used to quantify the mechanistic contribution of altered mitochondrial enzymes and transporters to cardiomyocyte bioenergetics in Prdm16-null hearts. Contents 1 Introduction 1 1.1 Computational Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Experimental Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Link Between Experimentation and Computational Model . . . . . . . . . . 4 2 Methods 2.1 2.2 7 Experimental . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Mitochondrial Isolation . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 Mitochondrial O2 Consumption Measurement . . . . . . . . . . . . . 8 Computational Model Development . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.1 Enzyme and Transport Reactions . . . . . . . . . . . . . . . . . . . . 9 2.2.2 Derivation of Enzyme Flux Equations . . . . . . . . . . . . . . . . . 10 2.2.2.1 Single Substrate Single Product Reactions . . . . . . . . . . 10 2.2.2.2 Two Substrate Two Product Reactions . . . . . . . . . . . . 11 2.2.2.3 Multi Substrate Multi Product Reactions . . . . . . . . . . 12 2.2.3 Derivation of Transport Flux Equations . . . . . . . . . . . . . . . . 14 2.2.4 Derivation of Ordinary Di↵erential Equations . . . . . . . . . . . . . 16 2.2.5 Estimation of Model Parameters . . . . . . . . . . . . . . . . . . . . . 18 3 Results 22 3.1 Experimental . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Computational . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4 Discussion 25 5 Next Steps and Future Work 26 6 Acknowledgements 28 iv 7 Appendix A 32 8 Appendix B 37 9 Appendix C 41 v 1 Introduction 1.1 Computational Modeling The use of computational models to supplement experimentation has become increasingly popular over recent years. The traditional experimental approach is to develop a hypothesis based on theoretical scientific background, run experiments, analyze data, adapt experiments to address findings and further questions, and then eventually end with a result. The idea behind using computational models is that it reduces the number of experiments required and streamlines the process. Computational models use principles of mathematics, physics, and biology to simulate a complex system. Experimental results are compared to simulated results and the di↵erences identified. Individual model parameters/reaction equations can be altered to match the experimental result. The specific parameters that are adjusted can provide great insight into what in the system might be changing, this then guides what experiments should be run to confirm the results from the simulation. This allows for a much more systematic and specific approach to understanding biological changes. 1.2 Experimental Background Heart failure (HF) is a complex clinical syndrome with symptoms and signs that result from structural or functional impairment of ventricular filling or ejection of blood [1]. 6.2 million adults are currently living with heart failure, and heart disease is the leading cause of death for both men and women in the United States and worldwide [2, 3]. HF can arise from diseases with environmental origins such as ischemic heart disease or metabolic diseases and/or from genetic/familial cardiomyopathies. Chromosome 1p36 deletion is the most common human chromosomal deletion syndrome and an estimated 30,000 - 200,000 people in the U.S with 1p36 deletion su↵er from HF [4]. PRDM16 is one of the genes frequently deleted in 1p36 syndrome, and loss of PRDM16 has been associated with pediatric left ventricular non-compaction and dilated cardiomyopathies 1 [5, 6, 7, 8, 9]. These findings have motivated further study of PRDM16 and its implications on cardiac development and HF. While PRDM16 is known to act as a transcriptional co-activator regulating hematopoiesis, palatogenesis, neuronal specification, and adipogenesis, the mechanisms by which PRDM16 loss leads to heart failure remain largely unknown. However numerous studies have found that altered cardiomyocyte mitochondrial function is a hallmark of many heart failure phenotypes. Mitochondria play a central role in energy production, homeostasis, and regulation of cell death. Mitochondria are the principal site of cellular energy production, as they are responsible for converting free energy from primary substrates into adenosine triphosphate (ATP) through oxidative phosphorylation (OXPHOS). The heart is the most energy consuming organ and mitochondria account for one third of the total volume of cardiomyocytes [10]. Therefore, a deeper understanding of the mechanistic relationship between cardiac mitochondrial bioenergetics and heart failure in the context of PRDM16 deletion may provide a target for therapies aimed at mitigating the degenerative cardiac changes associated with PRDM16 loss in 1p36 deletion syndrome. The production of ATP in mitochondria can be broken down into two main processes: tricarboxylic acid/citric acid cycle (TCA cycle) and electron transport chain (ETC) as seen in Figure 1. Glucose is metabolized into pyruvate through glycolysis in the cytoplasm and the pyruvate is shuttled into the mitochondria and processed by the TCA cycle to produce reducing equivalents used for OXPHOS. This is a complex process that is composed of several enzyme and transporter reactions. The TCA cycle reduces the coenzymes nicotinamide adenine dinucleotide (NAD+) and flavin adenine dinucleotide (FAD) into NADH and FADH2. These electron carriers are then passed to the ETC where they are reduced by the five complexes along with the transporters ubiquinone and cytochrome C. An electrochemical gradient potential ( p) is generated by the serial reduction of these electron carriers by the ETC. The reductive transfer of electrons through ETC complexes I–IV in the inner mitochondria membrane provides the energy to drive protons against the concentration gradient 2 across the inner mitochondrial membrane (out of the mitochondrial cytoplasm). This results in a net accumulation of H + outside the membrane, which then flows back into the mitochondria through the ATP synthase (Complex V), thus producing ATP and completing the ETC. The force driving the protons into the mitochondria ( p) is a combination of both the mitochondrial membrane potential ( a charge or electrical gradient) and the mitochon- drial pH gradient (H + concentration gradient). A simplified version of the Nernst equation can be used to represent p in terms of and H + . 3 Figure 1: Diagram illustrating the enzyme and transport equations occurring in the Mitochondria and the reactions that are modeled in the computational model. The model consists of three regions: the mitochondrial matrix, the intermembrane space, and the extracellular region or bu↵er region. The major biochemical species included in this model are: PYR: pyruvate, CoA: coenzyme-A, ACoA: acetyl-CoA, OXA: oxaloacetate, CIT: citrate, AKG: a-ketogluterate, SCoA: succinyl-CoA, SUC: succinate, FUM: fumarate, MAL: malate, GLU: glutamate, and ASP: aspartate, NAD and NADH: oxidized and reduced form of nicotinamide adenine dinucleotide, respectively, FAD and FADH2: oxidized and reduced form of Flavin adenine dinucleotide, respectively, ADP and ATP: adenosine triphosphate and adenosine diphosphate, respectively, UQ and UQH2: oxidized and reduced form of ubiquinone, respectively, CtyC3+ and CytC2+: oxidized and reduced form of cytochrome c, respectively 1.3 Link Between Experimentation and Computational Model A well-established method of assessing mitochondrial function is high resolution oxidative respirometry. Oxidative respirometry uses the fact that one of the last steps of cellular 4 respiration is the oxidation of cytochrome c in complex IV, which reduces oxygen to form water. While respirometry does provide a good insight into mitochondrial function, OXPHOS is a complex process that involves many simultaneous reactions, therefore identifying and isolating the specific mechanistic changes contributing to heart failure is difficult to do experimentally. The use of a computational model allows for greater insight into specific reactions by observing and identifying the changes in the parameters. Thus, the aim of this project was to develop a computational model of cardiac mitochondria and use it to simulate high resolution oxidative experiments so that the bioenergetic changes occurring with loss of PRDM16 could be identified with sufficient granularity. In order to have a useful and valid model of mitochondria bioenergetics it is important that it is both consistent with the available experimental data and that it follows the constraints of relevant physics/biophysics concepts. Modeling biochemical processes is complicated by the unique role enzymes play in the reaction kinetics of biochemical processes as compared to conventional chemical kinetics. To account for the di↵erence in kinetics, the foundation of most biochemical computational models are based on the well-established Michaelis Menten enzyme kinetic principles. These principles account for the dependency the reaction rate has on substrate concentration for reversible substrate product reactions. These principles can be used to derive single-substrate single-product reactions which then can be generalized into a multi-substrate multi-product reaction. In order to determine the rate of these reactions, they can be used in combination with the Haldane relation to derive thermodynamically closed flux equations. From the flux equations ordinary di↵erential equations (ODEs) can be derived which then can be solved to determine the concentration of all chemical species in a system. These ODEs are used as the basis of the computational model as the chemical species concentrations can be used to estimate the amount of respiration/oxygen consumption giving insight into the functionality of the mitochondria. This project was done in collaboration with the Boudina Lab at the University of Utah and the Dash Lab from the Medical College of Wisconsin. The project can be broken down 5 into three main steps: the development and parameterization of the computational model, the collection of experimental data, and the comparison of experimental data to simulated results to analyze the bioenergetic changes. The Boudina Lab which I have been a part of since 2021 provided the resources to collect the experimental respiration data and our collaborators from the Dash Lab, who specialize in computational modeling, provided the model framework to build our model of mitochondrial bioenergetics. For this project a computational model of isolated mitochondria respiration data was formulated based on 24 transporter/enzymatic equations that were derived based on the principles discussed above. From the flux equations 37 governing ODEs were derived. Then to parameterize the model to our groups of interest, intrinsic parameters were established based on published literature values and extrinsic (tissue dependent) parameters were determined experimentally with high-resolution oxidative respirometry (see methods section) and fitted using a least squares approximation. To examine the di↵erence between normal hearts and Prdm16 null hearts our lab generated cardiac-specific Prdm16 knockout mice (Prdm16 KO) with a Cre-Lox mouse model. Mitochondrial respiration data was collected from these Prdm16 KO mice as well as control mice (WT mice). The model was then parameterized for the control group and then reparameterized for the Prdm16 KO group. The specific parameters that changed were identified and analyzed allowing for a greater understanding of the mitochandrial bioenergetic changes that cause HF with respect to loss of PRDM16. 6 2 Methods 2.1 Experimental A Cre-Lox system was used to create a germline, cardiac-specific knockouts of Prdm16 on a C57B6 mouse background. Our lab has found that both male and female Prdm16 KO mice develop cardiac fibrosis, cardiomyocyte hypertrophy, and impaired cardiac function. Through echocardiography we have shown that female Prdm16 KO mice develop impaired systolic function between 8 and 20 weeks of age and die around 29 weeks. It is unclear which comes first, the mitochondrial bioenergetic changes or the structural remodeling. Thus, this project focuses on two groups, 14 week females and 20 week females, to help elucidate when and where the bioenergenic changes occur. Female controls and homozygous knockouts at 14 and 20 weeks were euthanized. Upon sacrifice, the heart were removed and the mitochondria was isolated by di↵erential centrifugation as described below with the tissue maintained at (4 °C) throughout the procedures. 2.1.1 Mitochondrial Isolation After the heart was removed it was briefly minced in ice-cold mitochondrial isolation bu↵er (MIB). The minced heart was then suspended in 1.5 mL MIB with 5 U/mL protease and homogenized until no intact tissue pieces were visible (⇠ 5 seconds). Next 1.5 mL cold MIB was added to the mixture and briefly mixed with a pipette, and transfered into two 2 mL round-bottom centrifuge tubes. The mixture was centrifuged at 8,000g for 10 minutes. The supernatant was discarded and 1.5 mL of clean MIB was added to each tube, and the pellet was re-suspended by vortexing and aspirating with a pipette. The mixture was centrifuged again at 8,000g for 10 minutes and the supernatant was again discarded and 1.5mL of clean cold MIB was added to each tube, and pellet re-suspended. The mixture was then centrifuged at 700g for 10 minutes and supernatant was transfered to a clean spin tube. It was then centrifuged at 8,000g for 10 minutes to pellet the mitochondria. The supernatant 7 was discarded and the pellet re-suspended in 100 uL of clean MIB. The mitochondrial protein content was determined using Bovine Serum Albumin (BSA) as the standard with the BioRad Quick Start Bradford Assay Kit. 2.1.2 Mitochondrial O2 Consumption Measurement Mitochondrial O2 consumption (respiration) was measured using an Oxygraph-2k (O2k) system (Oroboros Instruments, Innsbruck, Austria) with its DatLab 7 software used for data acquisition and analysis. Respiration was measured through single ADP addition experiments. 0.025 mg/ml of mitochondrial protein was suspended in respiration bu↵er (RB) within the O2k chamber (2 ml). Before each experiment, O2 concentration in RB was equilibrated for 2–3 minutes with air in the O2k chambers at 37 °C until a stable signal was obtained at an O2 concentration of approximately 160 µM. In all cases, experiments were initiated at t = 0 minutes with isolated mitochondria suspended in RB to reach state 1. The mitochondrial suspension was constantly stirred (750 rpm) inside the O2k chambers, and chemicals were added through the titanium injection port of the stoppers using Hamilton syringes. State 2 respiration was initiated with di↵erent substrate combinations at saturating concentrations, including pyruvate+malate (PM; 5 + 2.5 mM), glutamate+malate (GM; 5 + 2.5 mM), alpha-ketoglutarate+malate (AM; 5 + 2.5 mM), and succinate±rotenone (Suc + Rot; 10 mM + 0.5 µM). State 3 respiration was measured by adding 300 µM ADP. State 4 respiration was measured after the phosphorylation of added ADP to ATP. Oxygen consumption rates (OCR; JO2) during di↵erent states of respiration were calculated as the negative time derivative of O2 concentration measured in the closed O2k chambers and expressed as nmol of O2/min/mg protein. Data was recorded every 2 seconds and 40 data points were averaged to calculate the slope of the O2 concentration data using the DatLab 7 software. Cytochrome c was used to access the functional integrity of isolated mitochondria. 8 2.2 Computational Model Development The development of the model can be broken down into four main steps; the formulation of the enzyme and transport reactions, the derivation of the corresponding flux equations, the derivation of governing ODEs, and the parameter estimation and validation of the model. 2.2.1 Enzyme and Transport Reactions First, the enzymatic and transporter equations were formulated to represent the reactions occurring in the TCA cycle and ETC through principles of Michaelis Menten kinetics. As shown in Figure 1, the model consists of three regions (extracellular, intermembrane space, and mitochondrial matrix). In total 24 reactions were formulated, 14 enzyme reactions and 10 transporter reactions (Figure 2). Of the 14 enzyme reactions, 9 represent the reactions in the TCA cycle and the other 5 represent the 5 respiratory complexes in the ETC. The 10 transport equations account for the exchange of key metabolic species between the mitochondrial matrix and bu↵er region/cytoplasm. (a) Transport flux equations (b) Enzyme flux equations Figure 2: This figure contains two tables outlining the transport and enzyme reactions that the model uses as a basis.The subscript m indicates the mitochondrial matrix region and the subscript e indicates the extracellular region/bu↵er region. 9 2.2.2 2.2.2.1 Derivation of Enzyme Flux Equations Single Substrate Single Product Reactions After the enzymatic and transport reactions were formulated, they were used to derive flux equations to determine the rate of reactions. A single-substrate single-product reaction is illustrated below in Figure 3. Figure 3: E being enzyme, S being substrate, P being product and ESP being the hypothetical enzyme-substrate-product complex that binds both substrate and product. KS is the dissociation constant of ES complex and KP is the dissociation constant of EP. Reprinted “Integrated computational model of the bioenergetics of isolated lung mitochondria,” by X. Zhang, 2018, plos one, Volume 13, Supplemental Figure A1. To derive the flux equations it was assumed that all reactions followed generalized random-ordered rapid equilibrium kinetic mechanisms. This means all the enzymatic complex dissociation reactions are considered extremely rapid compared to the interconversion of ES and EP, which is the only rate limiting reaction step. It is also assumed that the binding of S does not a↵ect the binding of P. Under these two assumptions the following generalized flux equation can be derived: kp [ES] k p [EP ] kp [ES] k p [EP ] J = = Et Et [E] + [ES] + [EP ] + [ESP ] (1) Et is the total concentration of the enzyme. This equation can be further simplified by the fact that Vmaxf = kp Et and Vmaxr = k p Et are the maximum forward and reverse reaction velocities. The assumption that both the substrate and product do not bind to the enzyme 10 at the same time ([ESP] = 0) allows for the denominator to be simplified resulting in the following equation: J= [S] Vmaxf K S [P ] Vmaxr K P [S] [P ] ] 1+ K +K + K[S][P S P S KP = [S] Vmaxf K S [P ] Vmaxr K P [S] [P ] (1 + K )(1 + K ) S P = [S] Vmaxf K S [P ] Vmaxr K P [S] [P ] 1+ K S KP (2) For a reversible one-substrate one-product enzymatic reaction there are two intrinsic model parameters (enzymatic complex dissociation constants KS and KP ; also referred to as binding constants), and two extrinsic model parameters (forward and reverse reaction maximal velocities Vmaxf and Vmaxr ). The intrinsic parameters are assumed to be independent of the tissue that the enzymatic reaction is occurring in, while the extrinsic model parameters are dependent on this tissue. Specifically, the extrinsic model parameters are assumed to be dependent on the activity of the enzyme in a particular tissue, which is dependent on several cellular factors, such as enzyme expression, concentration, post-translational modifications, and physiological state of the tissue[11]. 2.2.2.2 Two Substrate Two Product Reactions A similar equation can be derived for a two-substrate two-product reaction (S1 + S2 E ! P1 + P2 ). J= 1 ][S2 ] (Vmaxf K[SS1 KS2 2] Vmaxr K[PP11][P ) KP 2 [S1 ] [S2 ] 2] (1 + K )(1 + K )(1 + K[PP11] )(1 + P[PS2 ) S1 S2 (3) It can be assumed that the S1 and P1 share the same binding site on the enzyme, and the same for S2 and P2 , as this is true for co-factor pairs like NADH and NAD+, and ATP and ADP. Therefore we can reduce the equation so that it excludes the contributions of [S1 ][P1 ] or [S2 ][P2 ] product terms: 11 J= 1 ][S2 ] (Vmaxf K[SS1 KS2 2] Vmaxr [PK1 ][P KP 2 ) P1 [S1 ] [S2 ] (1 + K + K[PP11] )(1 + K + K[PP22] ) S1 S2 (4) Like the one-substrate one-product there are again both intrinsic and extrinsic model parameter but in this case there are four intrinsic and two extrinsic parameters instead of two of each. 2.2.2.3 Multi Substrate Multi Product Reactions There are also a lot of reactions that are comprised of more than two substrates or products. So in a similar fashion we can derive a multi-substrate multi-product reactions, represented by NS X ↵ i Si i=1 ! NP X j Pj (5) i=1 With Si is ith substrate, Pj is jth product, NS and NP are the number of substrates and products in the reaction, respectively, and ↵i and j are the corresponding stoichiometric coefficients. The method for deriving flux equations for these multi-substrate multi-product reactions is very similar to the process for single-substrate single-product reactions. The multi-substrate multi-product reaction is below: QNS Vmaxf Qi=1 NS J= Q [S]↵i ↵i i=1 KSi QNP Vmaxr Qj=1 NP [P ] j j j=1 KPj NS [S]↵i QNP [P ] j ) i=1 (1 + K ↵i ) j=1 (1 + KP j S (6) To maintain the thermodynamic consistency of the enzymatic reactions the Haldane relationship is used. The general Haldane equation is below: 0 = Keq Vmaxf · KP Vmaxr · Ks 12 (7) This relationship constrains the forward and reverse reaction velocities to maintain the thermodynamic consistency. We can apply the Haldane relationship to a multi-substrate multiproduct reaction and use the fact that at equilibrium the net reaction flux is zero (J = 0), to get the following thermodynamic constraint: Vmaxf Q NP j j=1 KPj 0 = Keq Q NS ↵ i Vmaxr i=1 KSi = Q NP j j=1 [Pj ]eq Q NS ↵i i=1 [Si ]eq (8) 0 Keq is the equilibrium constant for a given reaction, which is the reaction equilibrium state 0 quotient at specified thermodynamic conditions. This constant is for pH=7 (Keq0 (pH = 7)) and can be calculated as 0 Keq0 = e Where rG 00 rG 00 /RT (9) , R and T are the standard Gibbs free energy of the reaction at pH = 7, gas constant, and temperature, respectively. The standard Gibbs free energy of reactions in the mitochondrial bioenergetics model was obtained from literature values [12, 13, 14]. The changing temperature can be accounted for using the Q10 temperature coefficient, T2 T1 R2 = Q10 10 R1 (10) where T is the temperature, R2 /R1 is the correction coefficient (i.e. the ratio of maximal reaction and transport velocities at temperature T2 and T1 ), and Q10 is the temperature coefficient which was set to 2. Using the Haldane constraint, Vmaxr can be calculated from 0 Vmaxf and Keq0 to reduce the number of unknown model parameters. Thus, the reaction flux expression for multi-substrate and multi-product enzymatic reaction becomes: 13 QNS Vmaxf ↵ QNS ↵i ( i=1 [Si ] i K i=1 S J= Q i QNP j j=1 [Pj ] 0 Keq ) NS [Si ]↵i QNP [P ] j (1 + ) ↵ ) i=1 j=1 (1 + KS i KP j i (11) j Again there are both intrinsic and extrinsic model parameter unknowns. And again there are just two extrinsic parameters, the forward and reverse Vmax , but there are i and j binding constants. 2.2.3 Derivation of Transport Flux Equations Transport equations are derived in a similar fashion as the enzyme equations. It can again be assume that all transporters have a random-ordered rapid-equilibrium binding mechanism, i.e the transporter can bind to the two substrates in an arbitrary order. In general, there are two di↵erent types of transporters in this system: a co-transporter (symporter) that transports two di↵erent substances from one side to the other side of the inner mitochondrial membrane, and an antiporter (exchanger) that transports two di↵erent substances across the inner mitochondrial membrane in the opposite directions. The general scheme for a co-transporter is shown in the diagram below. Figure 4: Ai and Bi are two distinct substrates at side i of the membrane (i =1, 2). The free transporter (E) binds to the two substrates in an arbitrary order before undergoing conformational changes. Reprinted “Integrated computational model of the bioenergetics of isolated lung mitochondria,” by X. Zhang, 2018, plos one, Volume 13, Supplemental Figure A3. 14 All association and dissociation processes are assumed to be reversible and rapid-equilibrium reactions. The rate limiting steps are the conformational changes in the substrate-bound transporters resulting in the transport of the substrates across the membrane. The general flux equations is below: Jco = 1] Tmaxf K[AA11 ][B KB1 2] Tmaxr K[AA22 ][B KB2 [A1 ] [A1 ] [B1 ] [B1 ] 1] 2] +K +K +K + K[AA11 ][B + K[AA22 ][B 1+ K KB1 KB2 A1 A2 B1 B2 (12) where Tmaxf and Tmaxr are the maximum forward and reverse transport rates, respectively. The general scheme for an anti-porter is shown below. Figure 5: Ai and Bi are two distinct substrates at side i of the membrane (i =1, 2). The free transporter (E) binds to the two substrates in an arbitrary order before undergoing conformational changes.Reprinted “Integrated computational model of the bioenergetics of isolated lung mitochondria,” by X. Zhang, 2018, plos one, Volume 13, Supplemental Figure A4. The anti-porters bind to the substrates from both sides of the membrane. The general flux equation for an anti-porter is similarly given by Janti = 2] Tmaxf K[AA11 ][B KB2 1] Tmaxr K[AA22 ][B KB1 [A1 ] [A1 ] [B1 ] [B1 ] 2] 1] 1+ K +K +K +K + K[AA11 ][B + K[AA22 ][B KB2 KB1 A1 A2 B1 B2 (13) For all the metabolite transporters, the internal and external binding constants are assumed 15 to be equal (KA1 = KA2 and KB1 = KB2 ). Therefore: Tmaxf 2 ][B2 ] ([A1 ][B1 ] [AK KA KB eq Jco = [A1 ] [A1 ] [B1 ] [B1 ] [A1 ][B1 ] 2 ][B2 ] 1 + KA + KA + KB + KB + KA KB + [A KA KB Tmaxf 2 ][B1 ] ([A1 ][B2 ] [AK KA KB eq Janti = [A1 ] [A1 ] [B1 ] [B1 ] [A1 ][B2 ] 2 ][B1 ] 1 + KA + KA + KB + KB + KA KB + [A KA KB (14) (15) For non-electrogenic transporters, Tmaxf = Tmaxr and Keq = 1. For electrogenic transporters, such as the glutamate-aspartate exchanger (GAE) and the adenine nucleotide translocase (ANT), equilibrium constants are dependent on membrane potential, determined by the Nernst equation. All flux equations for transporters and enzymes are listed in Appendix A[15, 11]. 2.2.4 Derivation of Ordinary Di↵erential Equations Ordinary di↵erential equations were derived from the flux equations to describe the dynamic changes in concentration of the 37 state variables, 34 biochemical species and three ODEs to account for ion concentration and the membrane potential relationship. The ODEs were derived based on principles of mass balance. The change in the concentration of a given biochemical species within the mitochondrial matrix region and extra-mitochondria (bu↵er) region are given by: Vm dCm,j X = ↵m,j Jm,j dt Ve X dCe,j X = Jm e,j dt Jm e,j (16) (17) where Jmj is the jth reaction flux in the mitochondrial matrix and Jm e,j is the jth transport flux between the mitochondrial matrix region and the extra-mitochondrial bu↵er region, Vm 16 and Ve are the volumes of the mitochondrial matrix region and the extra-mitochondrial region, respectively. For the extra-mitochondrial region which does not include chemical reactions, the right hand side of the equation contains just transport fluxes. The following ODE was used to account for the role of membrane potential in the mitochondria: d m dt = 1 Cimm (4JCI + 2JCIII + 4JCIV 3JCV JHLEAK JAN T ) (18) where Cimm is the capacitance of the inner mitochondrial membrane; JCI , JCIII , JCIV and JCV are the reaction fluxes of complex I, complex III, complex IV and complex V, respectively;JHleak is the proton leak between the IMS and the mitochondrial matrix regions; and JAN T is the transport flux characterizing ATP/ADP exchange via the ANT. Accounting for membrane potential is essential, as it plays such a large role in ATP production. There are several existing models of mitochondria bioenergetics that fail to properly account for membrane potential leading to inconsistencies when being compared with both the available experimental data and the constraints of relevant physics/biophysics concepts. For example, the Korzeniewski and Zoladz model is widely used and has been validated experimentally with several studies[16, 17, 18]. However, it relies on the assumption that there is a linear relationship between the di↵erence in pH across the mitochondrial inter membrane and the magnitude of the inner membrane potential. A thorough study of isolated mitochondria showed that this linear relationship is not obeyed, and the model breaks down because it cannot account for the fact that the matrix can become more acidic as the magnitude of the membrane potential increases[19]. Another widely used model is the Magnus–Keizer model. The Magnus–Keizer model treats pH as constant and does not account for pH bu↵ering by potassium ion exchange, so again the model does not hold up in regard to pH and membrane potential. These issues in membrane potential are accounted for by incorporated phosphate dependent control of oxidative phosphorylation[20]. 17 As the model framework was adapted from our collaborators at The Medical College of Wisconsin, they validated this method of accounting for membrane potential by running membrane potential assays and comparing it to simulations with these equations. The detailed mass balance equations for each chemical species are listed in Appendix B [15, 11]. MATLAB function ’ode15s’ was used to solve the system of governing ordinary di↵erential equations. 2.2.5 Estimation of Model Parameters As previously mentioned there are both intrinsic (non-tissue specific) and extrinsic (tissue specific) model parameters. The intrinsic parameters are the binding constants KS and KP associated with a given substrate/product. They were estimated based on published literature values [14, 11]. The binding constants and other model constants are listed in Appendix A along with the flux equations. The extrinsic model parameters (Vmax ) were estimated using a combination of experimental OCR data and least squares fitting. Initial Vmax approximations were adapted from the literature and values our collaborators at the Medical College of Wisconsin collected for di↵erent rodent heart tissue. These initial approximations were then refined by solving the following system of linear equations (of the form Ax = b) with a least squares approximation: 2 6 6 6 6 6 6 6 6 6 6 6 4 32 OCRt1 Vmax1 OCRt1 Vmax2 OCRt1 Vmax24 7 6 ... : : ... : : ... : : ... : OCRtn Vmax1 OCRtn Vmax1 ... OCRtn Vmax24 : : 76 76 76 76 76 76 76 76 76 76 54 3 2 3 Vmax1 7 6 OCRt1 7 7 6 7 7 6 7 : : 7 6 7 7 6 7 7 6 7 = 7 6 7 : : 7 6 7 7 6 7 7 6 7 : : 7 6 7 5 4 5 Vmax24 OCRtn The A matrix is a derivative approximation, with OCR representing the di↵erence be- tween the simulated data with the initial approximated Vmax values OCR0 and the simulated data after each Vmax value was adjusted one at a time by a 10% reduction, which is the 18 Vmax value. In matrix B OCR represents the di↵erence in our experimental data that we are trying to fit (target data) and the simulated values produced by our initially approximated Vmax values. The least squares approximation was used to solve for the x matrix ( Vmax ) producing the change in the Vmax values needed to improve the simulation’s fit. To further improve upon the estimated values, the least squares approximation was repeated but with the initial approximated Vmax and simulated values being updated each time by the previous least squares result. Meaning that for both matrix A and B OCR0 becomes OCR1 and so on for each least squares iteration. This process of repeatedly replacing the initial approximation was used until the least squares approximation was no longer improving the fit. This was assessed by checking the mean squared error for each iteration, once Vmax values were the error reached its minimum the iterations were stopped and the final recorded. To assesses the contribution of each of the extrinsic model parameters to the overall model solution a sensitivity analysis was performed on the 24 Vmax parameters. The sensitivity matrix was constructed as follows: 2 OCRt1 Vmax1 6 OCRt1 · Vmax1 OCRt1 · VVmax2 OCRt1 max2 6 6 6 6 6 6 6 6 6 6 4 ... 3 OCRt1 · VVmax24 OCRt1 max24 7 : : : ... : : ... : : ... : OCRtn · VVmax1 OCRtn max1 OCRtn · VVmax2 OCRtn max2 ... OCRtn · VVmax24 OCRtn max24 : 7 7 7 7 7 7 7 7 7 7 5 OCR is the di↵erence between OCR produced by the simulation with a set of initial Vmax values and the simulated data produced when each Vmax is changed by a given factor, in this case a 10% decrease. OCRt1 . . . OCRtn represents the OCR data produced by the simulation at each time point. OCR OCR is then multiplied by change in Vmax divided by the original Vmax . This produces a n X 24 matrix for n time points. Small values indicate that a parameter is insensitive and large values indicate sensitivity. This analysis indicated that the 19 Vmax parameters with the greatest sensitivity were the ones associated with the pyruvate dehydrogenase reaction, the citrate synthase reaction, the glutamate oxaloacetate reaction, the succinate dicarboxylate carrier, the malate dicarboxylate carrier and the inorganic phosphate carrier. This information about parameter sensitivity was taken into account for improving the least squares solution. Relatively small changes in certain Vmax values produced no notable change as they had low sensitivity. To confirm that the changes of the insensitive parameters had little/no e↵ect on the simulation out come, I wrote a MATLAB loop to increase and decrease the Vmax value by several orders of magnitudes. If there was indeed no e↵ect on the simulation output the that Vmax was dropped. Figure 6 is the plotted simulations from this loop. Figure 6a is the plot of a parameter that was indicated to be sensitive and Figure 6b is a plot of a parameter that was indicated to be less sensitive. (a) Inorganic phosphate carrier Vmax (b) Isocitrate dehydrogenases Vmax Figure 6: This is a plot of 12 iterations of one Vmax being decreased/increased by 10% each time. On the left in Figure (a) the highly sensitive inorganic phosphate carrier Vmax is being adjusted and the simulation plotted for each substrate. On the right in Figure (b) the insensitive Isocitrate dehydrogenases Vmax is plotted. This loop confirms the information from the sensitivity matrix as these adjustments cause a large change in the simulation for the inorganic phosphate carrier Vmax but very little change in the simulation for the Isocitrate dehydrogenases Vmax . There is a slight change in the pyruvate peak but no other changes for the Isocitrate dehydrogenases Vmax 20 The simulation with the parameterized Vmax values obtained from the least squares solution and the experimental data for the control group are plotted below in Figure 7: (a) (b) (c) (d) Figure 7: Plots of experimental data from the control mouse (WT) and the simulation parameterized to this group for each substrate. Experimental data is in blue and simulation is in orange. Plot (a) is pyruvate and malate, plot (b) is glutatmate and malate, plot (c) is alpha ketoglutarate and malate and plot (d) is succinate and rotenone. 21 3 Results 3.1 Experimental Respiratory states and rates were measured for 14 and 20 week female Prdm16 wildtype mice and Prdm16 KO mice (n = 1). In the presence of various respiratory substrate combinations. Figure 8 shows oxygen consumption rates for each group for each respiratory state (state 2-4). The 14 week females did not display a notable change in any of the respiratory states for any of the substrates. However the 20 week Prdm16 KO mice exhibited a 25% decrease in state 3 respiration for all substrates and no significant change in state 2 or state 4 respiration. (a) (b) (c) (d) Figure 8: Rate of oxygen consumption in each respiratory state for Prdm16 14 and 20 week control mice and Prdm16 KO females measured in nmol/min/mg of protein. 22 3.2 Computational In order to interrogate the bioenergetic changes that might be causing the decrease in state 3 respiration for the 20 week group, the di↵erences in Vmax values between Prdm16 KO mice and control mice were assessed. In order to do this the same fitting process was used for the Prdm16 KO data as was used for the initial parameterization of the model. There was no di↵erence in the Vmax values for the 14 week group Prdm16 KO and control group however there were di↵erences in the Vmax values between the control and Prdm16 KO mice in the 20 week group. The simulation parameterized for the Prdm16 KO group and experimental data form the Prdm16 KO can be seen in Figure 9. The model parameters for a control mouse and parameters for a Prdm16 KO mouse can be seen in Figure 10. (a) (b) (c) (d) Figure 9: Plots of experimental data and simulated data for the Prdm16 KO mice for each substrate (pyruvate and malate, glutatmate and malate, alpha ketoglutarate and malate, and succinate and rotenone). 23 Figure 10: Table of Vmax values by reaction, column 2 are the Vmax values for the control mice and column 3 is the values for the Prdm16 KO mice. The major changes are indicated by highlighting in red. The Vmax values for complex II (CII) and the inorganic phosphate carrier (PIC) were both decreased. It can be seen in Figure 10 there are two parameters that changed significantly. The Vmax for complex II (CII) and the inorganic phosphate carrier (PIC) both decreased. CII had a 28% decrease and PIC had a 51% decrease. A decrease in Vmax indicates that the reaction is occurring less rapidly. Indicating that there is some impairment of complex II and the transport of inorganic phosphate, that is causing the mitochondrial function in Prdm16 null hearts to be less e↵ective at metabolizing substrates. The change being in complex II and the inorganic phosphate carrier indicates that the bioenergetic changes have to do with the ETC and not an direct impairment of the TCA cycle. 24 4 Discussion The purpose of this project was to develop a thermodynamically closed computational model to elucidate the mitochondrial bioenergetic changes occurring with loss of PRDM16. The interest in PRDM16 in cardiac function has become more prominent over recent years but a lot about it is still unknown. Before this paper no one had used a combined computational and experimental approach to assess mitochondrial function with respect to PRDM16. The value of developing and using a computational model as a tool to assesses mitochondrial bioenergeic changes, is that it allows for a significant reduction in the number of experiments that need to be run. The changes in model parameters, gives focus to future experiments. It also allows for the isolation and manipulation of particular mechanisms rather easily to explore potential options for further experimentation. I was able to successfully develop and parameterize a computational model for cardiac mitochondria in female mice. The model used principles of Michaelis Menten kinetcs, the Haldane relation, the Nernst equation, and Gibbs free energy to derive governing ODEs. Parameterizing the model for our specific experimental groups required a great deal of trial and error and was a beneficial learning experience. Initially, I attempted to use a least squares solution with the simulation parameters provided by our collaborators. However as their parameters were derived based on rat tissue it turned out they were not close enough to the true values to converge properly. I then performed a sensitivity analysis on the 24 parameters to identify which parameters would influence the model outcomes the most. I then wrote a MATLAB loop to identify the specific behaviors that result from adjusting the parameters by orders of magnitude. From this I developed a new initial guess that was closer to my target values and then used the same set up for the least squares approximation, this resulted in a much better fit for the model parameters. Then to further improve the model fit this initial guess was replaced continuously and the least squares approximation repeated. This method seemed relatively e↵ective and is what I used to re-parameterize from the control group to the Prdm16 KO group. 25 The initial results from the model indicate that Prdm16 KO leads to mitochondrial bioenergetic changes at 20 weeks. The experimental data showed that respiration was down in the Prdm16 null hearts. However by analyzing just the experimental data it was unclear what would be causing the decrease in respiration. By pairing the data with the computational model simulation we were able to see that the loss of PRDM16 seems to lead to the impairment of the ETC. This is indicated by the decrease in the complex II and inorganic phosphate carrier Vmax values. Both of these reactions are part of the ETC. This directs further experimentation to focusing on these two particular aspects of the ETC. Specifically, the next experimental steps would be to run protein blots for the di↵erent complexes of the ETC to determine if the issue is one of frequency rather than dysfunction of the complexes. If the blots indicate it is not an issue of frequency the dysfunction can be assessed with target enzyme kinetic experiments and inhibition of the di↵erent complexes one at a time. Inhibitors such as Antimycin A (inhibits complex III), Oligomycin and Citreoviridin (inhibits complex V), Mersalyl (inhibits inorganic phospahte transport), Phenylsuccinate (inhibits succinate transporter e↵ecting complex 2), and Sodium azide (inhibits complex IV) can be used. This project was successful in developing the model to help elucidate potential changes the next steps are just to validate these changes with further experimentation. 5 Next Steps and Future Work Though the model is useful, there are a variety of improvements that can be made. It would be beneficial to run more experimental data to compare to the model simulation for the Prdm16 KO group to confirm that the changes are consistent and occurring for all mice of these genotypes and ages. To expand the utility of the model in general it would also be useful to parameterize this model for a more standard control group like C57 BL/J6 male mice of di↵erent age ranges as this is the standard control group for most mouse model experiments. That way 26 this model can be used to assesses a variety of other cardiac phenotypes with mouse models with C57 backgrounds. There is also room to expand on the model itself as at the moment it does not take into account -oxidation. -oxidation has been identified as extremely important in the heart with it accounting for 70% of ATP production. Adapting this model to incorporate -oxidation would be useful in general but to our lab’s research specifically as we have identified that oxidation might play an important role in changes with PRDM16 loss. We have collected experimental respirometry data with palmitoyl-l-carnitine (long chain fatty acid) as a substrate. Its state 3 respiration also appears to be down indicating that there is a change in oxidation that would be useful to explore. I am excited to see how we can further adapt and utilize this model now that we have a solid framework. 27 6 Acknowledgements I would like to thank the Dash Lab at Wisconsin Medical College’s for collaborating on this project and helping provide the framework for the computational model. As well as the Boudina lab here at the University of Utah and Postdoctoral fellow Benjamin Werbner for assisting in collecting experimental data. 28 References [1] Biykem Bozkurt, Ray E Hershberger, Javed Butler, Kathleen L Grady, Paul A Heidenreich, Maria Lizza Isler, James K Kirklin, and William S Weintraub. 2021 acc/aha key data elements and definitions for heart failure: a report of the american college of cardiology/american heart association task force on clinical data standards (writing committee to develop clinical data standards for heart failure). Journal of the American College of Cardiology, 77(16):2053–2150, 2021. [2] V SS et al. on behalf of the american heart association council on epidemiology and prevention statistics committee and stroke statistics subcommittee. Heart disease and stroke statistics-2020 update: a report from the American Heart Association. Circulation, 141:e1–e458, 2020. [3] SS Virani, A Alonso, EJ Benjamin, et al. Coronary heart illness and stroke statistics-2020 replace: a report from the american coronary heart affiliation. Circulation, 141:e139– e596, 2020. [4] Chromosome 1p36 deletion syndrome - about the disease. https://rarediseases.info.nih. gov/diseases/6082/chromosome-1p36-deletion-syndrome/. [5] Anne-Karin Arndt, Sebastian Schafer, Jorg-Detlef Drenckhahn, M Khaled Sabeh, Eva R Plovie, Almuth Caliebe, Eva Klopocki, Gabriel Musso, Andreas A Werdich, Hermann Kalwa, et al. Fine mapping of the 1p36 deletion syndrome identifies mutation of prdm16 as a cause of cardiomyopathy. The American Journal of Human Genetics, 93(1):67–77, 2013. [6] Tongbin Wu, Zhengyu Liang, Zengming Zhang, Canzhao Liu, Lunfeng Zhang, Yusu Gu, Kirk L Peterson, Sylvia M Evans, Xiang-Dong Fu, and Ju Chen. Prdm16 is a compact myocardium-enriched transcription factor required to maintain compact myocardial cardiomyocyte identity in left ventricle. Circulation, 145(8):586–602, 2022. 29 [7] Josef Finsterer, Claudia Stoellberger, and Je↵rey A Towbin. Left ventricular noncompaction cardiomyopathy: cardiac, neuromuscular, and genetic factors. Nature Reviews Cardiology, 14(4):224–237, 2017. [8] Je↵rey A Towbin, Angela Lorts, and John Lynn Je↵eries. Left ventricular non- compaction cardiomyopathy. The Lancet, 386(9995):813–825, 2015. [9] Pamela A Long, Jared M Evans, and Timothy M Olson. Diagnostic yield of whole exome sequencing in pediatric dilated cardiomyopathy. Journal of Cardiovascular Development and Disease, 4(3):11, 2017. [10] Namrata Tomar, Xiao Zhang, Sunil M Kandel, Shima Sadri, Chun Yang, Mingyu Liang, Said H Audi, Allen W Cowley Jr, and Ranjan K Dash. Substrate-dependent di↵erential regulation of mitochondrial bioenergetics in the heart and kidney cortex and outer medulla. Biochimica et Biophysica Acta (BBA)-Bioenergetics, 1863(2):148518, 2022. [11] Xiao Zhang, Ranjan K Dash, Elizabeth R Jacobs, Amadou KS Camara, Anne V Clough, and Said H Audi. Integrated computational model of the bioenergetics of isolated lung mitochondria. PLoS One, 13(6):e0197921, 2018. [12] Yanjun Li, Nicola Lai, John P Kirwan, and Gerald M Saidel. Computational model of cellular metabolic dynamics in skeletal muscle fibers during moderate intensity exercise. Cellular and molecular bioengineering, 5:92–112, 2012. [13] Ranjan K Dash, Yanjun Li, Jaeyeon Kim, Daniel A Beard, Gerald M Saidel, and Marco E Cabrera. Metabolic dynamics in skeletal muscle during acute reduction in blood flow and oxygen supply to mitochondria: in-silico studies using a multi-scale, top-down integrated model. PloS one, 3(9):e3168, 2008. [14] Fan Wu, Feng Yang, Kalyan C Vinnakota, and Daniel A Beard. Computer modeling of mitochondrial tricarboxylic acid cycle, oxidative phosphorylation, metabolite transport, and electrophysiology. Journal of Biological Chemistry, 282(34):24525–24537, 2007. 30 [15] Shima Sardina and Ranjan K Dash. Untitled manuscript. unpublished/personal communication. [16] Valdur A Saks, O Kongas, Marko Vendelin, and Laurence Kay. Role of the creatine/phosphocreatine system in the regulation of mitochondrial respiration. Acta Physiologica Scandinavica, 168(4):635–641, 2000. [17] Marko Vendelin, Olav Kongas, and Valdur Saks. Regulation of mitochondrial respiration in heart cells analyzed by reaction-di↵usion model of energy transfer. American Journal of Physiology-Cell Physiology, 278(4):C747–C764, 2000. [18] Satoshi Matsuoka, Nobuaki Sarai, Hikari Jo, and Akinori Noma. Simulation of atp metabolism in cardiac excitation–contraction coupling. Progress in Biophysics and Molecular Biology, 85(2-3):279–299, 2004. [19] Salil Bose, Stephanie French, Frank J Evans, Fredric Joubert, and Robert S Balaban. Metabolic network control of oxidative phosphorylation: multiple roles of inorganic phosphate. Journal of Biological Chemistry, 278(40):39155–39165, 2003. [20] Daniel A Beard. A biophysical model of the mitochondrial respiratory system and oxidative phosphorylation. PLoS computational biology, 1(4):e36, 2005. 31 7 Appendix A Below are the 24 flux equations for each of the enzymatic/transport reactions. The Gibbs free energy and binding constants for each reaction are listed directly after the equation. The free energy barrier for all the complexes is 0.55 J. The free energy barrier for Adenine nucleotide translocase (ANT) is 0.6 J. Vmaxf (ABC DEF 0 ) KA KB KC Keq JP DH = (1 + KAA )(1 + KBB + KDD )(1 + KCC + KFF )(1 + KEE ) (1) A = [P Y Rm ], B = [CoAm ], C = [N ADm ], D = [ACoAm ], E = [CO2,m ], F = [N ADHm ] The binding constants in molar are KA = 5.08⇤10 3 , KB = 1.49⇤10 5 , KC = 3.5⇤10 5 , KD = 1.49 ⇤ 10 5 , KE = 3.5 ⇤ 10 5 , KF = 1 ⇤ 10 5 and the Gibbs free energy at pH = 7 for this reaction is -39.64 kJ/mol. Vmaxf CD (AB K 0 ) KA KB eq JCIT S = (1 + KAA + KCC )(1 + KBB )(1 + KDD ) (2) A = [ACOAm ], B = [OXAm ], C = [CoAm ], D = [CITm ] The binding constants in molar are KA = 3.9⇤10 6 , KB = 4.53⇤10 6 , KC = 57.9⇤10 6 , KD = 4.3 ⇤ 10 3 and the Gibbs free energy at pH = 7 for this reaction is -36.61 kJ/mol. Vmaxf (AB CDE 0 ) KA KB Keq JCIT DH = (1 + KAA )(1 + KCC )(1 + KBB + KDD )(1 + KEE ) (3) A = [CITm ], B = [N ADm ], C = [AKGm ], D = [CITm ], E = [CO2,m ] The binding constants in molar are KA = 1.3 ⇤ 10 3 , KB = 500 ⇤ 10 6 , KC = 3.5 ⇤ 10 6 , KD = 4.7 ⇤ 10 6 , KE = 1 ⇤ 10 5 and the Gibbs free energy at pH = 7 for this reaction is 3.29 kJ/mol. Vmaxf (ABC DEF 0 ) KA KB KC Keq JAKGDH = (1 + KCC + KEE )(1 + KBB + KDD )(1 + KAA ) 32 (4) A = [AKGm ], B = [CoAm ], C = [N ADm ], D = [SCoAm ], E = [N ADHm ], F = [CO2 ] The binding constants in molar are KA = 85.87 ⇤ 10 6 , KB = 1.634 ⇤ 10 6 , KC = 46.6 ⇤ 10 6 , KD = 1.63 ⇤ 10 6 , KE = 7.14 ⇤ 10 6 , KF = 1 ⇤ 10 5 and the Gibbs free energy at pH = 7 for this reaction is -37.08 kJ/mol. Vmaxf (ABC DEF 0 ) KA KB KC Keq JSCAS = (1 + KAA + KDD + KFF + KDF )(1 + KBB + KCC + KEE + KBC ) D KF B KC (5) A = [SCoAm ], B = [GDPm ], C = [P im ], D = [SU Cm ], E = [GT Pm ], F = [CoAm ] The binding constants in molar are KA = 1.2538 ⇤ 10 5 , KB = 5.6977 ⇤ 10 6 , KC = 2.1 ⇤ 10 3 , KD = 4.51 ⇤ 10 4 , KE = 2.3143 ⇤ 10 5 , KF = 1.7735 ⇤ 10 5 and the Gibbs free energy at pH = 7 for this reaction is 1.26 kJ/mol. Vmaxf CD (AB K 0 ) KA KB eq JN DK = (1 + KAA + KCC )(1 + KBB + KDD ) (6) A = [GT Pm ], B = [ADPm ], C = [GDPm ], D = [AT Pm ] The binding constants in molar are KA = 111⇤10 6 , KB = 100⇤10 6 , KC = 260⇤10 6 , KD = 278 ⇤ 10 6 and the Gibbs free energy at pH = 7 for this reaction is -.56 kJ/mol. Vmaxf (A KB0 ) KA KB eq JF U M = A (1 + KA )(1 + KBB ) (7) A = [F U Mm ], B = [M ALm ] The binding constants in molar are KA = 1800 ⇤ 10 6 , KB = 1800 ⇤ 10 6 and the Gibbs free energy at pH = 7 for this reaction is -3.6 kJ/mol. Vmaxf CD (AB K 0 ) KA KB eq JM DH = (1 + KBB + KCC )(1 + KAA + KDD ) (8) A = [N ADm ], B = [M ALm ], C = [OXAm ], D = [N ADHm ] The binding constants in molar are KA = 1.55⇤10 4 , KB = 1.1⇤10 3 , KC = 3.38⇤10 6 , KD = 33 3.61 ⇤ 10 5 and the Gibbs free energy at pH = 7 for this reaction is 28.83 kJ/mol. Vmaxf CD (AB K 0 ) KA KB eq JGOT = (1 + KAA )(1 + KBB )(1 + KCC )(1 + KDD ) (9) A = [ASPm ], B = [AKGm ], C = [GLUm ], D = [OXAm ] The binding constants in molar are KA = 3.9 ⇤ 10 3 , KB = 430 ⇤ 10 6 , KC = 8.9 ⇤ 10 3 , KD = 88 ⇤ 10 6 the Gibbs free energy at pH = 7 for this reaction is -1.31 kJ/mol. JCI = Vmaxf (exp( KA KB CI (4 GH + RT 00 r GCI ) ( CI 1)(4 GH + RT A C B D (1 + KA + KC )(1 + KB + KD ) )AB exp( 00 r GCI ) )CD) (10) A = [N ADHm ], B = [U Qm ], C = [N ADm ], D = [U QH2m ] The binding constants in molar are KA = 1.5⇤10 6 , KB = 58.1⇤10 6 , KC = 428⇤10 6 , KD = 520 ⇤ 10 6 the Gibbs free energy at pH = 7 for this reaction is -69.37 kJ/mol. Vmaxf CD (AB K 0 ) KA KB eq JCII = (1 + KBB + KCC )(1 + KAA + KDD ) (11) A = [F ADH2,m ], B = [U Qm ], C = [F ADm ], D = [U QHm ] The binding constants in molar are KA = 1800 ⇤ 10 6 , KB = 140 ⇤ 10 6 , KC = 1800 ⇤ 10 6 , KD = 2.45 ⇤ 10 6 the Gibbs free energy at pH = 7 for this reaction is -2.41 kJ/mol. JCIII = For Y = ( 4RTGH + 2FRT Vmaxf 2 2 (exp( CIII (Y )AB KA KB exp( CIII (1 + KAA + KCC )(1 + 00 r GCIII RT B2 2 KB + D2 2 KD 1)(Y )) ) (12) ) and A = [U QH2,m ], B = [CytCo], C = [U Qm ], D = [CytCr2m ] The binding constants in molar are KA = 4.66⇤10 6 , KB = 3.76⇤10 6 , KC = 4.08⇤10 6 , KD = 4.91 ⇤ 10 6 the Gibbs free energy at pH = 7 for this reaction is -32.53 kJ/mol. JCIV = Vmaxf 2 .5 2 K .5 (exp( IV (Y )A B KA B 2 exp( CIV 2 2 .5 B (1 + KA2 + KCC )(1 + K .5 ) A 34 B 1)(Y )C 2 ) (13) For Y = 2 GH RT 00 r GCIV 2F RT RT ) and A = [CytCr], B = [O2 ], C = [CytCo] The binding constants in molar are KA = 680 ⇤ 10 6 , KB = 5.4 ⇤ 10 6 , KC = 79.2 ⇤ 10 6 and the Gibbs free energy at pH = 7 for this reaction is -122.94 kJ/mol. JCV = 3 ( GH Vmaxf (exp( V KA KB RT 00 r GCV ) )AB exp( 3( CV 1)( GH RT 00 r GCV ) )C (1 + KAA + KCC )(1 + KBB ) (14) A = [ADPm ], B = [P im ], C = [AT Pm ] The binding constants in molar are KA = 50 ⇤ 10 6 , KB = 3 ⇤ 10 3 , KC = 50 ⇤ 10 6 and the Gibbs free energy at pH = 7 for this reaction is 36.03 kJ/mol. e] Tmax ( [PKYP YRRe ][H KH [P Y Rm ][Hm ] ) KP Y R KH JP Y RH = e] m] e] m] 1 + [PKYP YRRe ] + [PKYP RY Rm ] + [H + [H + [PKYP YRRe ][H + [PKYPRYmR][H KH KH KH KH (15) The binding constants in molar are KP Y R = 0.24 ⇤ 10 3 , KH = 1 ⇤ 10 7 . JGLU H = e ][He ] Tmax ( [GLU KGLU KH [GLUm ][Hm ] ) KGLU KH e] m] e] m] e ][He ] m ][Hm ] 1 + [GLU + [GLU + [H + [H + [GLU + [GLU KGLU KGLU KH KH KGLU KH KGLU KH (16) The binding constants in molar are KGLU = 1.4 ⇤ 10 3 , KH = 1 ⇤ 10 7 JDCC(M AL) = m ][P ie ] Tmax ( [MKAL M AL KP i [P im ][M ALe ] ) KM AL KP i Y (17) ie ] ALe ] ALm ] CCe ] CCm ] e ][P im ] m ][P ie ] For Y = 1 + [P + [PKiPmi ] + [M + [M + [SU + [SU + [MKAL + [MKAL + KP i KM AL KM AL KSU CC KSU CC M AL KP i M AL KP i [SU CCe ][P im ] CCm ][P ie ] + [SU KSU CC KP i KSU CC KP i The binding constants in molar are KM AL = 1.17 ⇤ 10 2 , KP i = 0.93 ⇤ 10 3 JDCC(SU CC) = CCm ][P ie ] Tmax ( [SU KSU CC KP i Y [P im ][SU CCe ] ) KSU CC KP i (18) ie ] CCe ] CCm ] ALe ] ALm ] CCe ][P im ] CCm ][P ie ] For Y = 1 + [P + [PKiPmi ] + [SU + [SU + [M + [M + [SU + [SU + KP i KSU CC KSU CC KM AL KM AL KSU CC KP i KSU CC KP i [M ALe ][P im ] m ][P ie ] + [MKAL KM AL KP i M AL KP i 35 The binding constants in molar are KSU CC = 1 ⇤ 10 3 , KP i = .93 ⇤ 10 4 JT CC = e ][Hm ][CITm ] Tmax ( [MKAL M AL KH KCIT [CITe ][M ALm ][He ] ) KM AL KH KCIT (19) Y ALm ] [M ALe ] [CITm ][Hm ][M ALe ] [CITe ][He ][M ALm ] m ][Hm ] e ][He ] m ][Hm ] For Y = 1+ [CIT + [CIT + [CIT + [M + KM AL + KM AL KH KCIT + KM AL KH KCIT KCIT KH KCIT KH KCIT KH KM AL The binding constants in molar are KCIT = 1 ⇤ 10 3 , KM AL = 0.25 ⇤ 10 3 , KH = 1 ⇤ 10 7 e ][AKGm ] Tmax ( [MKAL M AL KAKG [AKGe ][M ALm ] ) KM AL KAKG JOM E = ALm ] ALe ] AKGe m ][M ALe ] m ][AKGe ] m 1 + [M + [M + AKG +K + [AKG + [MKAL KM AL KM AL KAKG KM AL KAKG AKG AKG KM AL (20) The binding constants in molar are KAKG = 0.24 ⇤ 10 3 , KM AL = 1 ⇤ 10 3 JGAE = e ][Hm ][GLUm ] Tmax (exp( 0.5F ) [ASP RT KASP KH KGLU m ][He ][GLUe ] exp( 0.5F ) [ASP )) RT KASP KH KGLU ASPe e ][GLUe ] m] e ][He ][ASPm ] m ][Hm ][ASPe ] m 1 + [H + [HKmH][GLU + ASP +K + [GLU + [GLU KH KGLU KGLU KASP KGLU KH KASP KGLU KH KASP ASP (21) The binding constants in Molar are KASP = 0.12 ⇤ 10 3 , KGLU = 0.25 ⇤ 10 3 , KH = 1 ⇤ 10 7 JAN T = TF Tmax (exp( ANRT e ][AT Pm ] ) [ADP KADP KAT P Pm ] m] TF (1 + [ADP + [AT exp( ANRT KADP KAT P exp( ( AN TRT1)F m ][AT Pe ] )) [ADP ) KADP KAT P Pe ] e] )(1 + [ADP + [AT exp( ( AN TRT1)F KADP KAT P )) (22) The binding constants in molar are KADP = 4.08 ⇤ 10 6 , KAT P = 5.7 ⇤ 10 6 + JP IC = + [P im ][Hm ] ) KP i KH e ] Tmax ( [PKiPe ][H i KH + + + + ie ] ][Hm ] e ] m] e ] 1 + [P + [PKiPmi ] + [H + [H + [PKiPe ][H + [PKim KP i KH KH i KH P i KH (23) The binding constants in molar are KP I = 9.4 ⇤ 10 3 , KH = 1 ⇤ 10 7 JHLEAK = Tmax F ) ([Hi+ ]exp( KH 2RT The binding constant in molar is KH = 1 ⇤ 10 7 36 + [Hm ])exp( F ) 2RT (24) 8 Appendix B Below are the governing ODEs used in the computational model. Bu↵er Region ODEs Nucleotides d[ADPe ] = JAN T dt d[AT Pe ] = JAN T Ve dt Ve (8.1) (8.2) Substrates d[P ie ] = JP IC + JDCC(M AL) + JDCC(SU CC) dt d[P Y Re ] = JP Y RH Ve dt d[M ALe ] = JOM E JDCC(M AL) JT CC Ve dt d[CITe ] = JT CC Ve dt d[AKGe ] = JOM E Ve dt d[SU CCe ] = JDCC(SU CC) Ve dt d[GLUe ] = JGLU H + JGAE Ve dt d[ASPe ] = JGAE Ve dt Ve (8.3) (8.4) (8.5) (8.6) (8.7) (8.8) (8.9) (8.10) Ions d[He+ ] =0 dt (8.11) 37 Mitochondrial Matrix Nucleotides d[P im ] = JP IC JSCAS JCV JDCC(M AL) JDCC(SU CC) dt d[AT Pm ] = JAN T + JN DK + JCV Vm dt d[ADPm ] d[AT Pm ] = Vm Vm dt dt d[GT Pm ] = JSCAS JN DK Vm dt d[GDPm ] d[GT Pm ] = Vm Vm dt dt Vm (8.12) (8.13) (8.14) (8.15) (8.16) Respiratory Substrates d[N ADHm ] = JP DH + JCIT DH + JAKGDH + JM DH dt + d[N ADm ] d[N ADHm ] Vm = Vm dt dt d[U Qm ] = JCI JCII + JCIII Vm dt d[U QH2m ] d[U Qm ] = Vm Vm dt dt d[CY T Com ] = 2( JCIII + JCIV ) Vm dt d[CY T Crm ] d[CY T Com ] = Vm Vm dt dt Vm 38 JCI (8.17) (8.18) (8.19) (8.20) (8.21) (8.22) TCA Cycle Substrates d[GLUm ] = JGOT JGAE + JGLU H dt d[ASPm ] = JGOT + JGAE Vm dt d[P Y Rm ] = JP DH + JP Y RH Vm dt d[OXAm ] = JM DH JCIT S + JGOT Vm dt d[CITm ] = JCIT S JCIT DH JT CC Vm dt d[AKGm ] = JCIT S JAKGDH JGOT JOM E Vm dt d[SCAm ] = JAKGDH JSCAS Vm dt d[SU CCm ] = JSCAS JSU CDH + JDCC(SU C) Vm dt d[F U Mm ] = JSDH JF H Vm dt d[M ALm ] = JF H JM DH + JDCC(DCC(M AL) + JT CC + JOM E Vm dt d[ACOAm ] = JP DH JCIT S Vm dt d[COAm ] = JP DH + JCIT S + JSCAS JAKGDH Vm dt Vm (8.23) (8.24) (8.25) (8.26) (8.27) (8.28) (8.29) (8.30) (8.31) (8.32) (8.33) (8.34) Oxygen Consumption (Vm + Ve + Vi ) d[O2 ] = dt 0.5JCIV (8.35) Membrane Potential d m dt = 1 Cimm (4JCI + 2JCIII + 4JCIV 3JCV 39 JHLEAK JAN T ) (8.36) Proton Leak d m dt = 1 Cimm (4JCI + 2JCIII + 4JCIV 3JCV 40 JHLEAK JAN T ) (8.37) 9 Appendix C Below is the MATLAB code for the model simulation. The first section is the main body code adapted from our collaborators used to produce the simulations. Additional code needed to run the simulation from our collaborators was used but not altered [15]. The sensitivity analysis code that was used in combination with the main body code is listed next and then lastly is the least squares approximation section of the code. % Load parameters for OCR simulation . clear all close all clc addpath ( ' ./ enzymes ') lsqV = readmatrix ( ' lsqVmax . xlsx ') ; % load delta Vmax values from least squares soultion . test_data2 = readmatrix ( ' 030923 _WT_JO2_cropped . xlsx ') ; % Load UOFU experimental data WT . test_data3 = readmatrix ( ' 030923 _KO_JO2_cropped . xlsx ') ; % Load UOFU experimental data Prdm16 KO . load Model_Params % Loads constants for model . run Read_Condt run Read_IndexSize % lsq delta vmax , different lsq for KO vs WT . p . pest = p . pest + lsqV ; % start of simulation information and loop . 41 p . Ve =.5*.4*20 e -3*4; % buffer_Vol / mito_Vol for mito .025 mg / mL p . time =[6 , 3 , 1 , 4 , 10] '; % addition sequence : Mito , Sub , Rot , ADP , S4 in minutes . %% load param and data % data coloumns indicate : [ PM GM AM Suc SucRot ] JO2_1ADPd =[ test_data2 (: ,1) test_data2 (: ,2) test_data2 (: ,3) test_data2 (: ,4) test_data2 (: ,5) ]; % nM JO2_1ADPKO =[ test_data3 (: ,1) test_data3 (: ,2) test_data3 (: ,3) test_data3 (: ,4) test_data3 (: ,5) ]; %% substrates and ADP additions PYR_index =[1 0 0 0 0]; GLU_index =[0 1 0 0 0]; aKG_index =[0 0 1 0 0]; MAL_index =[1 1 1 0 0]; SUC_index =[0 0 0 1 1]; ADP_add =[0 300]*1 e -6; % uM p . ADPL = length ( ADP_add ) ; options = odeset ( ' NonNegative ' ,[1: p . NOde ]) ; % concentrations of the state variables should be positive for i = p . ISub :1: p . NSub X0 = ICs ( p ) ; p_tem = p ; 42 % %% solving ODEs and calculating state variables T0 =0; jj =1; cc =1; p . Es =3; % extra states for ii =1:1: p . ADPL + p . Es % %% Substrate addition if ii == p . Es -1 X0 ( p_tem . iMALe ) = cc * MAL_index ( i ) *2.5 e -3; % mM X0 ( p_tem . iPYRe ) = cc * PYR_index ( i ) *5 e -3; % mM X0 ( p_tem . iGLUe ) = cc * GLU_index ( i ) *5 e -3; % mM X0 ( p_tem . iaKGe ) = cc * aKG_index ( i ) *5 e -3; % mM X0 ( p_tem . iSUCe ) = cc * SUC_index ( i ) *10 e -3; % mM end if i ==5 && ii == p . Es % Rot addition p_tem . ini_VTmax ( p . iCI ) =0* p_tem . ini_VTmax ( p . iCI ) ; % inhibit CI end % %% ADP addition if ii >= p . Es && ii ~= p . ADPL + p . Es X0 ( p_tem . iADPe ) = X0 ( p . iADPe ) + ADP_add ( jj ) ; jj = jj +1; end % %% Solving ODEs tspan =[ T0 : p . tstep :( T0 + p . time ( ii ,1) ) ]; [T , X ] = ode15s (@ ODEs , tspan , X0 , options , p_tem ) ; 43 T0 = T ( end ,:) ; X0 = X ( end ,:) ; % redefining initial values for the next time periode Tc ( ii , i ) ={ T }; Xc ( ii , i ) ={ X }; CO2 ( ii , i ) =1 e3 * X ( end , p . iO2m ) ; % %% Calculating fluxes for zz =1:1: length ( T ) % length ( tspan ) % length ( T ) J ( zz ,: , i ) = Fluxes ( X ( zz ,:) , p_tem ) ; end Jc ( ii , i ) ={ J (1: zz ,: , i ) }; % %% finding state 3 peaks of all states variables ( ODEs ) if ii < p . Es || ii == p . ADPL + p . Es XPk ( ii ,: , i ) = mean ( Xc { ii , i }) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es XPk ( ii ,: , i ) = max ( Xc { ii , i }) ; Xmin ( ii ,: , i ) = min ( Xc { ii , i }) ; end % %% finding state 3 peaks of fluxes st =2; for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es JPk ( ii ,s , i ) = mean ( J (2: end ,s , i ) ) ; each column HPk ( ii , i ) = mean ( J (2: end ,s , i ) ) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es 44 % returns mean of JPk ( ii ,s , i ) = max ( J ( st : end ,s , i ) ) ; % returns max of each column HPk ( ii , i ) =.5* max ( J ( st : end ,s , i ) ) ; end end % %% finding length of state 3 for fluxes for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es Ix ( ii , i ) =0; Ixb ( ii , i ) =0; Ixf ( ii , i ) =0; LS3 ( ii , i ) =0; elseif ii >= p . Es || ii ~= p . ADPL + p . Es if isempty ( find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ' ) ) ; Ix ( ii , i ) =0; wp = sprintf ( ' cannot find Ix for org =% d ii =% d flx =% d ' ,i , ii , s ) ; disp ( wp ) ; else Ix ( ii , i ) = find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ') ; end % find the max index end end end % ii time protocol for - loop % storing cell variables in vectors Tv (: , i ) = }; [ Tc {1 , i }; Tc {3 , i }; Tc {4 , i Xc {2 , i }; Xc {3 , i }; Xc {4 , i Tc {5 , i };]; Xv (: ,: , i ) =1 e0 * }; Tc {2 , i }; [ Xc {1 , i }; Xc {5 , i };]; % M 45 Jv (: ,: , i ) =1 e9 * }; [ Jc {1 , i }; Jc {2 , i }; Jc {3 , i }; Jc {4 , i Jc {5 , i };]; % nmol / min / mg mito p . sp =3; % start of plot OCR (: , i ) = movmean ( Jv ( p . sp : end , p . iCIV , i ) ,40) ; % 5 point average slope ??? end %% Plot % %% OCR fluxes % %%% PYR + MAL Plot %%%%%%% figure (100) ; subplot (3 ,2 ,1) ; % OCR dynamic p . sp =3; % start of plot OCR (: ,1) = movmean ( Jv ( p . sp : end , p . iCIV ,1) ,40) ; plot ( timev2 (1: length ( JO2_1ADPd ( p . sp : end ,1) ) ) , JO2_1ADPd ( p . sp : end ,1) , ' DisplayName ' , ' Experimental WT ' , ' linewidth ' , linewidth ) ; hold on plot ( Tv ( p . sp : end ,1) , OCR (: ,1) , ' DisplayName ' , ' Simulation ' , ' linewidth ' , linewidth ) ; hold off title ( ' PYR + MAL ') ; xlabel ( ' Time ( min ) ') ; ylabel ({ ' J_O_2 ( nmol / min / mg ) ' }) ; set ( gca , ' Fontsize ' ,15) ; lgd1 = legend ; legend boxoff ; lgd1 . NumColumns = 1; xlim ([0 40]) ; xticks ([0:5:40]) 46 ylim ([0 800]) ; yticks ([0:200:800]) % %%% GLU + MAL Plot %%%%%%% subplot (3 ,2 ,2) ; p . sp =3; % start of plot OCR (: ,2) = movmean ( Jv ( p . sp : end , p . iCIV ,2) ,40) ; plot ( timev2 (1: length ( JO2_1ADPd ( p . sp : end ,1) ) ) , JO2_1ADPd ( p . sp : end ,2) , ' DisplayName ' , ' Experimental WT ' , ' linewidth ' , linewidth ) ; hold on plot ( Tv ( p . sp : end ,2) , OCR (: ,2) , ' DisplayName ' , ' Simulation ' , ' linewidth ' , linewidth ) ; title ( ' GLU + MAL ') ; xlabel ( ' Time ( min ) ') ; ylabel ({ ' J_O_2 ( nmol / min / mg ) ' }) ; set ( gca , ' Fontsize ' , 15) ; lgd1 = legend ; legend boxoff ; lgd1 . NumColumns = 1; xlim ([0 40]) ; xticks ([0:5:40]) ylim ([0 800]) ; yticks ([0:200:800]) % %%% aKG + MAL Plot %%%%%%% subplot (3 ,2 ,3) ; p . sp =3; % start of plot OCR (: ,3) = movmean ( Jv ( p . sp : end , p . iCIV ,3) ,40) ; plot ( timev2 (1: length ( JO2_1ADPd ( p . sp : end ,1) ) ) , JO2_1ADPd ( p . sp : end ,3) , ' DisplayName ' , ' Experimental WT ' , ' linewidth ' , 47 linewidth ) ; hold on plot ( Tv ( p . sp : end ,3) , OCR (: ,3) , ' DisplayName ' , ' Simulation ' , ' linewidth ' , linewidth ) ; title ( ' aKG + MAL ') ; xlabel ( ' Time ( min ) ') ; ylabel ({ ' J_O_2 ( nmol / min / mg ) ' }) ; set ( gca , ' Fontsize ' , 15) ; lgd1 = legend ; legend boxoff ; lgd1 . NumColumns = 1; xlim ([0 40]) ; xticks ([0:5:40]) ylim ([0 800]) ; yticks ([0:200:800]) % %%% Succ + Rot Plot %%%%%%% subplot (3 ,2 ,4) ; p . sp =3; % start of plot OCR (: ,5) = movmean ( Jv ( p . sp : end , p . iCIV ,5) ,40) ; plot ( timev2 (1: length ( JO2_1ADPd ( p . sp : end ,1) ) ) , JO2_1ADPd ( p . sp : end ,5) , ' DisplayName ' , ' Experimental WT ' , ' linewidth ' , linewidth ) ; hold on plot ( Tv ( p . sp : end ,5) , OCR (: ,5) , ' DisplayName ' , ' Simulation ' , ' linewidth ' , linewidth ) ; title ( ' SUCC + ROT ') ; xlabel ( ' Time ( min ) ') ; ylabel ({ ' J_O_2 ( nmol / min / mg ) ' }) ; set ( gca , ' Fontsize ' , 15) ; lgd1 = legend ; legend boxoff ; 48 lgd1 . NumColumns = 1; xlim ([0 40]) ; xticks ([0:5:40]) ylim ([0 800]) ; yticks ([0:200:800]) sgtitle ( ' text ' , ' Fontsize ' , 22 , ' Linewidth ' , 2) ; sgtitle ( ' Simulated vs Experimental Data Parameterization ' , ' Fontsize ' , 22) ; This code below was used to generate the base of the sensitivity matrix. The sensitivity matrix was completed and further analyzed in excel. % %%%%%% Load Data and time matrcies %%%%%%%%%%% clear all close all addpath ( ' ./ enzymes ') timev2 = readmatrix ( ' timev2 . xlsx ') ; WT_exp_data = readmatrix ( ' WTcropextend . xlsx ') ; load Model_Params run Read_Condt run Read_IndexSize WMC_exp_data = data_HT . JO2_1ADP_nmol ; % %%%%%%% constructing A matrix %%%%%%%% % A matrix % get original simulation values p . Ve =.5*.4*20 e -3*4; % buffer_Vol / mito_Vol for mito .1 mg / mL Boudina uses .025 mg / mL 49 %% p . time =[6 , 3 , 1 , 4 , 10] '; % addition sequence : mito , Sub , Rot , ADP , S4 % %% For more info on the variables declaration see Read_Condt %% substrates and ADP additions PYR_index =[1 0 0 0 0]; GLU_index =[0 1 0 0 0]; aKG_index =[0 0 1 0 0]; MAL_index =[1 1 1 0 0]; SUC_index =[0 0 0 1 1]; ADP_add =[0 300]*1 e -6; % uM p . ADPL = length ( ADP_add ) ; options = odeset ( ' NonNegative ' ,[1: p . NOde ]) ; % concentrations of the state variables should be positive for i = p . ISub :1: p . NSub X0 = ICs ( p ) ; p_tem = p ; % %% solving ODEs and calculating state variables T0 =0; jj =1; cc =1; p . Es =3; % extra states for ii =1:1: p . ADPL + p . Es % %% Substrate addition if ii == p . Es -1 X0 ( p_tem . iMALe ) = cc * MAL_index ( i ) *2.5 e -3; % mM X0 ( p_tem . iPYRe ) = cc * PYR_index ( i ) *5 e -3; % mM 50 X0 ( p_tem . iGLUe ) = cc * GLU_index ( i ) *5 e -3; % mM X0 ( p_tem . iaKGe ) = cc * aKG_index ( i ) *5 e -3; % mM X0 ( p_tem . iSUCe ) = cc * SUC_index ( i ) *10 e -3; % mM end if i ==5 && ii == p . Es % Rot addition p_tem . ini_VTmax ( p . iCI ) =0* p_tem . ini_VTmax ( p . iCI ) ; % inhibit CI end % %% ADP addition if ii >= p . Es && ii ~= p . ADPL + p . Es X0 ( p_tem . iADPe ) = X0 ( p . iADPe ) + ADP_add ( jj ) ; jj = jj +1; end % %% Solving ODEs tspan =[ T0 : p . tstep :( T0 + p . time ( ii ,1) ) ]; [T , X ] = ode15s (@ ODEs , tspan , X0 , options , p_tem ) ; T0 = T ( end ,:) ; X0 = X ( end ,:) ; % redefining initial values for the next time periode Tc ( ii , i ) ={ T }; Xc ( ii , i ) ={ X }; CO2 ( ii , i ) =1 e3 * X ( end , p . iO2m ) ; % %% Calculating fluxes for zz =1:1: length ( T ) % length ( tspan ) % length ( T ) J ( zz ,: , i ) = Fluxes ( X ( zz ,:) , p_tem ) ; end Jc ( ii , i ) ={ J (1: zz ,: , i ) }; 51 % %% finding state 3 peaks of all states variables ( ODEs ) if ii < p . Es || ii == p . ADPL + p . Es XPk ( ii ,: , i ) = mean ( Xc { ii , i }) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es XPk ( ii ,: , i ) = max ( Xc { ii , i }) ; Xmin ( ii ,: , i ) = min ( Xc { ii , i }) ; end % %% finding state 3 peaks of fluxes st =2; for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es JPk ( ii ,s , i ) = mean ( J (2: end ,s , i ) ) ; % returns mean of each column HPk ( ii , i ) = mean ( J (2: end ,s , i ) ) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es JPk ( ii ,s , i ) = max ( J ( st : end ,s , i ) ) ; % returns max of each column HPk ( ii , i ) =.5* max ( J ( st : end ,s , i ) ) ; end end % %% finding length of state 3 for fluxes for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es Ix ( ii , i ) =0; Ixb ( ii , i ) =0; Ixf ( ii , i ) =0; LS3 ( ii , i ) =0; 52 elseif ii >= p . Es || ii ~= p . ADPL + p . Es if isempty ( find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ') ) ; Ix ( ii , i ) =0; wp = sprintf ( ' cannot find Ix for org =% d ii =% d flx =% d ' ,i , ii , s ) ; disp ( wp ) ; else Ix ( ii , i ) = find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ') ; end % find the max index end end end % ii time protocol for - loop % storing cell variables in vectors Tv (: , i ) = [ Tc {1 , i }; Tc {2 , i }; Tc {3 , i }; Tc {4 , i }; [ Xc {1 , i }; Xc {2 , i }; Xc {3 , i }; Xc {4 , i }; Jc {2 , i }; Jc {3 , i }; Jc {4 , i }; Tc {5 , i };]; Xv (: ,: , i ) =1 e0 * Xc {5 , i };]; % M Jv (: ,: , i ) =1 e9 * [ Jc {1 , i }; Jc {5 , i };]; % nmol / min / mg mito p . sp =3; % start of plot OCR (: , i ) = movmean ( Jv ( p . sp : end , p . iCIV , i ) ,40) ; end OCR_original = OCR ; OCR_original (: ,4) = []; stacker_OCR_orig = OCR_original (:) ; 53 original_Vmax = p . pest ; clearvars - except delta_OCR_exp original_Vmax stacker_OCR_orig stackerWT_exp % populate rest of A matrix with adjusted vmax for ss = 1:24 load Model_Params run Read_Condt run Read_IndexSize delta_vmax = 4; p . pest ( ss ) = p . pest ( ss ) * delta_vmax ; p . Ve =.5*.4*20 e -3*4; % buffer_Vol / mito_Vol for mito .025 mg / mL p . time =[6 , 3 , 1 , 4 , 10] '; % addition sequence : mito , Sub , Rot , ADP , S4 % %% For more info on the variables declaration see Read_Condt %% substrates and ADP additions PYR_index =[1 0 0 0 0]; GLU_index =[0 1 0 0 0]; 54 aKG_index =[0 0 1 0 0]; MAL_index =[1 1 1 0 0]; SUC_index =[0 0 0 1 1]; ADP_add =[0 300]*1 e -6; % uM p . ADPL = length ( ADP_add ) ; options = odeset ( ' NonNegative ' ,[1: p . NOde ]) ; % concentrations of the state variables should be positive for i = p . ISub :1: p . NSub X0 = ICs ( p ) ; p_tem = p ; % %% solving ODEs and calculating state variables T0 =0; jj =1; cc =1; p . Es =3; % extra states for ii =1:1: p . ADPL + p . Es % %% Substrate addition if ii == p . Es -1 X0 ( p_tem . iMALe ) = cc * MAL_index ( i ) *2.5 e -3; % mM X0 ( p_tem . iPYRe ) = cc * PYR_index ( i ) *5 e -3; % mM X0 ( p_tem . iGLUe ) = cc * GLU_index ( i ) *5 e -3; % mM X0 ( p_tem . iaKGe ) = cc * aKG_index ( i ) *5 e -3; % mM X0 ( p_tem . iSUCe ) = cc * SUC_index ( i ) *10 e -3; % mM end if i ==5 && ii == p . Es % Rot addition p_tem . ini_VTmax ( p . iCI ) =0* p_tem . ini_VTmax ( p . iCI ) ; % inhibit CI end 55 % %% ADP addition if ii >= p . Es && ii ~= p . ADPL + p . Es X0 ( p_tem . iADPe ) = X0 ( p . iADPe ) + ADP_add ( jj ) ; jj = jj +1; end % %% Solving ODEs tspan =[ T0 : p . tstep :( T0 + p . time ( ii ,1) ) ]; [T , X ] = ode15s (@ ODEs , tspan , X0 , options , p_tem ) ; T0 = T ( end ,:) ; X0 = X ( end ,:) ; % redefining initial values for the next time periode Tc ( ii , i ) ={ T }; Xc ( ii , i ) ={ X }; CO2 ( ii , i ) =1 e3 * X ( end , p . iO2m ) ; % %% Calculating fluxes for zz =1:1: length ( T ) % length ( tspan ) % length ( T ) J ( zz ,: , i ) = Fluxes ( X ( zz ,:) , p_tem ) ; end Jc ( ii , i ) ={ J (1: zz ,: , i ) }; % %% finding state 3 peaks of all states variables ( ODEs ) if ii < p . Es || ii == p . ADPL + p . Es XPk ( ii ,: , i ) = mean ( Xc { ii , i }) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es XPk ( ii ,: , i ) = max ( Xc { ii , i }) ; Xmin ( ii ,: , i ) = min ( Xc { ii , i }) ; end 56 % %% finding state 3 peaks of fluxes st =2; for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es JPk ( ii ,s , i ) = mean ( J (2: end ,s , i ) ) ; % returns mean of each column HPk ( ii , i ) = mean ( J (2: end ,s , i ) ) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es JPk ( ii ,s , i ) = max ( J ( st : end ,s , i ) ) ; % returns max of each column HPk ( ii , i ) =.5* max ( J ( st : end ,s , i ) ) ; end end % %% finding length of state 3 for fluxes for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es Ix ( ii , i ) =0; Ixb ( ii , i ) =0; Ixf ( ii , i ) =0; LS3 ( ii , i ) =0; elseif ii >= p . Es || ii ~= p . ADPL + p . Es if isempty ( find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ' ) ) ; Ix ( ii , i ) =0; wp = sprintf ( ' cannot find Ix for org =% d ii =% d flx =% d ' ,i , ii , s ) ; disp ( wp ) ; else Ix ( ii , i ) = find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ') ; end % find the max index end end 57 end % ii time protocol for - loop % storing cell variables in vectors Tv (: , i ) = }; [ Tc {1 , i }; Tc {4 , i [ Xc {1 , i }; Xc {2 , i }; Xc {3 , i }; Xc {4 , i Jc {2 , i }; Jc {3 , i }; Jc {4 , i Xc {5 , i };]; % M Jv (: ,: , i ) =1 e9 * }; Tc {3 , i }; Tc {5 , i };]; Xv (: ,: , i ) =1 e0 * }; Tc {2 , i }; [ Jc {1 , i }; Jc {5 , i };]; % nmol / min / mg mito p . sp =3; % start of plot OCR (: , i ) = movmean ( Jv ( p . sp : end , p . iCIV , i ) ,40) ; end OCR (: ,4) = []; B = ones (1443 - length ( OCR (: ,4) ) ,4) ; for r =1:4 if length ( OCR (: , r ) ) ~= 1443 k = OCR (: , r ) ; f = k ( end ) ; L (r ,1) = f ; e = 1443 - length ( OCR (: , r ) ) ; B (: , r ) = B (: , r ) * L (r ,1) ; end 58 end OCR2 = vertcat ( OCR , B ) ; stacker_OCR = OCR2 (:) ; A_ocr { ss }= stacker_OCR ; clearvars - except delta_vmax delta_OCR_exp original_Vmax ss A_ocr stacker_OCR_orig stackerWT_exp B end C = zeros ([5772 24]) ; for bb =1:24 C (: , bb ) = A_ocr {1 , bb }; end orig_stack_24 = zeros (5772 ,24) ; for bb =1:24 A (: , bb ) = stacker_OCR_orig ; end delta_OCR = A - C ; delta_OCR = delta_OCR ./ A ; delta_v = original_Vmax / delta_vmax ; Least squares approximation code. % %%%%%% Load Data and time matrcies %%%%%%%%%%% clear all close all addpath ( ' ./ enzymes ') timev2 = readmatrix ( ' timev2 . xlsx ') ; WT_exp_data = readmatrix ( ' 030923 _WT_JO2_cropped . xlsx ') ; load Model_Params 59 run Read_Condt run Read_IndexSize % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%% Interpolate data to match length for simulation and our data %%%%%%% OCR_size_holder = zeros (1443 ,5) ; q = size ( WT_exp_data ,1) ; w = size ( OCR_size_holder ,1) ; interp_WT_exp = zeros (w ,5) ; x2i = (1: w ) ' ; x2 = linspace (1 ,w , q ) ' ; for i = 1:5 A2i = interp1 ( x2 , WT_exp_data (: , i ) , x2i ) ; interp_WT_exp (: , i ) = A2i ; end % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % delta OCR matrix ( B matrix ) % remove succinate only and stack substrated . interp_WT_exp (: ,4) = []; stacked_WT_exp = interp_WT_exp (:) ; 60 % %%%%%%% constructing A matrix %%%%%%%% % A matrix % get initial simulation values based on aproximation from sensitivity analysis . p . pest (23) = p . pest (23) *.60*.60*.60*.8; p . pest (10) = p . pest (10) *.010; PM p . pest (13) = p . pest (13) *10; p . pest (4) = p . pest (4) *.015; p . pest (22) = p . pest (22) *.97; p . pest (18) = p . pest (18) *.96; p . Ve =.5*.4*20 e -3*4; % buffer_Vol / mito_Vol for mito .025 mg / mL p . time =[6 , 3 , 1 , 4 , 10] '; % addition sequence : mito , Sub , Rot , ADP , S4 % %% For more info on the variables declaration see Read_Condt %% substrates and ADP additions PYR_index =[1 0 0 0 0]; GLU_index =[0 1 0 0 0]; aKG_index =[0 0 1 0 0]; MAL_index =[1 1 1 0 0]; SUC_index =[0 0 0 1 1]; 61 ADP_add =[0 300]*1 e -6; % uM p . ADPL = length ( ADP_add ) ; options = odeset ( ' NonNegative ' ,[1: p . NOde ]) ; % concentrations of the state variables should be positive for i = p . ISub :1: p . NSub X0 = ICs ( p ) ; p_tem = p ; % %% solving ODEs and calculating state variables T0 =0; jj =1; cc =1; p . Es =3; % extra states for ii =1:1: p . ADPL + p . Es % %% Substrate addition if ii == p . Es -1 X0 ( p_tem . iMALe ) = cc * MAL_index ( i ) *2.5 e -3; % mM X0 ( p_tem . iPYRe ) = cc * PYR_index ( i ) *5 e -3; % mM X0 ( p_tem . iGLUe ) = cc * GLU_index ( i ) *5 e -3; % mM X0 ( p_tem . iaKGe ) = cc * aKG_index ( i ) *5 e -3; % mM X0 ( p_tem . iSUCe ) = cc * SUC_index ( i ) *10 e -3; % mM end if i ==5 && ii == p . Es % Rot addition p_tem . ini_VTmax ( p . iCI ) =0* p_tem . ini_VTmax ( p . iCI ) ; % inhibit CI end % %% ADP addition if ii >= p . Es && ii ~= p . ADPL + p . Es X0 ( p_tem . iADPe ) = X0 ( p . iADPe ) + ADP_add ( jj ) ; 62 jj = jj +1; end % %% Solving ODEs tspan =[ T0 : p . tstep :( T0 + p . time ( ii ,1) ) ]; [T , X ] = ode15s (@ ODEs , tspan , X0 , options , p_tem ) ; T0 = T ( end ,:) ; X0 = X ( end ,:) ; % redefining initial values for the next time periode Tc ( ii , i ) ={ T }; Xc ( ii , i ) ={ X }; CO2 ( ii , i ) =1 e3 * X ( end , p . iO2m ) ; % %% Calculating fluxes for zz =1:1: length ( T ) % length ( tspan ) % length ( T ) J ( zz ,: , i ) = Fluxes ( X ( zz ,:) , p_tem ) ; end Jc ( ii , i ) ={ J (1: zz ,: , i ) }; % %% finding state 3 peaks of all states variables ( ODEs ) if ii < p . Es || ii == p . ADPL + p . Es XPk ( ii ,: , i ) = mean ( Xc { ii , i }) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es XPk ( ii ,: , i ) = max ( Xc { ii , i }) ; Xmin ( ii ,: , i ) = min ( Xc { ii , i }) ; end % %% finding state 3 peaks of fluxes st =2; for s =1: p . NPar 63 if ii < p . Es || ii == p . ADPL + p . Es JPk ( ii ,s , i ) = mean ( J (2: end ,s , i ) ) ; % returns mean of each column HPk ( ii , i ) = mean ( J (2: end ,s , i ) ) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es JPk ( ii ,s , i ) = max ( J ( st : end ,s , i ) ) ; % returns max of each column HPk ( ii , i ) =.5* max ( J ( st : end ,s , i ) ) ; end end % %% finding length of state 3 for fluxes for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es Ix ( ii , i ) =0; Ixb ( ii , i ) =0; Ixf ( ii , i ) =0; LS3 ( ii , i ) =0; elseif ii >= p . Es || ii ~= p . ADPL + p . Es if isempty ( find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ') ) ; Ix ( ii , i ) =0; wp = sprintf ( ' cannot find Ix for org =% d ii =% d flx =% d ' ,i , ii , s ) ; disp ( wp ) ; else Ix ( ii , i ) = find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ') ; end % find the max index end end end % ii time protocol for - loop % storing cell variables in vectors 64 Tv (: , i ) = [ Tc {1 , i }; Tc {2 , i }; Tc {3 , i }; Tc {4 , i }; [ Xc {1 , i }; Xc {2 , i }; Xc {3 , i }; Xc {4 , i }; Jc {2 , i }; Jc {3 , i }; Jc {4 , i }; Tc {5 , i };]; Xv (: ,: , i ) =1 e0 * Xc {5 , i };]; % M Jv (: ,: , i ) =1 e9 * [ Jc {1 , i }; Jc {5 , i };]; % nmol / min / mg mito p . sp =3; % start of plot OCR (: , i ) = movmean ( Jv ( p . sp : end , p . iCIV , i ) ,40) ; end OCR_original = OCR ; OCR_original (: ,4) = []; stacked_OCR_orig = OCR_original (:) ; original_Vmax = p . pest ; clearvars - except delta_OCR_exp original_Vmax stacked_OCR_orig stacked_WT_exp % populate rest of A matrix with adjusted vmax for ss = 1:24 load Model_Params run Read_Condt run Read_IndexSize delta_vmax = .1; 65 p . pest ( ss ) = p . pest ( ss ) * delta_vmax ; p . Ve =.5*.4*20 e -3*4; % buffer_Vol / mito_Vol for mito .025 mg / mL p . time =[6 , 3 , 1 , 4 , 10] '; % addition sequence : mito , Sub , Rot , ADP , S4 % %% For more info on the variables declaration see Read_Condt %% substrates and ADP additions PYR_index =[1 0 0 0 0]; GLU_index =[0 1 0 0 0]; aKG_index =[0 0 1 0 0]; MAL_index =[1 1 1 0 0]; SUC_index =[0 0 0 1 1]; ADP_add =[0 300]*1 e -6; % uM p . ADPL = length ( ADP_add ) ; options = odeset ( ' NonNegative ' ,[1: p . NOde ]) ; % concentrations of the state variables should be positive for i = p . ISub :1: p . NSub X0 = ICs ( p ) ; p_tem = p ; % %% solving ODEs and calculating state variables T0 =0; jj =1; cc =1; 66 p . Es =3; % extra states for ii =1:1: p . ADPL + p . Es % %% Substrate addition if ii == p . Es -1 X0 ( p_tem . iMALe ) = cc * MAL_index ( i ) *2.5 e -3; % mM X0 ( p_tem . iPYRe ) = cc * PYR_index ( i ) *5 e -3; % mM X0 ( p_tem . iGLUe ) = cc * GLU_index ( i ) *5 e -3; % mM X0 ( p_tem . iaKGe ) = cc * aKG_index ( i ) *5 e -3; % mM X0 ( p_tem . iSUCe ) = cc * SUC_index ( i ) *10 e -3; % mM end if i ==5 && ii == p . Es % Rot addition p_tem . ini_VTmax ( p . iCI ) =0* p_tem . ini_VTmax ( p . iCI ) ; % inhibit CI end % %% ADP addition if ii >= p . Es && ii ~= p . ADPL + p . Es X0 ( p_tem . iADPe ) = X0 ( p . iADPe ) + ADP_add ( jj ) ; jj = jj +1; end % %% Solving ODEs tspan =[ T0 : p . tstep :( T0 + p . time ( ii ,1) ) ]; [T , X ] = ode15s (@ ODEs , tspan , X0 , options , p_tem ) ; T0 = T ( end ,:) ; X0 = X ( end ,:) ; % redefining initial values for the next time periode Tc ( ii , i ) ={ T }; Xc ( ii , i ) ={ X }; 67 CO2 ( ii , i ) =1 e3 * X ( end , p . iO2m ) ; % %% Calculating fluxes for zz =1:1: length ( T ) % length ( tspan ) % length ( T ) J ( zz ,: , i ) = Fluxes ( X ( zz ,:) , p_tem ) ; end Jc ( ii , i ) ={ J (1: zz ,: , i ) }; % %% finding state 3 peaks of all states variables ( ODEs ) if ii < p . Es || ii == p . ADPL + p . Es XPk ( ii ,: , i ) = mean ( Xc { ii , i }) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es XPk ( ii ,: , i ) = max ( Xc { ii , i }) ; Xmin ( ii ,: , i ) = min ( Xc { ii , i }) ; end % %% finding state 3 peaks of fluxes st =2; for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es JPk ( ii ,s , i ) = mean ( J (2: end ,s , i ) ) ; % returns mean of each column HPk ( ii , i ) = mean ( J (2: end ,s , i ) ) ; elseif ii >= p . Es || ii ~= p . ADPL + p . Es JPk ( ii ,s , i ) = max ( J ( st : end ,s , i ) ) ; % returns max of each column HPk ( ii , i ) =.5* max ( J ( st : end ,s , i ) ) ; end 68 end % %% finding length of state 3 for fluxes for s =1: p . NPar if ii < p . Es || ii == p . ADPL + p . Es Ix ( ii , i ) =0; Ixb ( ii , i ) =0; Ixf ( ii , i ) =0; LS3 ( ii , i ) =0; elseif ii >= p . Es || ii ~= p . ADPL + p . Es if isempty ( find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ' ) ) ; Ix ( ii , i ) =0; wp = sprintf ( ' cannot find Ix for org =% d ii =% d flx =% d ' ,i , ii , s ) ; disp ( wp ) ; else Ix ( ii , i ) = find ( J ( st : end ,s , i ) == JPk ( ii , i ) ,1 , ' first ') ; end % find the max index end end end % ii time protocol for - loop % storing cell variables in vectors Tv (: , i ) = }; [ Tc {1 , i }; Tc {4 , i [ Xc {1 , i }; Xc {2 , i }; Xc {3 , i }; Xc {4 , i Jc {2 , i }; Jc {3 , i }; Jc {4 , i Xc {5 , i };]; % M Jv (: ,: , i ) =1 e9 * }; Tc {3 , i }; Tc {5 , i };]; Xv (: ,: , i ) =1 e0 * }; Tc {2 , i }; [ Jc {1 , i }; Jc {5 , i };]; % nmol / min / mg mito p . sp =3; % start of plot OCR (: , i ) = movmean ( Jv ( p . sp : end , p . iCIV , i ) ,40) ; 69 end OCR (: ,4) = []; B = ones (1443 - length ( OCR (: ,4) ) ,4) ; for r =1:4 if length ( OCR (: , r ) ) ~= 1443 k = OCR (: , r ) ; f = k ( end ) ; L (r ,1) = f ; e = 1443 - length ( OCR (: , r ) ) ; B (: , r ) = B (: , r ) * L (r ,1) ; end end OCR2 = vertcat ( OCR , B ) ; stacker_OCR = OCR2 (:) ; A_ocr { ss }= stacker_OCR ; clearvars - except delta_vmax delta_OCR_exp original_Vmax ss A_ocr stacked_OCR_orig stacked_WT_exp B end C = zeros ([5772 24]) ; for bb =1:24 C (: , bb ) = A_ocr {1 , bb }; end 70 orig_stack_24 = zeros (5772 ,24) ; for bb =1:24 A (: , bb ) = stacked_OCR_orig ; end delta_OCR_B = stacked_WT_exp - stacked_OCR_orig ; sim_delta_OCR = C - A ; V = readmatrix ( ' originalV . xlsx ') ; sim_delta_OCR = sim_delta_OCR /.1; w = lsqr ( sim_delta_OCR , delta_OCR_B ,.0000001 , 1000000) ; 71 Name of Candidate: Sophie Stephens Date of Submission: April 24, 2023 72