Bioinformatics Advance Access originally published online on October 27, 2004
Bioinformatics 2005 21(6):774-780; doi:10.1093/bioinformatics/bti068
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
An enzyme mechanism language for the mathematical modeling of metabolic pathways
1Department of Microbiology and Molecular Genetics, College of Medicine, University of California Irvine, CA 92697, USA
2School of Information and Computer Science, University of California Irvine, CA 92697, USA
3Institute for Genomics and Bioinformatics, University of California Irvine, CA 92697, USA
4Jet Propulsion Laboratory, California Institute of Technology 4800 Oak Grove Drive, Pasadena, CA 91109, USA
*To whom correspondence should be addressed.
 |
Abstract |
|---|
 TOP Abstract INTRODUCTIONSYSTEMS AND METHODSIMPLEMENTATIONDISCUSSIONREFERENCES |
|---|
Motivation: As a first step toward the elucidation of the systems biology of complex biological systems, it was our goal to mathematically model common enzyme catalytic and regulatory mechanisms that repeatedly appear in biological processes such as signal transduction and metabolic pathways.
Results: We describe kMech, a Cellerator language extension that describes a suite of enzyme mechanisms. Each enzyme mechanism is parsed by kMech into a set of fundamental association–dissociation reactions that are translated by Cellerator into ordinary differential equations that are numerically solved by MathematicaTM. In addition, we present methods that use commonly available kinetic measurements to estimate rate constants required to solve these differential equations.
Availability: A MathematicaTM executable kMech.m file is available at the University of California, Irvine, Institute for Genomics and Bioinformatics website, http://www.igb.uci.edu/servers/sb.html. Cellerator, free of charge to academic, US government, and other non-profit organizations, can be obtained at the Cellerator website, (http://www-aig.jpl.nasa.gov/public/mls/cellerator/feedback.html).
Contact: Biology correspondence should be addressed to gwhatfie{at}uci.edu. Computation correspondence should be addressed to emj{at}uci.edu.
Supplementary information: http://www.igb.uci.edu/servers/sb.html.
|
INTRODUCTION |
|---|
|
TOPAbstractINTRODUCTIONSYSTEMS AND METHODSIMPLEMENTATIONDISCUSSIONREFERENCES |
|---|
Systems biology may be broadly defined as the integration of diverse data into useful biological models that allow scientists to easily observe complex cellular behaviors and to predict the outcomes of metabolic and genetic perturbations. In this report we present kMech, a Cellerator (Shapiro et al., 2003) language extension that describes a suite of enzyme mechanisms suitable for the mathematical modeling of enzyme-related pathways. Cellerator is a tool for generating reaction network models of cellular processes. It contains reaction mechanisms based on mass action kinetics for elementary catalytic and non-catalytic reactions, along with other regulatory relationships (Shapiro et al., 2003). kMech supplements Cellerator with a suite of catalytic reaction mechanisms involving multiple substrates and products and regulatory mechanisms. In kMech, each mechanism has been codified to generate a set of elementary reactions that can be translated by Cellerator into ordinary differential equations (ODEs) solvable by MathematicaTM, a widely used commercial computer algebra system that integrates numeric and symbolic computational engines with a graphical output and a programming language. Alternatively, Cellerator can generate ODEs in systems biology markup language (SBML) for other simulators (Hucka et al., 2003).
Traditional enzyme modeling approaches, for example, GEPASI (Mendes, 1997) or JARNAC (Sauro et al., 2003), use the Michaelis–Menten kinetic equation for one substrate/ one product reactions and the King–Altman method to derive equations for more complex multiple reactant reactions. These types of equations are called steady-state velocity equations since the derivatives of the concentration of each reactant in the model over time are set to zero to simplify a set of non-linear differential equations to linear algebric equations (Falk et al., 1998). Therefore, the kinetic model based on this approach has embedded in it the steady-state hypothesis. In contrast, the model generated by kMech/Cellerator consists of non-simplified, non-linear, differential equations that describe the rates of change of each reactant in the model over time. For the solution of these differential equations, users can input the rate constants (kf, kr) if they have been experimentally measured. Alternatively, these rate constants can be approximated from easily measured kinetic constants (Km, kcat) using the Lambda and Omega approximation methods presented in this paper. Moreover, if the assumption of equilibrium between substrate and enzyme is not favored, values of Lambda and Omega can be adjusted to fit the experimental data.
kMech is built in a modular manner such that complex enzyme mechanisms may be constructed from simpler ones. This design is sufficiently versatile to accommodate enzyme mechanisms such as single and multiple substrate (e.g. Ping Pong and Bi Bi) enzyme kinetic reactions, and feedback inhibition (e.g. allosteric, competitive and non-competitive) mechanisms.
|
SYSTEMS AND METHODS |
|---|
|
TOPAbstractINTRODUCTIONSYSTEMS AND METHODSIMPLEMENTATIONDISCUSSIONREFERENCES |
|---|
The development of kMech models for the simulation of enzyme reaction mechanisms
Simple catalytic model
The simple catalytic model describes a single-substrate, single-product, non-reversible catalytic mechanism that involves no cofactors or enzyme intermediates, is implemented in Cellerator.
 |
