Using Structural Kinetic Modeling To Identify Key Determinants of

Subscriber access provided by Binghamton University | Libraries

Article

Using Structural Kinetic Modeling to Identify Key Determinants of Stability in Reaction Networks Nicole Jean Carbonaro, and Ian F. Thorpe J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 09 Jun 2017 Downloaded from http://pubs.acs.org on June 11, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Title: Using Structural Kinetic Modeling to Identify Key Determinants of Stability in Reaction Networks Nicole J. Carbonaro and Ian F. Thorpe* Department of Chemistry and Biochemistry, University of Maryland, Baltimore County, Baltimore, MD 21250 USA Corresponding author address: [email protected] Abstract Kinetic modeling is increasingly used to understand the reaction dynamics of metabolic systems. However, one major drawback of kinetic modeling is that appropriate rate parameters required to implement such models are often unavailable. To circumvent this limitation, an approach known as structural kinetic modeling was developed as a way to understand the dynamics of reaction networks without explicitly requiring rate parameters. This study describes a novel approach to use structural kinetic modeling to identify reaction components that contribute most significantly to mediating network stability. We applied this method to analyze the metabolic pathway of glycolysis in yeast. As a result, we identified specific metabolic components that contribute most significantly to defining the stability properties of the glycolysis reaction network and predict the responses of these components to perturbations. These results were validated via comparison to a conventional kinetic model of glycolysis. Thus, applying our approach allows more detailed information about the stability and dynamics of the metabolic network to now be accessible without requiring rate parameters. We anticipate that this method can focus efforts of experimental studies by identifying the susceptibility of reaction components to metabolic engineering. The approach may be applied to a variety of complex reaction networks. 1. Introduction One focal point of chemical kinetics over the last half century has been the modeling of reaction networks. Such models can be used to track changes in concentrations of reactants over time, evaluate network stability, or calculate reaction rates.1-3 Moreover, by perturbing the concentrations of individual components one can see how the concentrations of other components are impacted. Since it is often difficult to measure concentrations in situ in real time, modeling provides a way to understand these relationships in silico. A specific application that has elicited significant interest in recent years is modeling the complex networks of chemical reactions that underlie metabolic processes.4-6 The main methods of analyzing metabolic reaction networks (listed in order of decreasing quantitative character) are kinetic modeling, structural kinetic modeling (SKM) and flux balance analysis (FBA) (see Figure 1). Below we briefly introduce the different approaches to highlight their strengths and limitations. In particular, the advantages of SKM over FBA and kinetic modeling are discussed. In the present study, we employed an approach based on SKM to understand the relative contributions of components of a reaction network to maintaining network stability.

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Methods of metabolic pathway modeling on a qualitative versus quantitative spectrum. Inspired by Figure 1 of Steuer et al.7 1.1. Kinetic modeling Kinetic modeling is one of the oldest and most frequently used methods for determining the dynamics of a system of chemical reactions.8 Detailed descriptions of chemical kinetics and kinetic modeling are widely available.9-14 Kinetic modeling allows one to directly view the concentration of metabolites at all points of time in a simulation. Constructing a kinetic model requires rate parameters dictating the various reaction rates, a reaction scheme to understand the flow of chemical components through the system and initial concentrations of these components.4 There are several examples of kinetic models describing microbe metabolism.4-6, 1516 The rules defining flow of metabolites through the system are recorded in the stoichiometric matrix, which is a key component of all three approaches to modeling described in Figure 1. Unfortunately, for many systems of interest the kinetic rate parameters are not known, preventing the application of kinetic modeling approaches. 1.2. Flux balance analysis More recently flux balance analysis has emerged as an efficient way to describe metabolic flux in large reaction networks.17-22 Flux is a measure of the flow of metabolites through a system and flux balance applies the law of conservation of mass to the reactions in the network. FBA applies steady state conditions, where the flux of each metabolite into the system must equal the flux of that metabolite leaving the system This is not the case in kinetic modeling where steady state is not required. An advantage of FBA is its computational efficiency, which allows the method to readily be applied to an entire cell of an organism. However, a disadvantage of FBA is that it is not possible to determine actual concentrations of metabolites or reaction dynamics since kinetic parameters are not employed.23 FBA approaches have been widely used to study metabolic networks in microorganisms.18-19, 24-26 1.3. Structural Kinetic Modeling There are examples of techniques described as structural kinetic modeling that have been previously reported in the literature. These studies investigated the kinetics of flow of different materials with regard to shear rates and viscosity and had no link to describing metabolic reactions. Malik et al. used the term in their study of the structure and flow of solder pastes,26 while Koocheki and Razavi investigated the relationship between seed gum viscosity and reaction rate.27 These studies are based on applying integrated kinetic rate equations to describe these processes. In contrast, our work will focus on structural kinetic modeling (SKM) used to describe metabolic networks as articulated by Steuer et al.28 in which distinct flux values are

