As a first step toward the elucidation of the systems biology of …

• Abstract
• Introduction
• Experimental Procedures
• Results
• Discussion
• References
• Figures
• Reactions

The Mathematical Model—Here we use kMech (10), a Cellerator (11) language extension that describes a suite of enzyme reaction 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 (ODEs) that are numerically solved by Mathematica^TM (10). To build a model for a metabolic pathway, users need only to string together appropriate kMech enzyme mechanism models and provide the physical and kinetic parameters for each enzyme. The development of approximation methods for estimating unavailable model parameters such as forward and reverse rate constants (k_f and k_r) from kinetic measurements (K_m and k_cat) is described elsewhere (10).

The detailed kMech models for each of the pathway enzymes (Fig. 1), a Mathematica^TM-executable kMech.m file, and a Mathematica^TM notebook file with detailed kMech inputs and the corresponding ODEs, kinetic rate constants, and initial conditions for solving the ODEs (or its Systems Biology Markup Language (SBML) version) are available at the University of California, Irvine, Institute for Genomics and the Bioinformatics web site (www.igb.uci.edu/servers/sb.html). The PDF version of the Mathematica^TM notebook and a list of reported and optimized enzyme kinetic and physical parameters used in simulations are available in the supplemental data in the on-line version of this article. Cellerator, available at www.aig.jpl.nasa.gov/public/mls/cellerator/feedback.html, is free of charge to academic, United States government, and other nonprofit organizations.

Carbon Flow Channeling—Traditional modeling approaches 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 focus on conversion between metabolites (metabolic flux) rather than enzyme mechanisms. Although metabolic flux provides valuable information about biomass conversions (12), it cannot simulate, for example, the pathway-specific regulation patterns that control carbon flow channeling through the three AHAS isozymes of the parallel L-isoleucine and L-valine pathways and the final transamination reactions. This level of mathematical modeling requires a detailed understanding of enzyme kinetic mechanisms and regulatory circuits (Fig. 2) as described below.