is the Cellerator notation for an enzyme-catalyzed reaction.) where S is the substrate, P is the product, En is the free enzyme, EnS is the enzyme–substrate complex, kf is the rate constant for enzyme–substrate association, kr is the rate constant for enzyme–substrate dissociation and kcat is the rate constant for the formation and dissociation of product P.
Cellerator converts this reaction input into its basic association–dissociation reactions (Note: Mathematica syntax for all association–dissociation reactions and differential equations is shown in Supplementary Materials.):
 |
 |
For a reversible catalytic reaction, Cellerator employs two simple catalytic reactions, one for the forward reaction and one for the reverse reaction. For example, the Cellerator input for a reversible simple catalytic reaction is
 |
Bi Bi model
The Bi Bi model describes a generalized two-substrate two-product (Bi Bi) reactions model where product formation occurs only after the formation of an enzyme–two substrate complex.
 |
kMech converts this input into association–dissociation reactions defined in Cellerator as follows:
 |
The first reaction represents the formation of the EnS1S2. The second reaction represents the release of products and free enzyme from the complex. Given these reactions, the Cellerator ‘interpret’ function generates the following differential equations for each reactant:
 |
Ping Pong Bi Bi model
The Ping Pong Bi Bi model describes a specialized Bi Bi mechanism in which the binding of substrates and release of products is ordered. It is a Ping Pong mechanism because the enzyme shuttles between a free and a substrate-modified intermediate state:
 |
 |
The first reaction represents the formation of the enzyme–substrate complex with S1 (EnS1), the release of product P1 from the complex and the formation of the enzyme intermediate Enx. The second reaction represents the formation of the intermediate enzyme–substrate complex with S2(EnxS2) and the release of product P2 and free enzyme, En.
Given these reactions, Cellerator generates the following differential equations for each reactant:
 |
Extensions of enzyme reaction models with regulatory circuits
The enzyme that catalyzes the first step in a biological pathway is usually feedback inhibited by its end product (Umbarger, 1992). Non-competitive, competitive, and allosteric inhibition mechanisms are implemented here.
Non-competitive inhibition (NCI) indicates an inhibitor can bind all enzyme states including free enzyme, intermediate enzyme and enzyme–substrate complex (Segel, 1993).
 |
Competitive inhibition indicates an inhibitor can only bind to the free enzyme state (Segel, 1993).
 |
Allosteric regulation
Allosteric regulation is simulated by the two-state, concerted transition, allosteric model of Monod, Wyman and Changeux (MWC) where the enzyme can exist in either an active state, R, or an inactive state, T (Monod et al., 1965). This model is described by two equations:
 |
 |
 |