2


Page 2 of 29

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


employed for different network components. SKM allows one to learn about the dynamics of reaction networks without requiring explicit definition of rate constants. This goal is accomplished by defining the rate constants in a parametric way that accurately describes dynamic properties of the network.7, 28-29 Expanding on the pioneering work of Steuer et al.28 in applying SKM to describe glycolysis, we present a novel approach to analyze the Jacobian matrices obtained in SKM in order to identify the components of metabolic networks with the dominant impact on determining network stability. SKM is more quantitative than flux balance analysis since it allows the study of network dynamics. However, it is less quantitative than kinetic modeling, as it does not allow measurement of metabolite concentrations as a function of time. Nonetheless, SKM has the advantage compared to conventional kinetic models that it also does not require rate parameters to be explicitly defined and thus can be applied to a larger variety of metabolic systems. The Jacobian matrix is a central element of SKM. Consider a simplified system of reactions containing three metabolites M1, M2, and M3: M1 M2 M3. The Jacobian matrix for this system contains the partial derivatives of the metabolite fluxes with respect to the concentration of each metabolite. This is defined as: =

(1)

where Vi is the total flux of the ith metabolite and Si the concentration of metabolite i where there are i=1,2,…,m number of metabolites. In the simple reaction scheme above, m=3. The derivatives are evaluated under the assumption that the steady state concentration of each metabolite has a defined value So. The total flux of the ith metabolite can be written in terms of the chemical reactions occurring in a system as: = =

(2) where there are l=1,2,…,r chemical reactions, Nli is an element of the r × m stoichiometric matrix and Rl(S) is the flux in the lth reaction. In general, the reaction fluxes Rl are dependent on the rates and concentrations of metabolites in the respective chemical reactions. Thus, the Jacobian element describing the rate of change of the flux of metabolite Mi with respect to the concentration of another metabolite Mj can be written as: = =

In SKM this equation is represented as follows:

3


(3)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

= Λ Λ = ! " = # " $ % ! " = " =

=

(4-8)

Information about reaction dynamics is obtained by diagonalizing the Jacobian in order to obtain eigenvalues. All eigenvalues with a negative real component denote stable solutions to the combination of fluxes in the reaction network, while any eigenvalue with a positive real component represents an unstable solution. Moreover, eigenvalues with imaginary components indicate solutions that give rise to sustained oscillations. Fundamental concepts describing stability in non-linear systems are delineated in a seminal text by Gray and Scott.30 Many of these ideas form a basis for the SKM studies of Steuer et al.28 The variable in the equations above is defined as the saturation parameter. In Michaelis-Menten enzyme kinetics, a reaction initially is subject to first-order kinetics at small substrate concentrations. As the enzyme becomes saturated with increasing concentration of substrate, the reaction becomes subject to zero-order kinetics. Thus, ranges from one to zero, where of one indicates the first order kinetics regime while of zero indicates that a reaction is in the zero order kinetics regime. Note that our notation differs slightly from that presented in the original paper by Steuer et al. Specifically, we use R to represent individual reaction fluxes whereas Steuer et al. represented reaction fluxes using the symbol & . In contrast, we use V to represent the flux of individual metabolites. Thus, in our notation V is determined by the sum of reaction fluxes R. Consequently, the definition of ! above (equation 7) contains R while Steuer et al. defined ! in terms of & . This difference in notation is relevant to the definitions of J, μ and Λ above. However, the definitions of and x are otherwise identical. 28