-Acetohydroxy Acid Synthase (AHAS) Isozymes—The AHAS isozymes are controllers of carbon flow into either the L-isoleucine or L-valine biosynthetic pathway. The Ping Pong Bi Bi enzyme mechanism of these isozymes describes a specialized two-substrate, two-product (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 indicated as white and shaded ovals, respectively, in Fig. 2.

Carbon flow through these isozymes is controlled by the affinities (K_m) of the enzyme intermediates for their second substrates as shown in Fig. 2 (13). For example, the AHAS II enzyme reactions shown in Fig. 1 are described by the two reaction sets designated Reaction 1 and Reaction 2.

Reaction 1 is for the condensation of two pyruvate (Pyr) molecules for the biosynthesis of the -acetolactate of the L-valine and L-leucine pathways. Reaction 2 is for the condensation of one Pyr molecule and one KB molecule for the biosynthesis of the -acetohydroxybutyrate of the L-isoleucine pathway. As written in the reactions, the initial reaction of Pyr with free AHAS II to form the activated enzyme intermediate is represented twice. Therefore, if these reactions were modeled, two molecules of Pyr would be consumed instead of just one for each molecule of Pyr or KB produced. This redundancy can be resolved by rewriting these reactions as shown in Reaction 3.

In Mathematica^TM, the Union operator is used in conjunction with the kMech PingPong models to solve this redundancy of pathways problem. The user kMech inputs in Mathematica^TM syntax for the channeling of pyruvate through the AHAS II isozyme into L-valine or L-isoleucine are shown in Reaction 4. In this reaction, Pyr and aKB (-acetolactate) are substrates, aAL (-acetolactate), aAHB (-acetohydroxybutyrate), and CO₂ are products, AHASII is free enzyme, AHASIICH₃CO is the modified enzyme intermediate (it is underlined in Reactions 1, 2, 3 to indicate that it is in the intermediate state after reacting with the first substrate), Enz[...] denotes a kMech enzyme model that provides additional capabilities to Cellerator, PingPong indicates that the enzyme model is Ping Pong Bi Bi, variable names with a kf prefix are rate constants of the enzyme-substrate associations, variable names with a kr prefix are rate constants of the enzyme substrate dissociations, and variable names with a kcat prefix are catalytic rate constants for the formation of products.

kMech parses the three non-redundant AHAS II reactions shown above into elementary association-dissociation reactions and produces the output in Cellerator/Mathematica^TM syntax (11) shown in Reaction 5. This output is passed to Cellerator where the differential equations that describe the rate of change for each reactant involved in the AHAS II isozyme reaction are generated in Mathematica^TM syntax and presented collectively in Reaction 6. These differential equations and variable definitions are passed to Mathematica^TM, where they are solved by the numeric solver (NDSolve) function, and graphs of enzyme product versus time are generated.

The Union operator also was used for the modeling of the L-valine-inhibited AHAS I and AHAS III isozymes described in the supplemental data in the on-line version of this article. Detailed descriptions of other kMech models used in this simulation are published elsewhere (10).

Reversible Transamination Mechanism—The pyridoxal 5'-phosphate-dependent transaminase B (TB) enzyme catalyzes the final, reversible step of the biosynthetic pathways of all three of the branched chain amino acids (Figs. 1 and 2). The first step of each of these Ping Pong Bi Bi transamination reactions uses glutamate as an amino donor to form a pyridoxamine-bound enzyme intermediate (TBNH₂, shaded oval in Fig. 2) for the transamination of the three different -ketoacids of each pathway. Carbon flow through TB is controlled by the affinities (K_m) of the enzyme intermediates for their second -ketoacid substrates as shown in Fig. 2. The TB enzyme reactions of Fig. 1 are described by the chemical equations shown in Reactions 7, 8, 9. Because the first substrate reaction with glutamate is the same for all three of the branched chain -ketoacid second substrates, the Mathematica^TM Union operator is once again used to eliminate this redundancy as shown in Reaction 10. Because transamination is reversible, kMech models must be entered in both reaction directions for each of the three branched chain amino acid transaminations and, again, the Union operator is used to eliminate the duplicated second substrate reactions (TBNH₂ + aKG TBNH₂·aKG TB + Glu; TBNH₂ is underlined to indicate that the TB enzyme is in the intermediate state after reacting with the first substrate) of each transamination (long gray arrows in Fig. 2) as shown in Reaction 11.

These reactions are parsed by kMech into elementary association-dissociation reactions and passed on to Cellerator, where they are processed as described above. The same method was used for modeling transaminase C, a reversible Ping Pong Bi Bi mechanism enzyme that uses alanine as the amino donor for the transamination of L-valine (Fig. 2).

Allosteric Regulation—Threonine deaminase is an allosteric enzyme whose kinetic behavior can be described by the concerted allosteric transition model of Monod, Wyman, Changeux known as the MWC model (7, 8). According to the MWC model, TDA can exist in an active state (R) or an inactive state (T) (8, 9). The fraction of enzyme in the R or T state is determined by the concentrations and relative affinities of the substrate (L-threonine), the inhibitor (L-isoleucine), and the activator (L-valine) for each state. This model is described by two equations, designated Equations 1 and 2,

(Eq. 1)

(Eq. 2)

in which S, I, and A are substrate, inhibitor, and activator concentrations, respectively, K_m, K_i, and K_a are their respective dissociation constants, n is the number of substrate and effector ligand binding sites, c is the ratio of the affinity of the substrate for the catalytically active R state and the inhibited T state, L₀ is the equilibrium constant (allosteric constant) for the R and T states in the absence of ligands, v_o is the initial reaction velocity, and V_max is the maximal reaction velocity.

Equation 1 describes the fraction of the enzyme in the catalytically active state (R) as a function of substrate and effector concentrations. Equation 2 describes the fractional saturation (Y_f = v_o/V_max) of the enzyme occupied by substrate as a function of substrate and effector concentrations (7).

We have recently described implementation of the MWC model in Cellerator (10). Experimental values of the kinetic parameters and ligand concentrations listed above are most often available in the literature. However, values of c and L₀ are often not available. These values can be calculated by fitting substrate saturation curves in the presence and absence of several inhibitor concentrations (10, 14, 15).

Approximation of Intracellular Enzyme Concentrations—With few exceptions, intracellular enzyme concentrations are not available. However, with careful experimental documentation, these concentrations can be approximated from the yields and specific activities of purified enzymes. For example, calculations based on purification tables in the literature suggest that the intracellular concentration of TDA is 4 µM (16). Furthermore, recent experiments have shown a positive correlation between mRNA levels measured with DNA microarrays and protein abundance in both E. coli (17) and yeast cells (18, 19). Thus, the intracellular levels of the remaining enzymes of the branched chain amino acid biosynthetic pathway can be inferred from the calculated intracellular level of TDA and the relative mRNA levels of the other branched chain amino acid biosynthetic enzymes using DNA microarray data (20). The data in supplemental Table I in the on-line version of this article demonstrate that this is a reasonable method. Indeed, simulations using intracellular enzyme concentrations inferred in this manner produce experimentally observed steady-state pathway intermediate and end product levels (21, 22), usually within 2-fold to one-half adjustments of these inferred values.

Optimization of Model Parameters—A list of reported enzyme kinetic and physical parameters needed to solve the differential equations for the simulations reported here, and their literature sources are available in supplemental Table I in the on-line version of this article. The optimized values to simulate known steady-state intracellular levels of pathway substrates, intermediates, and end products are also listed for comparison. In brief, for each enzyme there are at least three parameters needed, namely the total enzyme concentration (E_T), the K_m for each substrate, and the k_cat for each enzyme reaction. For enzymes with additional regulatory mechanism, extra parameters such as the K_i for each inhibitor and the K_a for each activator also are required.

In the initial simulation, E_T values were inferred from microarray data as described above, and the K_m and k_cat values were obtained from in vitro enzyme kinetic data of purified enzymes with the exception of transaminase C and -isopropylmalate isomerase, where empirical values were used because of a lack of experimental data. These values were manually adjusted to match the published in vivo steady-state levels of intermediate and end product metabolites (21, 22). Interestingly, the inferred E_T and in vitro K_m values work quit well, because the adjustments are usually within 2-fold to one-half of the initial values. However, because many variables can influence in vitro measurements, including the relative activities of purified enzymes, larger corrections were sometimes necessary for the estimation of k_cat values (5 of 9 enzymes). Once the mathematical model was optimized with the parameters reported in supplemental Table I in the on-line version of this article, it was used without further adjustment for the simulations of the metabolic and genetic perturbations reported below.