In order to simulate the MWC model, several parameters such as substrate (S), activator (A) and inhibitor (I) concentrations and their respective dissociation constants Km, Ka and Ki must be known. With these values, usually available in the literature, the values of 
, ß and 
are calculated. The value of n, the number of substrate and effector ligand binding sites, is also usually available (e.g. n = 4 for Escherichia coli threonine deaminase); however, with few exceptions the values of c, the ratio of the affinity of the substrate for the catalytically active R state and the inhibited T state, and L0, the equilibrium constant (allosteric constant) for the R and T states in the absence of ligands, are unavailable. Nevertheless, these values can be readily derived by fitting generally available substrate saturation curves generated in the presence of several inhibitor concentrations (Hatfield, 1970) as described in the Discussion section.
Lambda approximation method for enzyme rate constants
In order to solve the Cellerator-generated ODEs, values for forward and reverse kinetic rate constants kf, kr are required that often are not available. To deal with this problem, we developed the Lambda (
) approximation method. Briefly, this method calculates these constants from the easily determined, and usually available, enzyme parameters, Km (Michaelis–Menten constant) and kcat (catalytic constant or enzyme turnover number). This approach results in simple mathematical relations that can be used to estimate rate constants of differential equations for enzyme reactions. For example, the single substrate–single product ‘catalytic’ enzyme reaction considered earlier can be represented as:
 | (1) |
 |
 | (2) |
 | (3) |
In the special condition (denoted by *) where [A]* = Km, 50% of En is saturated with A at equilibrium. In other words, under this condition, the amounts of free and substrate-bound enzyme are equal,
 | (4) |
 | (5) |
as a special case of Q
 | (6) |
 | (7) |
is a large time-invariant constant number, and the larger is, the faster En, A and EnA reach equilibrium. By rearrangement of (6), we can approximate kf as,
 | (8) |
 | (9) |
approximation is that two unknown parameters, kf and kr, are reduced to a single unknown parameter .
By substitution of Equations (3) and (9) into the differential equation describing the state of EnA, we obtain,
 | (10) |
is large and close to Q, that is /Q 
1, then Equation (10) (EnA reaches steady state). In other words, when is large, there is a very fast approach to steady state, which is the same as the Michaelis–Menten pseudo steady-state assumption (Murray, 1993). Under these conditions, when is large, kf/kf approaches 1/Km [Equation (4)] which is observable. To determine the usable range of values of that satisfy this criterion, we varied from 10 to 1 000 000 in simulations with no significant changes in the steady levels of enzyme intermediates and end products. Therefore, the value of can be reasonably set to a value of 100. Otherwise, the values of rate constants directly measured or empirically adjusted to fit experimental data can be used in the model. Km values can be obtained from the literature, and kcat (µmol/min/µmol) values can be calculated from specific activities (µmol/min/mg) and molecular weights of purified enzymes. However, because of uncertainties about the percentages of purified enzymes that are active, kcat values often require adjustments to fit experimental data.
The same approximation methods, (8) and (9) are also used for multiple substrate and product enzymes that bind and release reactants in a Ping Pong fashion. In the general case of two-substrate, two-product enzymes that bind and release reactants in a Bi Bi fashion, kf and kr are approximated according to Equations (11) and (12).
 | (11) |
 | (12) |
is the Km for substrate A, and 
is the Km for substrate B.
Omega approximation method for enzyme inhibition rate constants
In addition to needing methods to approximate forward and reverse reaction rates of enzyme substrate binding reactions, kf, kr, we need to be able to estimate the forward and reverse rate constants for enzyme–inhibitor binding reactions, kfi, kri if their values are not known experimentally. To accomplish this we define a value 
that approximates the rate that an enzyme binds to its substrate relative to its inhibitor. Similar to the derivation of , the derivation of is based on a rapid equilibrium assumption. The kfi is approximated by
 | (13) |
 | (14) |
