Dimerization and Bifunctionality Confer Robustness to the Isocitrate Dehydrogenase Regulatory System in Escherichia coli*

Background: Robustness of the Escherichia coli isocitrate dehydrogenase (IDH) regulatory system has been experimentally demonstrated. Results: A biochemically and structurally realistic mathematical model of the IDH system was developed. Conclusion: The model identifies features of the system, including homodimerization of IDH and bifunctionality of its regulatory enzyme, that lead to robustness. Significance: Mathematical modeling can connect enzymology and systems biology. An important goal of systems biology is to develop quantitative models that explain how specific molecular features give rise to systems-level properties. Metabolic and regulatory pathways that contain multifunctional proteins are especially interesting to study from this perspective because they have frequently been observed to exhibit robustness: the ability for a system to perform its proper function even as levels of its components change. In this study, we use extensive biochemical data and algebraic modeling to develop and analyze a model that shows how robust behavior arises in the isocitrate dehydrogenase (IDH) regulatory system of Escherichia coli, which was shown in 1985 to experimentally exhibit robustness. E. coli IDH is regulated by reversible phosphorylation catalyzed by the bifunctional isocitrate dehydrogenase kinase/phosphatase (IDHKP), and the level of IDH activity determines whether carbon flux is directed through the glyoxylate bypass (for growth on two-carbon substrates) or the full tricarboxylic acid cycle. Our model, which incorporates recent structural data on IDHKP, identifies several specific biochemical features of the system (including homodimerization of IDH and bifunctionality of IDHKP) that provide a potential explanation for robustness. Using algebraic techniques, we derive an invariant that summarizes the steady-state relationship between the phospho-forms of IDH. We use the invariant in combination with kinetic data on IDHKP to calculate IDH activity at a range of total IDH levels and find that our model predicts robustness. Our work unifies much of the known biochemistry of the IDH regulatory system into a single quantitative framework and highlights the importance of constructing biochemically realistic models in systems biology.