The parametric description of J described above allows the propensity of a network to possess certain dynamical properties to be assessed. Note that to date such SKMs have not allowed conclusions to be made about the specific dynamic behavior of a network. However, they have allowed an examination of the dynamic capabilities of a system. There are several key properties of the dynamics that can be evaluated: i) stability of the fluxes with regard to small perturbations ii) damped oscillation of the fluxes iii) instability of the fluxes with regard to small perturbations iv) oscillatory instability of the fluxes. For a system that is able to assume a steady state, either i) or ii) should apply. These dynamic properties are associated with well-defined values for the eigenvalues of J respectively.7 4


Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For a stable node, all real parts of the eigenvalues are negative with no imaginary () components. By considering the structure of J, one can deduce that the linear combination of (* that diagonalizes J in this case acts to oppose changes to the steady state So. In the case of a stable focus, the real parts of the eigenvalues are negative and eigenvalues determining the slowest timescale have complex conjugate imaginary components corresponding to stable but oscillatory dynamics. In contrast to stable systems, a saddle point contains a maximum eigenvalue that is real and positive. Thus, any small perturbation is exponentially amplified. In () this case one can deduce that the linear combination of (* that diagonalizes J acts to amplify changes to So. Finally, for an unstable focus the maximum eigenvalue is positive but contains complex conjugate imaginary components so that the dynamic behavior is oscillatory. Hopf or saddle-type bifurcations are present in between each of the regimes ii through iv.7 In addition, allosteric regulation can be readily incorporated via modifying the matrix + . This can be done by increasing the value denoting the role of a metabolite in a reaction that it allosterically activates to be greater than 1 or decreasing the value to be less than zero to denote allosteric inhibition. 2. Methods 2.1. Analysis of Contributions to the Jacobian We have developed a novel approach based on Principal Component Analysis (PCA) that allows one to identify how chemical components contribute to the stability of reaction networks described via SKM (see below). This is new information that was not available in the original applications of SKM.7, 28-29 The method involves performing PCA of the Jacobian matrices obtained from SKM. Our objective in doing so is to identify elements of the Jacobian that play a dominant role in determining the eigenvalues of the matrix and thus the dynamic properties of the network. In particular, Jacobian elements that consistently make large contributions to large, negative eigenvalues will be most likely to play a significant role in determining stable dynamics of the reaction network. The central basis of this effort rests on the following matrix equality associated with diagonalizing the Jacobian: XLT*J*XR= XLT*XR*D (9) In this equation is the transpose of the left eigenvector, J is the original Jacobian matrix, XR is the right eigenvector, and D is the diagonal matrix of eigenvalues. The left side of the equality in which the Jacobian matrix is multiplied by each eigenvector is of particular interest. We compute the values obtained by multiplying a given Jacobian element by the two eigenvectors and keep track of the indices of Jacobian elements contributing to each eigenvalue. The products of the Jacobian elements and eigenvectors are then normalized by dividing each by the component from the right side of the equality above that corresponds to the appropriate eigenvalue. The relative magnitudes of these contributions are then recorded. To illustrate this procedure, below is a schematic representation of the collection of Jacobian elements that contribute to a given eigenvalue DA. XLT

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29

αJ11 + β J12 + γ J13 + ... ω Jmm = κADA

(10) The normalized contribution of Jacobian index J11 to the solution corresponding to eigenvalue DA is thus given by , ⁄κ. /0 , where κ0 is a constant. The manipulations described in equations 9 and 10 are explicitly delineated in the supporting material for a simplified Jacobian matrix. In carrying out this analysis, it is essential to understand the significance of applying PCA to the Jacobian matrix. PCA is a widely used and well understood technique that has been employed in numerous contexts to identify the components of a data set that contribute the most variability to that data.31 In general applications, PCA is used to reduce the dimensions of multidimensional data sets. When there are many more variables compared to the number of samples, the individual samples can often be characterized by only a few variables and correlations between the variables can be identified within the data set via PCA.31 When applied to the covariance matrix generated from molecular dynamics simulations, PCA is used to extract the principal components of motion (also known as the essential dynamics) for those simulations.32 This reveals correlations between individual groups involved in collective motions as well the relative amplitude of these motions. In the color technology industry, spectral properties describing a large data set can be reduced using spectra obtained via PCA.33 These applications are similar in principle to our method since we are trying to identify the essential elements of the Jacobian that play a dominant role in determining the eigenvalues and thus stability in the system. In a similar manner to other applications of PCA, it is reasonable to expect that applying PCA to the Jacobian will allow us identify those components of the Jacobian that make the largest contributions to defining the eigenvalues. To investigate this possibility, we first applied PCA to symmetric Jacobians for model reaction systems. PCA is readily interpretable when applied to symmetric matrices such as a covariance matrix. In such cases, the left eigenvector is the transpose of the right eigenvector in equation 9 and the eigenvectors are orthonormal. Consequently, we reasoned that applying PCA to symmetric Jacobians describing reaction networks should also be readily interpretable. This analysis should reveal the elements of the Jacobian that play the largest role in defining solutions in a similar manner to the way that PCA of covariance matrices reveals the components of a system that are responsible for the largest correlated variability. We did find this to be the case for several simplified reaction schemes representing key steps in glycolysis (see supporting material). One would anticipate that components of the system with a large role in defining stable solutions will demonstrate the most recalcitrance to being perturbed when the system is disrupted, such as it is when the enzyme associated with that matrix element is knocked down. We further demonstrated that PCA of asymmetric Jacobians results in a similar interpretation as that of PCA for symmetric Jacobians (see supporting material). This result indicates the general applicability of performing PCA for decomposing the contributions of arbitrary Jacobians to eigenvalues obtained via equation 9. This issue is important because Jacobians in general are asymmetric. In fact, symmetric Jacobian tend to be the exception rather than the rule for realistic reaction systems. Given these findings, it is straightforward to extend what is known about PCA of other systems to analysis of Jacobians by analogy.

6


Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2. Application to Glycolysis and Pentose Phosphate Pathway Consequently we applied our analysis to our primary system of interest: the well-studied metabolic pathway of glycolysis in yeast also modeled by Steuer et al.28 A schematic depiction of this pathway is shown below. The extensive information available for glycolysis allows us to validate our findings using previous computational and experimental studies.4-6, 24, 34-36 We also studied the pentose phosphate pathway in yeast individually. The pentose phosphate pathway is the primary source of NADPH in Saccharomyces cerevisiae.6 NADPH is essential for life and takes part in redox reactions used to generate important molecules like DNA and proteins.6 Finally, we combined glycolysis with the pentose phosphate pathway to better understand how connections between the two pathways impact their dynamic properties. We did this by using all the reactions from both pathways to make up a new system of 20 reactions. These reactions are shown in figure 3. The first 13 reactions make up glycolysis while the remaining 7 are from the pentose phosphate pathway. To combine the pathways, reactions and metabolites from both the pathways were used to generate a single stoichiometric matrix. Thus, the metabolites that are shared by both reaction networks allow for communication between them. These shared metabolites include glucose-6-phosphate, fructose-6-phosphate and glyceraldehyde-3-phosphate (colored purple, light blue and orange respectively in Figure 3 below). The black block arrows in Figure 3 represent each enzymatic step in the pathway and are numbered sequentially for later reference (see Table 1). The circles containing A, N, and NP represent points where the cofactor pairs ATP/ADP, NADH/NAD+, and NADPH/NADP+ are involved in a reaction. Dashed lines represent instances where the products of one reaction are not the reactants of the following reaction. In such cases the next reactants in the pathway are drawn below the line. Finally, the triangles represent points of regulation. In the second step of glycolysis, hexokinase is inhibited by its product glucose-6-phosphate.37 In the fourth step, phosphofructokinase (PFK) is allosterically inhibited by ATP.34 In the tenth step in the pathway pyruvate kinase (PYK) is allosterically activated by fructose bisphosphate (FBP) and allosterically inhibited by ATP.34 We focused our study of regulation on inhibition of PFK by ATP, inhibition of PYK by ATP, and activation of PYK by FBP via changes in the saturation parameters () as will be described below. These sites have been suggested to be the primary points of regulation in glycolysis.34 Table 1. The Metabolites and their Abbreviations used in the text are shown here. Metabolite Abbreviation Glycolysis Glucose Glucose-6-phosphate G-6-P Fructose-6-phosphate F-6-P Fructose-1,6-bisphosphate FBP Dihydroxyacetone phosphate DHAP Glyceraldehyde-3-phosphate GAP Triose Phosphates DHAP and GAP 1,3-Bisphosphoglycerate BPG 3-phosphoglycerate 3-PG 2-phosphoglycerate 2-PG Phosphoenolpyruvate PEP

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Pyruvate Ethanol Pentose Phosphate Pathway 6-phospho-gluconolactone Ribulose-5-phosphate Ribose-5-phosphate Xyulose-5-phosphate Sedoheptulose-7-phosphate Erythrose-4-phosphate Cofactors Adenosine triphosphate Adenosine diphosphate Nicotinamide adenine dinucleotide Nicotinamide adenine dinucleotide phosphate

PYR EtOH 6-PG Ribulose-5-P Ribose-5-P Xyulose-5-P Sedoheptulose-7-P Erythrose-4-P ATP ADP NADH/ NAD+ NADPH/ NADP+

Table 2: List of Reactions and Corresponding Reaction Numbers Shown in Figure 3. Index Reaction Abbreviation 1 Glucose Intake into pathway 2 Hexokinase 3 Phosphoglucoisomerase 4 Phosphofructokinase PFK 5 Fructose-bisphosphate aldolase 6 Triose phosphate isomerase 7 Glyceraldehyde dehydrogenase 8 Phosphoglycerate kinase 9 Phosphoglycerate mutase 10 Enolase 11 Pyruvate Kinase PYK 12 Ethanol production 13 ATPases 14 Glucose-6-phosphate dehydrogenase 15 Gluconolactonase/6-phsophogluconate dehydrogenase 16 Phosphopentose isomerase 17 Phosphopentose epimerase 18 Transketolase 19 Transaldolase 20 Transketolase Only enzyme abbreviations utilized in this document are listed in the third column

8


Page 8 of 29

Page 9 of 29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. Reaction scheme for the interconnected glycolysis and pentose phosphate reaction network utilized in SKM models. The arrows corresponding to the enzymes catalyzing each reaction are numbered and also listed in Table 2. Models 1-6 used reactions 1-13 for glycolysis. Model 7 used reactions 14-20 for the pentose phosphate pathway. Models 8 and 9 employed all reactions 1-20 in combination. The metabolites and their abbreviations are listed in Table 1 and are shown in between the reaction arrows. The circled A represents a role for the cofactors ATP/ADP while the circled N represents NADH/NAD+ and the circled NP represents NADPH/NADP+. Metabolites common to both pathways (G-6-P, GAP and F-6-P) are shown in between the pathways in the purple, orange, and blue rounded rectangles respectively and allow for communication between the two pathways. Corresponding colored arrows show the connections between the two pathways and the common pools of these metabolites. Dashed lines in between reactions indicate that the product of the previous reaction is not the reactant of the following reaction, with the line dividing the preceding product from the new reactant. Points of regulation in glycolysis are shown by the blue and red triangles, which represent inhibition and activation respectively.

10


Page 10 of 29

Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3. Simplified Model of Glycolysis As described in the Supporting Material, our studies of simple reaction schemes representing key steps in glycolysis verified that applying PCA predicts similar contributors to stability when analyzing both symmetric and asymmetric Jacobians. We then applied the analysis to a more realistic reaction network describing glycolysis in greater detail. This network was originally described by Steuer et al.28 and is displayed in Figure 3 below. We also generated a reversible version of this network and verified that compatible observations are obtained for both cases (not shown). The Jacobian matrix for this reaction scheme without negative feedback as well as the concentration values used in the model are given in supporting information.

Figure 3. Steuer et al. Simplified Model of Glycolysis. This figure shows the reaction scheme from the model of glycolysis shown in Figure 2 of Steuer et al.28 Negative feedback on Phosphofructokinase in Reaction 1 by ATP is shown. Reaction 1 of this figure accounts for enzymes Hexokinase, Phosphoglucoisomerase and Phosphofructokinase. Reaction 2 represents Fructose bisphosphate aldolase. Reaction 3 incorporates Triose phosphate isomerase and Glyceraldehyde dehydrogenase. Reaction 4 accounts for Phosphoglycerate kinase, Phosphoglycerate mutase, Enolase and Pyruvate kinase. Reaction 5 is Ethanol production. Reaction 6 describes the transport of Pyruvate to Citric Acid Cycle (this reaction is not included in our model). Reaction 7 describes the connection to the Pentose Phosphate Pathway. Finally, reaction 8 represents the general action of ATPases. All values were set to 1 and the feedback of reaction 1 set to 0 indicating no inhibition, resulting in an asymmetric Jacobian. The Jacobian matrices obtained were diagonalized as described previously in Methods. The θ position denoting the presence of feedback was kept constant at values of 0, 1 and 2 for different simulations. At each feedback value we completed 1000 trials while randomizing the distributions of the remaining θ entries. We created histograms of the normalized contributions for each of the 36 entries in the Jacobian. A spectrum of the eigenvalues for the three individual sets of trials was also generated that can be compared with the spectrum of eigenvalues presented in the work of Steuer et al. (see supporting material). Sets of 1000 simulations were performed on a local Windows PC with an Intel i5-2500K processor running Windows 10 Pro and using MATLAB 2016a. In these simulations construction of the Jacobian typically required less than 4 hours. Identifying contributors to the eigenvalues was significantly faster, requiring less than an hour to analyze 1000 simulations. Analysis of the reaction schemes described above reveals that some Jacobian elements have larger contributions to the eigenvalues. Moreover, the magnitude of these contributions is 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

consistent with intuitive expectations regarding the relative importance of Jacobian elements in determining overall reaction flux (see Results). The Jacobian elements with large contributions to the largest eigenvalues impact more of the reactions in the network, while Jacobian elements with small contributions have more of a localized effect on the reaction network. Whether the Jacobian matrix is symmetric or asymmetric was not observed to affect the nature of the contributions, indicating that this method can be applied to both types of matrices. 2.4. Fully Detailed Models of Glycolysis, Pentose Phosphate Pathway, and Combined Glycolysis/Pentose Phosphate Network We constructed a more detailed description of glycolysis than that shown in Figure 3 by explicitly including as individual steps the reactions that had been grouped together by Steuer et al.28 (see Table 3). We also examined the pentose phosphate pathway, explicitly including each of its constituent reactions. The individual glycolysis and pentose phosphate pathways were first studied in isolation. We then combined the glycolysis and pentose phosphate pathway metabolic networks to understand how interconnections between them influence the results of the analysis, resulting in 20 reactions for the combined pathway. Metabolite concentrations were taken from literature, please see supplemental information for a detailed list of data sources.4-6, 28, 34, 36, 38 Simulations were conducted with θ values varying between 0 and 1 for reactions not involved in regulation or affected by cofactor ratios. Values of θ describing regulation of PFK and PYK were set as follows. For PFK inhibition by ATP the θ for ATP in the PFK reaction (reaction 4 in figure 3) was set to -1 for models 2, 5 and 9. For PYK activation of FBP, the θ for FBP in the PYK reaction (reaction 11 in Figure 3) was set to 2 for models 3, 5 and 9. For PYK inhibition by ATP the θ for ATP in the PYK reaction was set to -1 for models 4, 5 and 9. Model 1 describes glycolysis without the impact of regulation. Models 2-4 of glycolysis each had only one point of regulation specified by the θ noted above. Model 5 combines all three types of regulation specified above. The SKM version of the kinetic model by Smallbone et al. as well as the Pentose Phosphate pathway on its own did not incorporate regulation. One combined glycolysis/pentose phosphate pathway model (model 8) did not include regulation while the other combined pathway (model 9) included all three forms of regulation noted above. The different forms of enzyme regulation replicate known regulatory mechanisms present in glycolysis and allow us to determine whether differential regulation impacts the fundamental dynamic properties of the model. The systems studied are listed in Table 4: Table 3: Description of Simulation Systems Studied Model System No Regulation PFK Model Inhibition by Number ATP 1 Glycolysis + 2 Glycolysis + 3 Glycolysis 4 Glycolysis 5 Glycolysis + 6 Smallbone + Kinetic Model SKM version 7 Pentose+

12

PYK Activation by FBP


PYK Inhibition by ATP

+ +

+ +

Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Phosphate 8 Combined + 9 Combined + + : indicates the presence of a given set of simulation conditions

+

+

Analysis of the 9 systems was done by constructing histograms of the Jacobian contributions and eigenvalue spectra. We also compared the results from our SKM model of glycolysis to the kinetic model of yeast glycolysis described by Smallbone et al. (i.e. Kinetic model 184) to test the SKM predictions. We accessed the model of Smallbone et al. via Mathematica. We identified metabolites predicted to induce stability or instability in our SKM models and used the kinetic model of Smallbone et al. to vary the steady state concentrations of these metabolites by decreasing and increasing their concentrations relative to their original steady state values. This was done in order to determine how resistant these metabolite concentrations are to perturbations in the kinetic model. To further compare our SKM models of glycolysis with the kinetic model of Smallbone et al., we constructed an SKM (model 6 in Table 4) using fluxes and initial concentrations identical to those reported by these authors (http://jjj.mib.ac.uk/models/smallbone18/simulate/) and computed the resulting Jacobian contributions and eigenvalue spectra. We performed 105 simulations for each of the 9 systems using MATLAB 2016a at the UMBC High Performance Computing Facility (HPCF).

3. Results 3.1. Simple Model of Glycolysis We assessed Jacobian contributions for the glycolysis model originally studied by Steuer et al. that is depicted in Figure 3. In doing the analysis we reproduced the same dynamic properties as observed by these authors. Our analysis revealed that there are specific Jacobian elements with consistently larger contributions to the eigenvalues than others, just as observed for the simplified models described in the supporting material (see section Analysis of Contributions to the Jacobian). Thus, these components of the reaction network are likely to be the primary determinants of network stability. As noted by Steuer et al., this system does not generate any positive eigenvalues. For the largest negative eigenvalue, the major contributor () was (*234 . The flux 567 is determined by the enzymes glyceraldehyde dehydrogenase, triose 234

phosphate isomerase, phosphoglycerate kinase, phosphoglycerate mutase, enolase and pyruvate kinase (see Figure 3 and Table 3). In the process, ATP is produced at two distinct points. The fact that the largest eigenvalue has greatest contributions from the two ATP producing enzymes seems significant. An important role of glycolysis is the production of energy from glucose and the enzymes listed above are important contributors to this process. The second largest () eigenvalue had the greatest contribution from (*83 .The flux 96 is governed by the enzyme 83

glyceraldehyde dehydrogenase. The third largest eigenvalue had the greatest contribution from ():;0?@ involves generation of BPG (see above) as well as the anaerobic (* :;

Using Structural Kinetic Modeling To Identify Key Determinants of

Recommend Documents