can be an arbitrary number if free enzyme, enzyme–inhibitor complex, and enzyme–substrate complex are at equilibrium. Consequently, the value of may be set to 1. Otherwise, the values of rate constants directly measured or empirically adjusted to fit experimental data can be used in the model. For the purpose of demonstration, a simplified metabolic pathway (Fig. 1) consisting of various enzyme catalytic and regulatory mechanisms described above was constructed using kMech.
|
|
IMPLEMENTATION |
|---|
|
TOPAbstractINTRODUCTIONSYSTEMS AND METHODSIMPLEMENTATIONDISCUSSIONREFERENCES |
|---|
kMech was written in the MathematicaTM language. It converts complex enzyme mechanisms into elementary Cellerator association–dissociation reactions. It was saved in a package format that can be read into MathematicaTM with Cellerator. kMech-modeled simulations can be executed by MathematicaTM software installed in a Microsoft Windows, MacOS or Linux operating system.
|
DISCUSSION |
|---|
|
TOPAbstractINTRODUCTIONSYSTEMS AND METHODSIMPLEMENTATIONDISCUSSIONREFERENCES |
|---|
As interest in systems biology grows, a major challenge is to develop user-friendly software tools for biologists to model biochemical pathways without knowledge of the underlying mathematics. The kMech/Cellerator/Mathematica package described here addresses this need. To mathematically model a biochemical pathway, users need only to invoke kMech models for the enzyme mechanisms of a pathway without entering the differential equations that describe these mechanisms. Because of this simple user input and the integration of kMech, Cellerator and MathematicaTM, human errors are greatly reduced.
MathematicaTM provides an excellent development environment that integrates symbolic user input (kMech) seamlessly with Cellerator, which converts these inputs into association–dissociation reactions and differential equations that are used by MathematicaTM to simulate the models and to supply graphic outputs. For demonstration purpose, a metabolic pathway that consists of five enzymes, ten metabolites, and one enzyme cofactor was constructed using kMech (Fig. 1). The mathematical model for this metabolic system consists of 27 ODEs, with 16 association and dissociation rate constants and 9 catalytic rate constants. The enzymes of these interacting pathways employ five distinct enzyme mechanisms (Enz1 is allosteric; Enz2 is Ping Pong Bi Bi; Enz3 is Bi Bi; Enz4 is simple catalytic and Enz5 is reversible Ping Pong Bi Bi). Enz1, feedback regulated by the pathway end product, was simulated with the MWC model employing physical parameters based on E.coli L-threonine deaminase (Calhoun et al., 1973). With the Michaelis–Menten pseudo steady-state assumption that substrate, enzyme and enzyme–substrate complex are under rapid steady-state equilibrium, the fractional saturation (Yf) values of the enzyme occupied by substrate are equivalent to measured v/Vmax values (Segel, 1993). These values of v/Vmax (Yf) are available from substrate (threonine) saturation data obtained in the presence and absence of inhibitor (isoleucine) reported in the literature (black dots) (Calhoun et al., 1973) (Fig. 2A). The error function (Err) in Figure 2B is the sum of squared differences of experimental substrate saturation data (black dots) that fit theoretical (Yf) curves with values of c from 0 to 1, and L0 from 0 to 40 as discussed in the Systems and methods section. The pair of c and L0 values that produce the best-fit curves, determined with the non-linear programming Mathematica function, FindMinimum, are c = 0.013, L0 = 1.05. The time-dependent approach to steady state for the six pathway intermediates and end products are shown in Figure 1. At t = 0, an initial burst of intermediate synthesis is observed that is inhibited as end product accumulates and feedback inhibits the first enzyme of the pathway. While the steady-state levels of the substrates for Enz3 and Enz4 (M5 and M6) are maintained at a low level, the level of the precursor of the end product, M7, reaches a high steady-state level comparable with that of the end product because of the reversible nature of this reaction catalyzed by Enz5. This pathway simulates the enzyme mechanisms of L-isoleucine biosynthesis in E.coli K12 (Umbarger, 1996). In this case, the last reversible reaction of the pathway is the transamination between L-isoleucine and L-glutamic acid. The detailed kMech inputs, corresponding ODEs, kinetic rate constants and initial conditions for solving the ODEs are presented in Supplementary Figure 3 online.
|
While deterministic, continuous, modeling of the type described here can provide valuable information such as predicted steady-state levels of metabolic substrates, intermediates and end products, and predict the outcomes of biochemical and genetic perturbations that require detailed enzyme kinetic and regulatory mechanisms, a major challenge of this approach is the measurement or estimation of rate constants required for the solution of rate equations. It nevertheless has a long history. For example, as early as 1943, Chance described methods to measure rate constants directly, and reports describing computer-based methods to fit rate constants to match experimental data date back to the 1960s (Barshop et al., 1983; Garfinkel et al., 1968; Schoeberl et al., 2002; Waas et al., 2001; Zimmerle and Frieden, 1989). However, since experimental results required for either the direct measurement or the indirect fitting of all required rate constants are not available for every enzyme, these methods have not been well-suited for large-scale simulations of metabolic networks. To address this problem we have described general methods for the approximation of some of the missing rate constants. These Lambda and Omega approximations, described in the systems and methods section, assume rapid equilibrium conditions between free enzyme and enzyme–substrate complexes, and between free enzyme and enzyme–inhibitor complexes, that allow us to estimate reaction rate constants from available kinetic data (Km, kcat). This facilitates the modeling of systems in which only these types of kinetic data are available.
Clearly, an ultimate goal of systems biology is the complete simulation of cellular systems (Tomita et al., 1999). However, given the uncertainty of rate constant measurements needed for solving differential equations, some investigators have turned their attention to more abstract methods. For example, Reed et al., (2003) have used metabolic flux balance analysis (FBA) methods to simulate steady-state metabolite flux through E.coli pathways representing more than 900 enzyme steps. An advantage of FBA is that it can make predictions about large metabolic networks without knowing anything about individual kinetic parameters. It assumes that for each node of the network the entrance and exit rates must equal one another. In this case, a simple set of equations can be generated stating that total flux into a node minus the total flux out of the node must equal zero. Thus, FBA can provide valuable information about biomass conversions (Edwards et al., 2001). However, since FBA does not consider pathway-specific regulation patterns it is less suited for the simulation of biochemical or genetic perturbations that require detailed enzyme kinetic and regulatory mechanisms. The kMech/Cellerator models described here represent the first step of a ‘bottom-up’ approach to model complex biological processes.
|
Acknowledgments |
|---|
We are appreciative of support from the UCI Institute for Genomics and Bioinformatics and helpful advice from Donald F. Senear. This work was supported in part by National Institutes of Health grants GM55073 and GM68903 to G.W.H. C.-R.Y. is the recipient of a Biomedical Informatics Training (BIT) Program postdoctoral fellowship from the National Library of Medicine T15 LM-07443.
Received on February 23, 2004; revised on September 16, 2004; accepted on September 16, 2004
|
REFERENCES |
|---|
|
TOPAbstractINTRODUCTIONSYSTEMS AND METHODSIMPLEMENTATIONDISCUSSIONREFERENCES |
|---|
Barshop, B.A., Wrenn, R.F., Frieden, C. (1983) Analysis of numerical methods for computer simulation of kinetic processes: development of KINSIM—a flexible, portable system. Anal. Biochem., 130, 134–145[CrossRef][ISI][Medline].
Calhoun, D.H., Rimerman, R.A., Hatfield, G.W. (1973) Threonine deaminase from Escherichia coli. I. Purification and properties. J. Biol. Chem., 248, 3511–3516
Chance, B. (1943) The kinetics of the enzyme–substrate compound of peroxidase. J. Biol. Chem., 151, 553–577
Edwards, J.S., Ibarra, R.U., Palsson, B.O. (2001) In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. Nat. Biotechnol., 19, 125–130[CrossRef][ISI][Medline].
Falk, S., Guay, A., Chenu, C., Patil, S.D., Berteloot, A. (1998) Reduction of an eight-state mechanism of cotransport to a six-state model using a new computer program. Biophys. J., 74, 816–830
Garfinkel, D., Frenkel, R.A., Garfinkel, L. (1968) Simulation of the detailed regulation of glycolysis in a heart supernatant preparation. Comput. Biomed. Res., 2, 68–91[CrossRef][ISI][Medline].
Hatfield, G.W. (1970) The regulation of L-threonine deaminase in Bacillus subtilis by repression and endproduct inhibition. , Purdue PhD Thesis Department of Biological Sciences, Purdue University 103–112.
Hucka, M., Finney, A., Sauro, H.M., Bolouri, H., Doyle, J.C., Kitano, H., Arkin, A.P., Bornstein, B.J., Bray, D., Cornish-Bowden, A., et al. (2003) The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics, 19, 524–531
Mendes, P. (1997) Biochemistry by numbers: simulation of biochemical pathways with Gepasi 3. Trends Biochem. Sci., 22, 361–363[CrossRef][ISI][Medline].
Monod, J., Wyman, J., Changeux, J.P. (1965) On the nature of allosteric transitions: a plausible model. J. Mol. Biol., 12, 88–118[ISI][Medline].
Murray, J.D. (1993) Michaelis–Menten theory: detailed analysis and the pseudo-steady state hypothesis. Math. Biol., 19, 116–117.
Reed, J.L., Vo, T.D., Schilling, C.H., Palsson, B.O. (2003) An expanded genome-scale model of Escherichia coli K-12 (iJR904 GSM/GPR). Genome Biol., 4, R54[CrossRef][Medline].
Sauro, H.M., Hucka, M., Finney, A., Wellock, C., Bolouri, H., Doyle, J., Kitano, H. (2003) Next generation simulation tools: the Systems Biology Workbench and BioSPICE integration. OMICS, 7, 355–372[CrossRef][Medline].
Schoeberl, B., Eichler-Jonsson, C., Gilles, E.D., Muller, G. (2002) Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat. Biotechnol., 20, 370–375[CrossRef][ISI][Medline].
Segel, I.H. Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady State Enzyme Systems, (1993) , NY Equation (VII-66) Wiley, pp. 427.
Shapiro, B.E., Levchenko, A., Meyerowitz, E.M., Wold, B.J., Mjolsness, E.D. (2003) Cellerator: extending a computer algebra system to include biochemical arrows for signal transduction simulations. Bioinformatics, 19, 677–678
Tomita, M., Hashimoto, K., Takahashi, K., Shimizu, T.S., Matsuzaki, Y., Miyoshi, F., Saito, K., Tanida, S., Yugi, K., Venter, J.C., Hutchison, C.A., III. (1999) E-CELL: software environment for whole-cell simulation. Bioinformatics, 15, 72–84
Umbarger, H.E. (1992) The origin of a useful concept–feedback inhibition. Protein Sci., 1, 1392–1395[ISI][Medline].
Umbarger, H.E. Biosynthesis of the Branched-Chain Amino Acids Escherichia coli and Salmonella: Cellular and Molecular Biology, (1996) , Washington, DC ASM Press.
Waas, W.F., Lo, H.H., Dalby, K.N. (2001) The kinetic mechanism of the dual phosphorylation of the ATF2 transcription factor by p38 mitogen-activated protein (MAP) kinase alpha. Implications for signal/response profiles of MAP kinase pathways. J. Biol. Chem., 276, , pp. 5676–5684
Zimmerle, C.T. and Frieden, C. (1989) Analysis of progress curves by simulations generated by numerical integration. Biochem. J., 258, 381–387[ISI][Medline].
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  W. S. Hlavacek, J. R. Faeder, M. L. Blinov, R. G. Posner, M. Hucka, and W. Fontana Rules for Modeling Signal-Transduction Systems Sci. Signal., July 18, 2006; 2006(344): re6 - re6. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||