An important goal of systems biology is to develop quantitative models that explain how specific molecular features give rise to systems-level properties. Metabolic and regulatory pathways that contain multifunctional proteins are especially interesting to study from this perspective because they have frequently been observed to exhibit robustness: the ability for a system to perform its proper function even as levels of its components change. In this study, we use extensive biochemical data and algebraic modeling to develop and analyze a model that shows how robust behavior arises in the isocitrate dehydrogenase (IDH) regulatory system of Escherichia coli, which was shown in 1985 to experimentally exhibit robustness. E. coli IDH is regulated by reversible phosphorylation catalyzed by the bifunctional isocitrate dehydrogenase kinase/phosphatase (IDHKP), and the level of IDH activity determines whether carbon flux is directed through the glyoxylate bypass (for growth on two-carbon substrates) or the full tricarboxylic acid cycle. Our model, which incorporates recent structural data on IDHKP, identifies several specific biochemical features of the system (including homodimerization of IDH and bifunctionality of IDHKP) that provide a potential explanation for robustness. Using algebraic techniques, we derive an invariant that summarizes the steady-state relationship between the phospho-forms of IDH. We use the invariant in combination with kinetic data on IDHKP to calculate IDH activity at a range of total IDH levels and find that our model predicts robustness. Our work unifies much of the known biochemistry of the IDH regulatory system into a single quantitative framework and highlights the importance of constructing biochemically realistic models in systems biology.
During growth on a two-carbon substrate, Escherichia coli diverts significant carbon flux to the glyoxylate bypass, an anapleurotic pathway that circumvents two decarboxylation steps in the tricarboxylic acid cycle (1). The flux through the glyoxylate bypass is determined by the level of active isocitrate dehydrogenase (IDH), 2 a key tricarboxylic acid cycle enzyme that is regulated by reversible phosphorylation by the bifunctional isocitrate dehydrogenase kinase/phosphatase (IDHKP). Phosphorylation by IDHKP inactivates IDH, limiting the synthesis of ␣-ketoglutarate and directing flux through the bypass. When E. coli is transferred to a substrate containing more than two carbons, IDHKP dephosphorylates IDH to reactivate the tricarboxylic acid cycle as the primary metabolic pathway. Maintaining precise control of IDH activity is thus essential for the ability of E. coli to switch between growth on two-carbon and other substrates. Phosphorylation of IDH by IDHKP has also served as a classic model system for studying the role of protein posttranslational modification in complex cellular processes, and a number of fundamental biochemical discoveries were made from examination of the IDH regulatory system. IDH was the first substrate identified for a prokaryotic protein kinase, establishing that protein phosphorylation is not restricted to eukaryotes (2). IDHKP was one of the earliest examples of a multifunctional enzyme to be discovered and characterized (3), and the IDH regulatory system provided the original experimental basis for enzymatic zero-order ultrasensitivity (4).
Because of the physiological importance of the glyoxylate bypass and their unusual biochemical characteristics, IDH and IDHKP have been extensively studied using methods from classical biochemistry and molecular biology, as reviewed in Refs. 5-7. In a study of particular interest, in 1985, LaPorte et al. (8) provided experimental evidence that the IDH regulatory system exhibits robustness: the ability of a biochemical network to perform its correct function even as levels of the chemical species in the network fluctuate. During growth on acetate, IDH activity varied minimally even as the total amount of IDH present fluctuated 1.8-fold (natural variation between strains) and more than 20-fold (between natural levels and overproduction from a plasmid) (8).
This type of robust behavior is not unique to the IDH regulatory system. Similar experimental observations have been made for the nitrogen assimilation system and the EnvZ/OmpR two-component osmoregulatory system of E. coli and the carbon fixation pathway of C 4 plants (9 -11). A distinctive feature of these systems is that (like the E. coli glyoxylate bypass) each contains a bifunctional enzyme that catalyzes opposing post-translational modifications. These experimental results strongly suggest that bifunctional enzymes tend to confer robustness to biochemical processes. There have been longstanding efforts to understand this connection between bifunctionality and robustness theoretically, beginning with the suggestion by Russo and Silhavy (12) that the robustness of the EnvZ/OmpR system is due specifically to the bifunctionality of EnvZ and the Batchelor and Goulian model (10) of that system. More recently, a desire to explore the bifunctionality-robustness connection from a quantitative, systems biology perspective has motivated the development of mathematical models of the EnvZ/OmpR, IDH regulatory, nitrogen assimilation, and carbon fixation systems (13)(14)(15)(16)(17). Each of these studies relies on mathematical and numerical analysis of the underlying chemical reaction networks to derive expressions that predict robust behavior.
The recent study by Shinar et al. (15) of IDH robustness is representative of these systems-level analyses of bifunctional enzymes. They proposed a detailed chemical reaction network to describe how IDHKP interacts with IDH to regulate the glyoxylate bypass (Fig. 1A). At the center of their network is the ternary complex EI p I, where E denotes IDHKP and I and I p denote active and inactive IDH, respectively. The putative formation of such a ternary complex is a direct consequence of having a bifunctional enzyme that can simultaneously bind two molecules of its substrate. Using mass-action kinetics, Shinar et al. (15) derived a mathematical expression that shows that the ternary complex gives rise to approximately robust behavior in the IDH regulatory system. On the basis of this result, they concluded that their model provides an explanation for the experimental observation of robustness by LaPorte et al. (8).
The reaction network of Shinar et al. (15) (Fig. 1A) makes two major assumptions about the biochemistry of the IDH regulatory system. Formation of the ternary complex EI p I requires that one molecule of IDHKP be able to simultaneously bind one phosphorylated and one unphosphorylated IDH, which implies that IDHKP must have separate active sites for its kinase and phosphatase activities. Additionally, although IDH is known to form an obligate homodimer in vivo (18), Shinar et al. (15) opted to model it as a monomer. When Shinar et al. proposed their model, the crystal structure of IDHKP had not been determined, and the exact configuration of its active sites was an open question. However, detailed structural information about IDHKP has since been published (19), revealing that the enzyme employs a single site for both its kinase and its phosphatase activities and switches between them using a conformational shift (Fig. 1B). It follows from these data that the ternary complex proposed by Shinar et al. (15) cannot form and that their model needs to be amended. In this study, we show that these structural data and other biochemical information about IDH and IDHKP lead to a different model of the E. coli IDH regulatory system. Our model assumes that the IDHKP kinase and phosphatase reactions are catalyzed from the same active site and that dimerization of IDH is an important feature of the system. From mass-action analysis of our model, we derive a polynomial expression, an invariant, that relates the steady-state levels of the phospho-forms of the IDH dimer. The key invariant is obtained using algebraic techniques that do not require assigning numerical values to any of the parameters (i.e. the rate constants in the system). This approach allows us to develop a model that accurately reflects the biochemical complexity of the IDH regulatory system. We then study the behavior of the invariant numerically, given several experimentally justified assumptions about the system. Under these conditions, we find that the invariant predicts the robustness of active IDH to changes in total IDH levels, in agreement with existing experimental results.

EXPERIMENTAL PROCEDURES
Evidence for a Shared Active Site on IDHKP and Dimerization of IDH-Here we outline the biochemical justification for the two major assumptions (shared active site on IDHKP and dimerization of IDH) of our model. Prior to the determination of the crystal structure of IDHKP, a number of mutagenesis studies suggested that the kinase and phosphatase activities might share a single active site. For instance, Stueland et al. (20) found that mutation of a lysine residue in the ATP binding domain severely compromised both the kinase and the phosphatase activity of IDHKP. This result suggested that the two  (15) proposed the above chemical reaction network to describe how IDHKP (E) interacts with active (I) and inactive (I p ) IDH to regulate the E. coli glyoxylate bypass. IDHKP can form two different binary complexes with IDH, EI, in which E acts as a kinase, and EI p , in which E acts as a phosphatase. Two molecules of IDH can also form a ternary complex (EI p E) with IDHKP. To simplify the analysis of their model, Shinar et al. (15) assumed that EI p E does not have phosphatase activity and can only form in an ordered fashion (i.e. EI p can bind to I, but EI cannot bind to I p ). B, structure of two IDHKP molecules (blue/cyan and red/orange) bound to an IDH dimer (green and yellow) obtained by Zheng and Jia (19). Their structural results indicate that a ternary complex can be formed between two IDHKP molecules and one IDH dimer, but not between one IDHKP and two IDH molecules, as in the above network.
activities rely on the same ATP binding site, making it unlikely that the active sties could be located in structurally distinct regions. (This hypothesis was also supported by a photoaffinity labeling study (21).) A kinetic study of two IDHKP mutants provided further evidence for the existence of a shared active site (22). The V max for the phosphatase reaction was shown to be inversely related to the kinase V max , suggesting that the activities are closely linked and possibly differentiated by a conformational shift (22). Supporting these earlier experiments, the IDHKP structure of Zheng and Jia (19) conclusively demonstrated that the active site is shared. IDHKP was found to have two domains, a regulatory domain and a catalytic domain where ATP binds and the kinase and phosphatase activities take place. AMP, a potent inhibitor of IDHKP kinase activity, binds between the domains and induces a conformational change to promote phosphatase activity. This structural information reveals that it is impossible for an IDHKP molecule to act simultaneously as a kinase and a phosphatase.
Additionally, E. coli IDHKP is known to bind only to intact, dimeric IDH (23). This absolute substrate specificity (which has recently been characterized on a structural level (24)) is so strong that replacing E. coli IDH with Bacillus subtilis IDH, which dimerizes and exhibits extensive sequence and structural homology with E. coli IDH, reduces the K m of IDHKP 60-fold for kinase activity and 3,450-fold for phosphatase activity (25). In addition to characterizing isolated IDHKP, Zheng and Jia (19) also reported the structure of an IDH homodimer with an IDHKP bound to each half of the dimer. This structure revealed that the substrate recognition loop of IDHKP interacts extensively with both halves of the dimer, providing an explanation for why IDHKP fails to recognize monomeric IDH. These data strongly suggest that dimerization of IDH is essential to its function in the E. coli glyoxylate bypass.
Model of IDH Regulation in E. coli-We now describe our model of the IDH regulatory system and the corresponding reaction network (see Fig. 2). The model makes several standard assumptions. We assume that all reactions follow a Michaelis-Menten scheme. Accordingly, we do not explicitly account for molecules of ATP or the ATPase activity of IDHKP. ATP is assumed to be maintained in excess at a constant concentration, and levels of ADP and P i are also assumed to be constant. These assumptions are discussed in detail in Ref. 26. We treat IDH as a homodimer that can take one of three forms, I oo (neither monomer phosphorylated), I op (one monomer phosphorylated), and I pp (both monomers phosphorylated), and we assume that all IDH is dimerized (18). Because no cooperative effects are known to occur between the two halves of the IDH homodimer (27), we consider each phosphorylation to be an independent event and do not distinguish between the monophosphorylated forms (i.e. I op ϭ I po ). Each IDH monomer can bind to one IDHKP (denoted as E), leading to the formation of three different ternary complexes: EI oo E, EI op E, and EI pp E. Thus, in our notation, the enzyme-substrate complex studied by Zheng and Jia (19) is EI oo E. Each of these ternary complexes has one (dimeric) IDH attached to two molecules of IDHKP, the opposite of the ternary complex in the model of Shinar et al. (15), which consisted of one IDHKP bound to two units of IDH. Our model includes all possible binary and ternary complexes that can form between IDHKP and the various states of IDH. Additionally, IDHKP is allowed to act as either a kinase or a phosphatase when bound in a ternary complex. Because the two halves of the IDH dimer are independent, we assume that simultaneous double modifications (e.g. EI oo E going to EI pp E ϩ E or EI pp E going to EI oo E ϩ E) are unlikely. We exclude such reactions from the model.
We also assume that the conformational change required for IDHKP to switch between kinase and phosphatase activity occurs on a much faster time scale than the catalytic steps (k cat ϭ 9.8 ϫ 10 Ϫ2 s Ϫ1 for kinase activity and k cat ϭ 6.5 ϫ 10 Ϫ2 s Ϫ1 for phosphatase activity, calculated from data in Ref. 22). Therefore, we do not include a reaction representing the shift between the IDHKP kinase and phosphatase activities. The two conformations are assumed to be in rapid equilibrium before the formation of each binary or ternary complex. Although the kinetics of the IDHKP shift have not been studied, similar conformational changes for other enzymes are known to occur many orders of magnitude faster than the IDHKP catalytic steps, supporting our assumption (28).
Mass-action Analysis and Algebraic Calculations-Mass-action analysis of the reaction network in Fig. 2 yields a system of 10 ordinary differential equations (listed in Table 1). At steady state, these 10 differential equations give 10 polynomial equations, f 1 ϭ 0, . . . , f 10 ϭ 0, that involve the 10 species. For clarity, the species are written as x 1 , . . . , x 10 (defined in Table 1) in Table 1 and in the supplemental material.
Conclusions are based on mathematical analysis of the polynomial system and calculations involving the derived algebraic expressions. Symbolic algebraic and numerical calculations were performed using Mathematica 8.0 (Wolfram Research). A printout of a Mathematica notebook detailing the calculations is provided in the supplemental material. The actual notebook is available on request from the authors. Computation of reaction network deficiency was performed in MATLAB 7.12.0 (The MathWorks, Inc.) using a custom program developed by Guy Shinar and Eran Eden (Reaction Network Analyzer V1.0).
Gröbner basis calculations were performed following the procedure described by Manrai and Gunawardena (29). To derive Equation 2, a lexicographically ordered Gröbner basis was computed for the 10 polynomial equations in Table 1 with the variables ordered as Parameter Estimation-Parameters for the numerical calculations were determined from kinetic data on IDHKP obtained by Miller et al. (22). The parameters in the model (k 1 , . . . , k 20 ) fall into six categories, on-rate, off-rate, and catalytic rate for IDHKP acting as a kinase, and the same for IDHKP acting as a phosphatase. We assumed that the on-rates are determined by the diffusion-limited interaction of IDH and IDHKP and set both at 1 ϫ 10 8 M Ϫ1 s Ϫ1 , a typical value for protein-protein diffusion-limited interactions (30). We set the catalytic rate constants to the k cat values measured for IDHKP kinase and IDHKP phosphatase (22). Finally, we calculated the off-rates from the on-rates, catalytic rates, and Michaelis-Menten constants of Miller et al. (22) according to the formula

RESULTS
Derivation of the Key Invariant-As described under "Experimental Procedures," we applied steady-state mass-action kinetics to the reaction network in Fig. 2 to obtain the system of polynomial equations in Table 1. From algebraic analysis of this polynomial system, we find that there is a polynomial expression, an invariant, involving only these three variables that holds in any positive steady state of the system (i.e. any steady state in which all 10 variables take strictly positive values). It has the form where a 2 , . . . , a 6 are algebraic combinations of the rate constants, normalized so that the coefficient of I oo 2 is 1. The full expressions for the coefficients are given in the supplemental material. The above invariant provides important constraints on the levels of the three IDH phospho-forms at steady state and sheds light on the experimentally observed robustness of the IDH regulatory system, as discussed in the following section. Equation 2 is calculated by a method previously described in Ref. 29, as shown in the supplemental material. Briefly, by algebraically combining the polynomial equations, f 1 ϭ 0, . . . , f 10 ϭ 0, it may be possible to eliminate some of the variables and derive equations that only depend on a specified subset of variables x i1 , . . . , x ik . Gröbner bases provide a systematic way to do this and can be thought of as a polynomial generalization of Gaussian elimination for linear systems (29). Crucially, such elimination can be done with the rate constants treated as uninterpreted algebraic symbols that do not need to be given numerical values. The coefficients of the equations are then rational algebraic expressions in the symbolic rate constants. For instance, it might be possible to obtain a linear equation of the form I oo ϭ K, where K is some rational algebraic expression in the rate constants. Such an equation would imply that the concentration of unphosphorylated IDH is constant at steady state, irrespective of how much total IDH (T) or total IDHKP is present. This is "absolute concentration robustness" (ACR), as defined by Shinar and Feinberg (31). Although a demonstration of ACR would account for the robustness found experimentally, we show using Gröbner bases that there is no such equation involving only I oo (supplemental material). Instead, the invariant in Equation 2 emerges as the simplest algebraic relationship between the three IDH phospho-forms (I oo , I op , and I pp ) that holds in any positive steady state.
Numerical Prediction of Robustness from the Main Invariant-In this section, we show numerically that, given several reasonable assumptions, our invariant predicts the robustness of I oo to T. For the simulations, we use parameters estimated from the kinetic data of Miller et al. (22) as described under "Experimental Procedures" (Table 2). Our assignment of parameter values assumes that reactions are not influenced by the phosphorylation state or bound IDHKP on the other half of the IDH dimer. As such, reactions in a particular class (e.g. all the kinase catalytic reactions) are assumed to be chemically similar and are given the same rate constant (k cat , for the kinase catalytic rate). This reduces the number of independent parameters from 20 to six. We first consider how the invariant constrains I oo at various I op and I pp values with the parameters fixed as defined in Table  2. Fig. 3A shows a plot of I oo as a function of the two other phospho-forms. I oo varies approximately linearly with I op but remains nearly constant at any fixed I op even as I pp changes. This behavior can be confirmed from examination of ∂I oo /∂I pp . As shown in Fig. 3B, ∂I oo /∂I pp is close to zero for all I op and I pp values, confirming that I oo remains constant when I op is fixed.
These results suggest that our model predicts the robustness of active IDH when I op remains constant. The limited experimental data on monophosphorylated IDH appear to support this requirement. Borthwick et al. (32) attempted a systematic characterization of IDH extracted from E. coli growing on glycerol. They readily isolated both fully active (doubly unphosphorylated) and fully inactive (doubly phosphorylated) IDH, but found no monophosphorylated IDH. The authors speculated that the absence of detectable I op might be due in part to disproportionation of the phosphate during purification. Such behavior, however, would imply that the phosphates in singly phosphorylated dimers are markedly more labile than those in doubly phosphorylated dimers, which survived purification in abundance. Given that the halves of the IDH dimer are known to be independent, it appears unlikely that the phospho-state of one half could influence the stability of a phosphate on the other half. We believe instead that the inability of Borthwick et al. (32) to isolate monophosphorylated IDH reflects that very low but nonzero levels of I op are present at steady state. I op must be present in a nonzero concentration because otherwise I pp could only be formed through simultaneous double phosphorylations, which is also unlikely. This view is further supported by the identification by LaPorte et al. (8) of only two phosphoforms during their overexpression experiments and the apparent absence of any physiological role for monophosphorylated IDH. On the basis of this evidence, we suggest that I op is maintained at a consistently low (and therefore approximately constant) level, which we denote as ⑀.
We can now make use of the conservation law for IDH species. Cellular IDHKP levels have been measured to be roughly 1000-fold lower than IDH levels (32, 33). As such, we can write the approximate conservation law FIGURE 2. Proposed model of the IDH regulatory system incorporating substrate dimerization. In our model, IDH is treated as a homodimer that can bind two molecules of IDHKP, whereas IDHKP is assumed to have a single active site and be incapable of simultaneously binding two molecules of IDH. IDH and IDHKP can form a variety of binary and ternary complexes, all of which are represented in the network. No restrictions are placed on the activity of IDHKP in the ternary complexes or the order of complex formation. Accordingly, this is the most complete network that is consistent with what is currently known about the biochemistry of IDH and IDHKP.
T ϭ I oo ϩ I op ϩ I pp (Eq. 3) in which the complexes between IDH and IDHKP have been neglected. Using Equation 3, we write I pp in Equation 2 in terms of total IDH (T), I oo , and I op , set all I op terms to ⑀, and solve for I oo . This procedure yields a long expression for I oo as a function of only T, given in the supplemental material. As shown in Fig.  3C, I oo varies minimally with T assuming the parameter values given in Table 2 and some fixed value for ⑀. In contrast, using the conservation law to eliminate I oo from the invariant and I pp as a function of T reveals that I pp is linearly dependent on total IDH. These computational results are in agreement with the experimental data on IDH robustness (Fig. 3D) and demonstrate that our model provides an explanation for the observed robustness.
The Role of the Rate Constants in the Robustness-Inspection of Equation 2 suggests that robustness arises for the following reason. When the rate constants are assigned the values given in Table 2, ͉a 5 ͉ and ͉a 6 ͉, are at least 30 times larger than ͉a 2 ͉, ͉a 3 ͉, and ͉a 4 ͉ ( indicates that I oo will remain constant at any fixed I op regardless of the value of I pp , in agreement with the behavior observed from inspection of Fig. 3A. The coefficients a 2 , . . . a 6 are complicated polynomials, normalized by a 1 , in the six independent parameters, k on , k off , and k cat for the kinase and on , off , and cat for the phosphatase. We recall that a polynomial is a sum of monomials. For a 2 , a 3 , and a 4 , we found on inspection that each monomial in the corresponding polynomial has a relatively small numerical value when evaluated at the reference parameter values in Table 2. All monomials lie in the range [Ϫ200, 200], and no one monomial is dominant. In contrast, the coefficients a 5 and a 6 each have a pair of monomials that are much larger (by a factor of nearly 50) than all the others. Setting aside the two monomials in each of these pairs, the remaining monomials in a 5 and a 6 all lie in the range [Ϫ200, 200]. Taking the two dominant monomials in a 5 yields the expression k cat ␣, and doing the same in a 6 yields the expression Ϫ cat ␣ where

TABLE 1
Mass-action analysis of the IDH regulatory system Table 1 lists the 10 ordinary differential equations that describe the behavior of the chemical reaction network underlying our model. At steady state, these differential equations all equal zero and become the system of algebraic equations that is studied in detail in this paper. For clarity, the 10 species and complexes appearing in the reaction network are written as x 1 , …, x 10 in the table and in the supplemental material, and the 10 steady-state polynomial equations are referred to as f 1 , …, f 10 . At the reference parameter values in Table 2, the expression k cat ␣ equals 10,786, and the expression Ϫ cat ␣ equals Ϫ7,154, which are within 1% of the exact values of a 5 (10,872) and a 6 (Ϫ7,174), respectively. We see from this that the dominance of a 5 and a 6 over the other coefficients a 2 , a 3 , and a 4 , which formed the basis for the explanation of robustness offered above (Equation 5), arises largely from the structure of ␣, in which the on-and off-rates play a particularly significant and highly multiplicative role. We have not been able to interpret ␣ further, but this analysis does relate the origin of the robustness to the biochemical rate constants. Sensitivity Analysis-We undertook a sensitivity analysis to investigate whether the robust behavior demonstrated in the previous sections is due to the exact numerical values chosen for the parameters. We found that the six parameters, if allowed to vary, do not have to be confined to a smaller dimensional space around the reference values in Table 2 for robustness to be preserved. Mathematically, this indicates that the region of parameter space in which robustness is found is an open subset of six-dimensional space, at least in the vicinity of the reference values.

TABLE 2 Parameter values used in the simulations
We independently and randomly chose values for the onrates, off-rates, and catalytic rates from the uniform probability distribution on the interval [v/2, 2v], where v is the reference value for the corresponding parameter given in Table 2. We did this 200,000 times. We then calculated for each set of parameter values the Michaelis-Menten constants of the kinase and the phosphatase using Equation 1. We found that if the chosen parameter sets were further restricted to those whose Michaelis-Menten constants lay in the interval [0.8v, 1.2v], which reduced the number of parameter sets from 200,000 to ϳ18,000, then robustness was preserved. The interval [0.8v, 1.2v] was the maximum range for the Michaelis-Menten constants for which robustness was obtained with all parameter combinations tested. For each parameter set selected as described, we calculated the mean ∂I oo /∂I pp over a grid of I op and I pp values and found it to be very close to zero in each case (i.e. on the order of 10 Ϫ4 ). This behavior corresponds to what we found at the reference parameter values.
Restricting parameter sets to a range of Michaelis-Menten constants still serves to establish that the region of robust parameters is six-dimensional and hence that the robustness of I oo to changes in T is not due to an accident in the choice of the reference values. The Michaelis-Menten constants are known from experimental measurements, and restricting them to have physiologically realistic values seems appropriate. However, we did find that there is less latitude for change in the Michaelis-Menten constants, which need to be maintained within a smaller interval around their reference values than do the on-, off-, and catalytic rates.
We note one further issue regarding the invariant in Equation 2. Because the invariant is quadratic, there could be two positive solutions for I oo as a function of I op and I pp , depending on the values of the coefficients a 2 , . . . , a 6 . This would be dif- The straight line that interpolates between the widely separated data points is only suggestive. However, the simulation in C shows straight line behavior in the model. IDH activity was determined by monitoring reduction of NADP; one unit of IDH corresponded to the reduction of 1 mol of NADP/min (8). Table 3 gives the numerical values for the six coefficients in Equation. 2 normalized to the value of a 1 when the rate constants have the values specified in Table 2. ficult to interpret physiologically for the IDH regulatory system. The phenomenon of bistability has been reported previously for systems with two or more phosphorylation sites (34,35). However, bistability usually requires three steady states, of which two are stable (hence "bistable") and one is unstable (36). A second positive solution for the invariant may correspond instead to a steady state in which one or more other components have negative values, which would be unphysiological. We found, however, only a single positive solution for I oo at the reference parameter values in Table 2 and in all the parameter sets that we tested in the sensitivity analysis. This behavior can also be attributed to the large negative value for a 6

TABLE 3 Normalized coefficient values for Equation 2
Large negative a 6 values for biochemically reasonable parameters ensure that the above inequality is satisfied and a single positive solution for I oo is always found.

DISCUSSION
Connection to the Shinar-Feinberg Theorem-Shinar and Feinberg (31) recently proved an important theorem for determining when a biochemical network exhibits robustness. Their result hinges on the deficiency of the network, a non-negative integer that is a measure of the dynamical complexity of the network and depends only on its structure and not on its rate constants. Their theorem states that any network that has deficiency one and meets several other conditions will exhibit ACR (described above). The deficiency of the network describing our model (Fig. 2) is two, so the Shinar-Feinberg Theorem cannot be applied to analyze its robustness. The Shinar-Feinberg Theorem is useful because it enables one to study the robustness of a biochemical system without undertaking extensive algebraic calculations of the sort performed in this investigation. However, the failure of our network to meet the required conditions highlights the need for additional theoretical techniques, like Gröbner basis calculations, to characterize robustness and other systems-level properties in complicated biochemical networks.
Comparison with Other Models-We also considered two hypothetical models of the IDH regulatory system that relax one of the major assumptions of our model (dimerization of IDH and bifunctionality of IDHKP). Algebraic analysis of these simplified reaction networks reveals that they are described by fundamentally different invariants than the invariant (Equation 2) derived from our model. Treating IDH as a monomer regulated by a bifunctional IDHKP leads to the simple reaction network shown in Fig. 4A for the interconversion of active (I) and inactive (I p ) IDH. A Gröbner basis calculation on the corresponding reaction network leads to the following invariant (see supplemental material) As such, the model with monomeric IDH predicts that the ratio of active to inactive IDH is kept robust at steady state (i.e. the ratio is set exclusively by parameter values). As seen in Fig. 3D, this behavior is the opposite of what is observed experimentally; the measured ratio of active to inactive IDH dropped from 1.2 to 0.02 when total IDH was increased from 1.3 to 34.5 units (8). Equation 8 suggests that a bifunctional enzyme (that like IDHKP has a shared active site) acting on a monomeric substrate can ensure that the substrate is partitioned proportionately between the modified and unmodified forms even as total enzyme or total substrate levels fluctuate. This type of robustness is potentially useful in certain biological contexts, but is clearly incompatible with the IDH regulatory system, which relies on phosphorylating more and more IDH to maintain constant IDH activity when total IDH is increased. A second hypothetical model treats IDH as a dimer but splits IDHKP into a monofunctional kinase (K) and monofunctional phosphatase (F). Suppose that the two catalytic domains K and F are placed on the same polypeptide to form K ϫ F, in such a way that the two domains do not influence each other (i.e. substrate may be bind to either domain with the same rates, irrespective of what is bound to the other). We have shown in a separate study 3 that the dynamical behavior of such a bifunctional K ϫ F is identical to that of the two separate, monofunctional enzymes K and F, with the single difference that there is an additional conservation law that K tot ϭ F tot . We show that this kind of bifunctionality leads to the property of "robust upper bounds," which does not correspond to the robustness observed experimentally in the IDH system. Furthermore, in this bifunctional scenario, nothing prevents the formation of a ternary complex in which both K and F are simultaneously bound to their respective substrates. Such behavior is not possible with IDHKP, as discussed earlier. Hence, we can rule out a model in which K and F operate as monofunctional enzymes.
We have already discussed the differences between our model and the model of Shinar et al. (15) in detail. Fig. 4B shows one other previously developed model of the IDH regulatory system. Miller et al. (22) proposed a kinetic model of IDHKP that describes the relationship between its kinase and phosphatase activities and the role of ATP, ADP, and pyrophosphate in both reactions. Their model is biochemically realistic and is useful for understanding certain aspects of the mechanism of IDHKP. However, their network, which has deficiency zero, cannot explain ACR of active IDH. This conclusion can be drawn from another recent theorem of Shinar and Feinberg (37), that no network with deficiency zero can exhibit ACR. Additionally, we were unable to find an invariant leading to approximate robustness for the network. These results provide further evidence that the dimerization of IDH, which is not accounted for in the model of Miller et al. (22), is important to the robustness of the system.
Should IDHKP Be Treated as a Dimer?-Purified IDHKP is known to dimerize in vitro through the formation of a disulfide bond between two cysteine residues (33,38). Nevertheless, we opted to model IDHKP as a monomer. IDHKP has not been shown to form a homodimer in vivo, and there is no clear indication that dimerization of IDHKP is physiologically important. This is in sharp contrast to IDH, whose dimerization is essential for the binding of IDHKP. The complex of two IDHKP monomers bound to an IDH dimer was crystallized in the presence of the strong reducing agent dithiothreitol, which was added to prevent the formation of any disulfide bonds that might impede crystallization (39,40). Accordingly, a functional EI oo E complex was able to form even when IDHKP dimerization was inhibited. We speculate that the observed IDHKP homodimers formed because large numbers of IDHKP monomers were brought into close contact under favorably oxidative conditions during purification and that dimerization is unlikely to occur inside cells. Formation of new disulfide bonds in vitro is not unusual; dimerization of hepatitis B surface antigen produced in yeast is one example (41).
Regulation of IDHKP-Our model focuses on the regulation of IDH by IDHKP to understand how robustness in the glyoxylate bypass is implemented. To focus the model on IDH regulation, we assume that the kinase and phosphatase conformations of IDHKP are in rapid equilibrium and that both are therefore available as needed for reactions. This approach, however, ignores the regulation of IDHKP (e.g. how the cell increases kinase or phosphatase activity in response to a change in carbon substrate availability). IDHKP is regulated by a variety of cellular metabolites (including AMP), and the structural basis for the action of many IDHKP regulators is now understood (19). An interesting avenue for future research is the development of an integrated, dynamical model that explains how total IDHKP is partitioned between kinase and phosphatase activity by small-molecule regulators and how this partitioning enables responsiveness to changing environmental conditions. Experiments in which cells are subjected to rapid, oscillatory changes in carbon substrate could provide additional information for such a model (42,43).
Experimental Tests of the IDH Regulatory Model-Further experimental work is needed to test all the predictions of the model and to fully characterize IDH robustness. Detection of a rare phospho-form such as I op should now be possible with the advent of mass spectrometry-based proteomics (44). Levels of all three phospho-forms should be accurately quantified using mass spectrometry in E. coli strains with and without IDH overexpressed, enabling confirmation that the phospho-form relationship and robustness predicted from the invariant are correct. These measurements should be made at more total IDH levels than were considered by LaPorte et al. (8). Other experimental tests could be used to demonstrate directly that robustness is dependent on the bifunctional nature of IDHKP. For instance, robustness could be studied in a strain in which IDHKP has been deleted and replaced with a monofunctional kinase and monofunctional phosphatase. A comparable experiment has been performed with the bifunctional adenylyltransferase that regulates glutamine synthetase in E. coli (16).
Conclusions-We have developed a model of the E. coli IDH regulatory system that is firmly grounded in the current structural, mechanistic, and kinetic data on IDH and IDHKP and that condenses much information about this classic biochemical system into a succinct quantitative framework. We have used new methods of algebraic analysis, together with experimentally supported assumptions, to show how the experimental observations of robustness by LaPorte et al. (8) can be explained. Our model suggests that two features of the system, the bifunctionality of IDHKP and the dimerization of IDH, are important to achieving robustness. This approach highlights the value of incorporating detailed enzymatic data, which are often neglected in modeling studies, into quantitative systems biology models.
Many complex biological models are studied by numerical simulation of the underlying ordinary differential equations. This approach leaves all conclusions completely dependent on the numerical values assigned to parameters. In the case of mass-action models, these parameters are the rate constants of the individual reactions, which unlike phenomenological kinetic constants (e.g. K m or k cat ) generally cannot be measured directly. We have previously termed the reliance of many systems biology models on hard-to-measure constants the "parameter problem" (45). Our approach to deriving the main invariant circumvents the parameter problem by treating all constants symbolically and studying the model through exact mathematical analysis, not simulations. This method enabled us to obtain a simple expression for the relationship between the three IDH phospho-forms without making assumptions about parameter values. We then used the invariant to make specific, numerical predictions about the robustness of active IDH and the requirements for maintaining robust behavior. These results would have been difficult or impossible to arrive at just by simulating the parent system. The exact characterization of a complicated biochemical model (through derivation of the main invariant) was made possible through the application of recently developed algebraic techniques. Future work should aim to identify other analytical tools that can be exploited to efficiently study systems of biological interest without simulations.