The aspartate‐derived amino‐acid pathway from plants is well suited for analysing the function of the allosteric network of interactions in branched pathways. For this purpose, a detailed kinetic model of the system in the plant model Arabidopsis was constructed on the basis of in vitro kinetic measurements. The data, assembled into a mathematical model, reproduce in vivo measurements and also provide non‐intuitive predictions. A crucial result is the identification of allosteric interactions whose function is not to couple demand and supply but to maintain a high independence between fluxes in competing pathways. In addition, the model shows that enzyme isoforms are not functionally redundant, because they contribute unequally to the flux and its regulation. Another result is the identification of the threonine concentration as the most sensitive variable in the system, suggesting a regulatory role for threonine at a higher level of integration.
A challenging goal of systems biology is to develop detailed kinetic models for simulating and predicting the dynamic responses of metabolic networks. However, very few published models have been based on kinetic measurements on enzymes in conditions relevant to those in vivo. In addition, the few kinetic models of real systems that have been published do not address the action of effectors in branched pathways. As feedback regulation often nullifies the effects of genetic manipulations, it is important to understand the molecular mechanism of robustness to genetic perturbations of biochemical networks. Until now the functioning of allosteric controls in a branched metabolic system has only been studied theoretically (e.g. Savageau, 1974; Cornish‐Bowden et al, 1995). In addition, branch‐point enzymes often exist as isoforms that respond unequally to cooperative inhibition (or sometimes activation) by allosteric effectors. The physiological role of isoforms has been discussed for many years (Ureta, 1978), and study of the flux distribution between isoforms in a real pathway should shed light on the need for them. The aspartate‐derived amino‐acid pathway from plants constitutes an excellent model system for understanding regulatory mechanisms in branched metabolic pathways. This pathway is responsible for the distribution of the carbon flux from aspartate into the branches for synthesis of lysine, threonine, methionine and isoleucine. There are several branch‐points, many enzyme isoforms and different allosteric control mechanisms (inhibition, activation, antagonism and synergism), as illustrated in Figure 1B. Much is known about the individual components of the systems but little about their role in the economy of the ensemble, as until now no large‐scale detailed kinetic model of this system has been available either in plants or in microorganisms, though small parts of it have been modelled in bacteria (Chassagnole et al, 2001; Yang et al, 2005) and in plants (Curien et al, 2003).
A mathematical model of the core of the pathway in the chloroplasts of the model plant Arabidopsis has been constructed with a view to understanding and quantifying the short‐term regulatory capabilities of the system. Flux values and metabolite concentration calculated with the model (Figure 3) were very close to in vivo values. In addition, the metabolic pattern of several mutants could be reproduced. The present work represents the first detailed kinetic model of a real branched metabolic pathway in which all the allosteric controls and isoforms are taken into account. Kinetic characterization in conditions that mimic the in vivo metabolic context, together with the use of purified recombinant enzymes, explains the predictive power of the model. A key result from this model is the identification of allosteric interactions whose function is not to couple demand and but to maintain a high independence between fluxes in competing pathways. Specifically, it shows that proximal controls of dihydrodipicolinate synthase (DHDPS) by Lys and of threonine synthase (TS1) by S‐adenosylmethionine (AdoMet) are involved in the coupling between demand and supply. Distal controls (synergistic inhibition of aspartate kinase 1, inhibition of aspartate kinase 2 by lysine) may also contribute to this function but are much less efficient. Their main function is rather to attenuate flux changes in the branches where demand remains constant, and they thus participate in the independence between pathways. The existence of isoforms with different regulatory patterns also contributes to this independence indicating that enzyme isoforms are not functionally redundant. Another result is the explanation of threonine concentration instability, suggesting a regulatory role for threonine at a higher level of integration. Finally, important information useful for later work is the estimation of the half‐time of the system (400 s). This metabolic system is a slow responsive one, but kinetic controls are still slightly faster than protein concentration control mechanisms.
A key result from the quantitative model of the aspartate pathway from Arabidopsis is the identification of feedback controls whose function is not to couple demand and supply but to maintain a high independence between fluxes in competing pathways.
The model shows that enzyme isoforms are not functionally redundant, because they contribute unequally to the flux and/or to its regulation
Another result is the identification of threonine concentration as the most sensitive variable in the system, suggesting a regulatory role for threonine at a higher level of integration
A challenging aim of systems biology is to develop detailed kinetic models for simulating and predicting the dynamic responses of metabolic networks. However, very few published models have been based on kinetic measurements on enzymes in conditions relevant to those in vivo. For example, the model database at http://jjj.biochem.sun.ac.za now contains nearly 90 models, but most of these are stoichiometric rather than kinetic models, and the others include only four systems that have been modeled kinetically with realistic measured parameters: glycolysis in the human erythrocyte (Mulquiney and Kuchel, 1999), glycolysis in Trypanosoma brucei (Helfert et al, 2001), methionine/threonine metabolism in Arabidopsis (Curien et al, 2003) and sucrose metabolism in sugarcane (Uys et al, 2007).
Furthermore, apart from our earlier study in Arabidopsis, the few kinetic models that have been published do not address the action of effectors in branched pathways. As feedback regulation often nullifies the effects of genetic manipulations, it is important to understand the molecular mechanism of robustness to genetic perturbations of biochemical networks. Until now the functioning of allosteric controls in a branched metabolic system has only been studied theoretically (Savageau, 1974; Cornish‐Bowden et al, 1995). In addition, branch‐point enzymes often exist as isoforms that respond unequally to cooperative inhibition (or sometimes activation) by allosteric effectors. The physiological role of isoforms has been discussed for many years (Ureta, 1978), and study of the flux distribution between isoforms in a real pathway should shed light on the need for them. The aspartate‐derived amino‐acid pathway from plants constitutes an excellent model system for understanding regulatory mechanisms in branched metabolic pathways. This pathway is responsible for the distribution of the carbon flux from Asp into the branches for synthesis of Lys, Thr, Met and Ile (Figure 1). There are several branch‐points, many enzyme isoforms and different allosteric control mechanisms (inhibition, activation, antagonism and synergism), as illustrated in Figure 1B. Until now, no large‐scale detailed kinetic model of this system has been available either in plants or in microorganisms, though small parts of it have been modeled in bacteria (Chassagnole et al, 2001; Yang et al, 2005) and in plants (Curien et al, 2003). The pathway has great biotechnological importance for improving the nutritional quality of plants (Galili and Hofgen, 2002; Galili et al, 2005; Azevedo et al, 2006), the synthesis of secondary metabolites (Brown et al, 2003) and applications to green chemistry (Slater et al, 1999). A quantitative understanding of flux control in the network of the Asp family of amino acids in plants is a prerequisite to a rational engineering of this pathway.
A mathematical model of the core of the pathway in the chloroplasts of the model plant Arabidopsis has been constructed with a view to understand and to quantify the short‐term regulatory capabilities of the system. Kinetic data from the literature are generally inadequate for kinetic modeling (Cornish‐Bowden and Hofmeyr, 2005), and to avoid the problems with unrealistic parameter values we have combined kinetic measurements, obtained in vitro with purified recombinant enzymes, with in vitro reconstitution and numerical simulation. Care was taken to measure enzyme activity in a reconstituted metabolic context, thus maintaining biological relevance. This model has allowed a detailed understanding of the roles of the isoforms and their regulation.
Description of the metabolic system
The core of the Asp‐derived amino‐acid pathway in Arabidopsis (Figure 1B) is a system of 13 enzymes that distribute the carbon flux from Asp into the branches leading to the synthesis of Lys, Thr, Met and Ile. These amino acids are used for protein synthesis, and may also be degraded or used as precursors for the synthesis of secondary metabolites. Met is also the precursor of AdoMet, the universal donor of methyl groups.
There are various complications for understanding the pathway. Monofunctional aspartate kinase (which exists as two isoforms, AK1 and AK2), bifunctional aspartate‐kinase–homoserine‐dehydrogenase (two isoforms, AKI–HSDH I and AKII–HSDH II), dihydrodipicolinate synthase (two isoforms, DHDPS1 and DHDPS2), threonine synthase and threonine deaminase are allosteric enzymes that are co‐expressed in leaf chloroplasts. The other enzymes present in the chloroplast—the non‐allosteric enzymes, aspartate semialdehyde dehydrogenase, homoserine kinase and cystathionine γ‐synthase and the allosteric threonine deaminase—do not have isoforms.
The activity of the allosteric enzymes is modified by various effectors (Figure 1B), and binding of these increases or decreases enzyme activity by modulating the substrate‐binding or catalytic parameters, or both (Curien et al, 2008). Thr, Lys and AdoMet are involved in dual controls, each controlling two different catalytic steps. Synergistic inhibition exists for one enzyme (AK1), inhibited by Lys and AdoMet (Curien et al, 2007). There are three different cases of activation: TS1 shows pure activation by AdoMet (Curien et al, 1998); activation of TD by Val results from displacement of Ile inhibitor, with no effect of Val by itself (Wessel et al, 2000); and the AK activity of the bifunctional enzymes AKI–HSDH I and AKII–HSDH II is stimulated by Ala, Ile, Cys, Ser and Val (Curien et al, 2005). These activators, which compete with Thr, are stimulatory even in the absence of the inhibitor Thr. The isoforms of the allosteric enzymes differ in their regulatory and kinetic properties (values for 50% activity of the effectors and kinetic constants are in Supplementary data, section S1).
The complexity of this metabolism prevents an easy understanding of the function played by any given enzyme isoform relative to another. In addition, the functional importance of the cross‐pathway activation by AdoMet and the dual controls by Lys, Thr and AdoMet cannot be clearly understood without a quantitative analysis. The mathematical model was thus intended to simulate the integrated behavior of the system described in Figure 1 in conditions mimicking the metabolic context of photosynthesizing chloroplasts in the leaf. The modeling procedure is summarized in Materials and methods. The experimental data for enzyme purification and kinetics are grouped in Supplementary data, together with validation of the enzyme rate equations using in vitro reconstituted enzyme systems. The simulations described focus on understanding the functional role of isoforms and allosteric controls in the whole system.
The only unknown parameters in the model are the cellular demands for the amino acids, and the following procedure was followed to adjust these. Lys, Thr and Ile were initially assumed to be consumed exclusively by aminoacyl tRNA synthetases, but without taking into account fluxes of amino‐acid degradation (see below). Aminoacyl tRNA synthetases show Michaelis–Menten behavior with respect to their cognate amino acids, and Km values for these are known experimentally (equations (S14)–(S16), Supplementary data). The only parameter to adjust is the apparent limiting rate for lysyl, threonyl and isoleucyl tRNA synthetases (VAaRS), the same value being used for all three, because the abundance of Lys, Thr and Ile in plant proteins is similar. Methionyl tRNA synthetase is not included, because the flux in the Met branch is controlled by AdoMet (not Met), and [AdoMet] was set at its physiological value.
In the absence of a demand for Lys, Thr and Ile there was no steady state, and these metabolites accumulated without limit: this is clearly visible in Figure 2A for Thr, but is also true for Lys and Ile. The feedback mechanisms shown in Figure 1B are thus insufficient to maintain the stability in these conditions. For VAaRS close to 0.11 μM s−1 a steady state was possible, but with amino‐acid concentrations high and non‐physiological (data not shown); increasing the demand decreased these concentrations. For VAaRS=0.43 μM s−1, the steady‐state variables reached values close to physiological ones (Figure 2B and C). This steady state, called the reference state, is described in detail below. The half‐time of the system in the reference state, i.e. the time for [Thr] to reach 0.5[Thr]ss, was about 400 s (6.7 min) and concentrations were practically independent of time after about 8000 s (2.2 h). Direct calculation of the steady state confirmed that the state reached after this time was a true steady state, satisfying the tests of stability carried out by COPASI.
Model steady‐state characteristics
As a first validating step, concentrations and fluxes were compared with known concentration and flux values measured in planta. The variables in the model in the reference state are summarized in Figure 3. The steady‐state concentrations, [Lys]ss=69 μM, [PHser]ss=45 μM, [Thr]ss=303 μM and [Ile]ss=59 μM, are consistent with the subcellular concentrations in the chloroplast compartment, [PHser]=80 μM (Gakiere et al, 2000), [Lys]=70 μM, [Thr]=300 μM and [Ile]=50 μM (Curien et al, 2005). [AspP]ss, [ASA]ss and [Hser]ss are in the micromolar range, as expected from the impossibility of detecting them in soluble plant extracts. The total flux is 1 μM s−1, in agreement with fluxes measured in vivo by NMR, such as a flux of Asp consumption of 37 nmol h−1 per 106 cells measured in cultured tomato cells (Rontein et al, 2002). This converts into 1.7 μM s−1, assuming a cellular volume of 65 pl (Winter et al, 1994a) and taking into account the chloroplast localization of the synthesis of Asp‐derived amino acids, as chloroplasts represent 9.5% of the cellular volume (Winter et al, 1994a). A similar flux value was obtained by radioactive measurements in Lemna plants (Giovanelli et al, 1989).
In the model, the fluxes of Lys, Thr and Ile toward protein synthesis are almost equal (about 0.32 μM s−1) and the flux in the Met branch is 5.5‐fold lower, again in close agreement with the in vivo data of Giovanelli et al (1989). (The conversion of flux values of Giovanelli et al (1989), expressed in nmol per frond per doubling time into μM s−1 values is detailed in Curien et al (2003): Lys, 0.5 μM s−1; Thr → protein, 0.5 μM s−1; Thr → Ile, 0.4 μM s−1; Met, 0.14 μM s−1.) The metabolic pattern of several plant mutants could be reproduced by the model (see Discussion), adding to its credibility. In the light of the current data, the model seems realistic and can be used to provide deeper understanding of the system.
Metabolic control analysis of the reference state
Elasticities in the reference state are shown in Table I. The numerically large elasticities for AspP and ASA in the ASADH reaction result from the fact that this reaction is close to equilibrium in the steady state. The large elasticity for Ile in the TD reaction has a different origin, because of the high cooperativity (with Hill coefficient of 3) with respect to Ile.
The flux control coefficients in the reference state, for the common flux for the entire pathway (ASADH) and for the four demand branches, are shown in Table II, and concentration control coefficients are shown for all intermediates in Table III. Most of the control of the common flux is shared between AK1 and the three aminoacyl tRNA synthetases, the other enzymes having very little control. Flux to methionine is controlled almost completely by its own reaction (CGS), but there is a large negative flux control coefficient for the competing reaction TS1: this is an example of the branch‐point effect (LaPorte et al, 1984) and is explained by the much larger flux through TS1, so a small variation in TS1 flux inevitably affects the CGS flux. Each of the three protein‐forming demand steps has a large control coefficient for its own flux, indicating that the regulatory mechanisms efficiently transfer control of flux from the supply steps (where it is undesirable) to the demand steps (where it is needed), and minimize unwanted competition between these steps. Note, however, the substantial contribution of AK1 to the Thr → protein flux, which is a significant cause of the Thr instability considered below. This is also visible in the large concentration control coefficient of AK1 for [Thr]. Note also that HSK has no detectable control over any flux or any concentration apart from that of Hser, its substrate. This exception is inevitable for a reaction that has no flux control, because its effects can only be exerted on the concentrations of its substrate and product.
Flux distribution between isoforms at the reference state
The results in Figure 3 indicate that enzyme isoforms at a given step contribute to different extents to the flux. AK1, which is only weakly inhibited by a concentration of Lys sufficient to give 50% inhibition of AK2, and has the highest activity among the four AKs (Curien et al, 2007), accounts for 73% of the flux at the AK step (15% for AK2). The Thr‐sensitive enzymes contribute less, because they have lower kcat and higher Km values than the Lys‐sensitive AKs. AKI contributes less (2.7%) than AKII (9%), because it requires higher concentrations of the activators and is less activated (Curien et al, 2005). It follows from these differences that the fraction of the total flux controlled by Thr at the AK step is much lower than that controlled by Lys, with important consequences for the variation of [Thr] in conditions of low demand for products, as discussed below.
At the ASA branch‐point, DHDPS1 and DHDPS2 differ in quantitative importance as they account for 10% and 90%, respectively, of the total steady‐state flux in the Lys branch. This is explained by inhibition of DHDPS1 at values of [Lys] lower than those needed to inhibit DHDPS2 (see Supplementary Figure S7B). In the competing Hser branch, HSDH I and II contribute nearly equally to Hser flux. However, HSDH II is virtually insensitive to physiological [Thr] (see Supplementary Figure S6B). Thus, half of the flux at this step escapes control by Thr. Here again, a functional difference exists between the two isoforms.
Analysis of allosteric control function
To determine the function of individual allosteric controls, the response of the system was examined by simulating biologically pertinent situations.
Modulation of the demand for Met and AdoMet.
Met demand can be modified in vivo for the synthesis of secondary metabolites (Brown et al, 2003), when it is degraded (Rebeille et al, 2006; Goyer et al, 2007) and presumably transiently when the flux in the Met–AdoMet methylation cycle (which does not result in a net consumption of Met) is modified. Changes in Met demand alone are thus biologically relevant, the sensor of Met metabolism in plants being AdoMet. Changing [AdoMet] in the model simulates modified demand for Met or AdoMet. The results in Figure 4A indicate that when [AdoMet] increases (thus simulating a decreased demand) flux in the Met branch decreases. Thus, the dual control by AdoMet (TS1 activation and AK1 inhibition) provides a very efficient negative feedback mechanism for the Met branch. Until now, the role of AdoMet in controlling the system was unclear (Kocsis et al, 2003). However, the model shows that in the range 2–30 μM for [AdoMet], flux in the Met branch can be decreased 10‐fold, with only an 8% decrease in Lys flux, a slight increase followed by a slight decrease in Thr flux, making a 15% decrease overall, and similar slight changes in Ile flux, making an 8% decrease overall. Plots in Figure 4 allow the observed changes to be understood. An increase in [AdoMet] decreases AK1 activity (Figure 4B) and thus the total flux in the system, and hence the concentrations of the downstream metabolites (ASA, Figure 4C and PHser, Figure 4D) decrease. The decrease in [PHser] leads to a lower flux in the Met branch (Figure 4A). The Met flux is highly sensitive to AdoMet, whereas the flux in the Thr branch remains nearly constant (Figure 4A), because the increase in [AdoMet] (TS1 activator) compensates, almost exactly, for the decrease in [PHser], PHser being the substrate of TS1. The origin of the bell‐shaped curve for [Thr] (Figure 4D) is explained in the legend to Figure 4. Flux changes in the Thr branch are attenuated because the AKI and AKII activities decrease or increase in response to changes in [Thr] (reinforcement or relaxation of the inhibition, Figure 4B). Changes in Thr → protein flux are buffered further compared with the changes in [Thr] because threonyl‐tRNA synthetase is nearly saturated by Thr.
The flux in the Lys branch decreases very little when [AdoMet] increases, because the effect of the decrease in [ASA] (Figure 4C) is almost exactly counteracted by a relaxation of the inhibition of DHDPS and AK2 by Lys in the Lys path (Figure 4D). Changes in the Lys → protein flux are further buffered by the near‐saturation of lysyl‐tRNA synthetase by Lys.
Note that [Lys] decreases (Figure 4D) when [AdoMet] increases. As Lys and AdoMet are co‐inhibitors of AK1, but their concentrations vary in opposite directions, the interaction of Lys and AdoMet with AK1 does not show synergy at the system level when the demand for Met–AdoMet is modified.
AdoMet targets two different enzymes, with both proximal control (TS1 activation) and distal control (AK1 inhibition). (The second isoform TS2 is omitted from the model because it is not expressed in chloroplasts.) The relative importance of these can be examined. Considering TS1 allosteric activation first, the coupling between TS1 activity and [AdoMet] can be abolished in the model by replacing [AdoMet] in the TS1 equation by a constant value of 20 μM, i.e. the value of [AdoMet] in the reference state, thus maintaining the TS1 activity at the same level as in the reference state. Fluxes and concentrations in the reference state are then identical in the original and modified versions of the model, allowing a comparison of the response to changes in [AdoMet]. When TS1 is uncoupled from AdoMet, simulations indicate that the negative feedback effect of AdoMet on flux in the Met branch is virtually abolished (data not shown). Thus, the allosteric activation of TS1 is central to the control of Met flux by AdoMet, and its loss cannot be compensated for by the synergistic effects of AdoMet and Lys on AK1. Indeed, when the coupling is suppressed, and [AdoMet] is increased, the decrease in [PHser] because of AK1 inhibition is no longer compensated for by the activation of TS1. The decrease in input flux is distributed between the Met and Thr branches, and as the flux in the Met branch is much smaller than that in the Thr branch, the flux change in the Met branch is strongly attenuated. If the control of TS1 by AdoMet is essential, then what is the function of the distal control (AK1 inhibition by AdoMet)? Somewhat unexpectedly, if the coupling between AK1 activity and AdoMet level is suppressed ([AdoMet]=20 μM in the AK1 equation), the relation between [AdoMet] and flux in the Met branch is virtually identical to that in Figure 4 (data not shown). However, uncoupling AK1 activity and [AdoMet] modifies the relation between [Thr] and [AdoMet]. Rather than the bell‐shaped curve of Figure 4D, with a small amplitude of [Thr] changes, there is a monotonic increase in [Thr]ss from 130 μM at 2 μM AdoMet to 345 μM at 40 μM AdoMet (data not shown). Thus, control of AK1 by AdoMet attenuates [Thr] changes when demand in the Met branch is modified; therefore, the control of AK1 by AdoMet maintains independence between the Met and Thr pathways, whereas control of TS1 by AdoMet is a negative feedback.
Modulation of the demand for Lys.
Lys can be degraded by LKR, and a 12‐fold increase in the abundance of this enzyme occurs in hormone‐treated Arabidopsis plants (Stepansky and Galili, 2003). Therefore, simulating a change in [LKR] (i.e. a change in demand for Lys with no change in demand for the other products) is biologically relevant. [LKR] can be increased as much as to 30 μM from its normal value in vivo in the range 0–1 μM, with very little change in the rest of the system (Figure 5A). At 25 μM LKR (a rather high protein concentration), the flux of Lys degradation is of the same order of magnitude as the rate of Lys consumption by lysyl‐tRNA synthetase. The decrease in [Lys]ss because of the increased consumption by LKR is attenuated by relaxation of the inhibition of DHDPS and AK by Lys. Flux in the Hser branch is affected to a much smaller extent when [LKR] is increased, changing by about 0.8% in the range 0–25 μM LKR. This insensitivity results from a good coordination of Lys‐sensitive AK and DHDPS activities with respect to Lys inhibition, leading to very small changes in [ASA]ss (Figure 5C). In the Thr branch, the bell‐shaped curve in Figure 5E ([Thr]ss) is a consequence of the same shape of curve for [ASA]ss (Figure 5C, inset), as explained in the legend to Figure 5.
The relative contribution of the control by Lys at the DHDPS step (proximal control) and at the AK step (distal control) after changes in the demand for Lys can now be examined. The activities of DHDPS1 and DHDPS2 can be uncoupled from [Lys] in the model (variable [Lys] in DHDPS1 and DHDPS2 equations replaced by a fixed factor of 69.4 μM, equal to the value of [Lys]ss in the reference state). If [LKR] is increased in that case, the amplitude of [Lys]ss changes is 33 μM (curve labeled ‘unregulated’ in Figure 5B), 11 μM more than the range of 22 μM for the normal model (‘regulated’ curve). The difference is not large because the relaxation of inhibition of Lys‐sensitive AK1 and AK2 by Lys compensates in part for the lack of control at the DHDPS step. The major difference concerns [Thr]ss, which then depends nearly linearly on [LKR] and increases three‐fold over the range of [LKR] used, in contrast to the very small changes seen with the unmodified model (Figure 5D). Thus, the control of DHDPS by Lys has two major functions, to adapt supply to demand in the Lys branch and to increase the independence between the Lys and Thr pathways. Conversely, when the activities of AK1 and AK2 are uncoupled from variations in [Lys] (in the same way as described above for the DHDPS1 and DHDPS2 equations) and [LKR] is increased, the quantitative changes in [Lys] are only slightly more pronounced (decreased from 69 to 44 μM instead of from 69 to 48 μM, Figure 5B). However, there was a large decrease (data not shown) in [Thr]ss, from 300 to 190 μM, very different from the bell‐shaped curve, with an amplitude of 15 μM as shown in Figure 5E. Distal controls of AKs by Lys are thus mainly needed for homeostasis, buffering flux changes in the Hser branch when demand for Lys is modified. This conclusion was definitely not obvious.
Modulation of the demand for Thr.
THA is an enzyme involved in the degradation of Thr in plants (Jander et al, 2004), and indirect evidence from the model and from in vivo data suggests that its concentration in vivo is in the range 0–1 μM. (Thr does accumulate in vivo, but could never accumulate if [THA] were higher than 1 μM. This was confirmed by the model.) An increase in [THA] was simulated to examine the consequences of a modulation of Thr degradation in vivo. As shown in Figure 6A, [THA] can be increased up to about 8 μM with little effect on the flux in the Lys branch (Figure 6A). At 7.5 μM THA, the flux of Thr degradation reaches the same value as the flux of Thr → protein synthesis. [Lys]ss increased by 5.2% and the Lys flux by 1.2%. This increase in Lys flux is due to the increase in total flux caused by relaxation of the inhibition of the Thr‐sensitive AK activities by Thr (as [Thr] decreases). The increase in flux in the Lys branch (Figure 6A) is modest because of the efficient inhibition of DHDPS by Lys. At 7.5 μM THA, [Thr]ss is 51% lower than in the absence of THA, and the flux of Thr → protein is 21% lower. The results indicate that [Thr] homeostasis is maintained less efficiently than that of [Lys] (compare Figure 5A with Figure 6A). We will return later to the higher sensitivity of [Thr]ss to these modulations.
Modulation of the demand for protein synthesis.
The above simulations indicate that allosteric controls allow a specific redirection of flux in the demand pathway while maintaining fluxes relatively unchanged in the branches where demands remain constant. What would happen if demand is modified in a simultaneous manner in all the branches? This would correspond to a change in the rate of protein synthesis, and, as Lys, Thr and Ile have similar abundance in plant proteins, it can be simulated by modifying the demand in the three branches in the same proportion, defining the same value for the limiting rates VAaRS of lysyl, threonyl and isoleucyl‐tRNA synthetases in equations (S14)–(S16). Figure 7A and B show that when demands for Lys, Thr and Ile are modified simultaneously, quantitative changes are much more important for Thr (Figure 7A) than for Lys and Ile (Figure 7B). A two‐fold decrease of VAaRS leads to a marked increase of [Thr]ss, whereas [Lys]ss increases to a much lesser extent and [Ile]ss only slightly. When VAaRS is less than about 0.11 μM s−1, a steady state is no longer possible in the system, as Thr accumulates without limit. The system is thus not very robust to decreases in the global demand, and [Thr]ss shows the greatest instability. The half‐time also increases considerably when demand decreases, from about 400 s (6.7 min) in the reference state (VAaRS=0.2 μM s−1) to about 33 500 s (9.3 h) when VAaRS=0.43 μM s−1 (Figure 7D), because although the increase in [Thr]ss is large (Figure 7A), it occurs rather slowly.
Origin of the instability of the steady‐state concentration of Thr
Instability of [Thr]ss is a property not only of the model but also of the natural system, as it has been observed in planta (see Discussion), and can be explained by the control coefficients in Table III: two demand steps (Thr → protein and Ile → protein) have substantial negative control over [Thr], so [Thr]ss must fall sharply if these are increased. In addition, AK1, which carries the bulk of the AK flux but is not inhibited by Thr, has a large positive control over [Thr], and so any small variations in AK1 activity have magnified effects on [Thr]ss. Lys has numerically large elasticities for the two DHDPS isoenzymes (Lys branch) than those of Thr for the HSDH isoenzymes (Table I), which means that regulation by Lys is much more stringent than regulation by Thr, and so excess flux is diverted preferentially into the Hser branch. Downstream, at the PHser branch‐point, flux in the Met branch is much smaller than Thr flux, and the excess flux is diverted into the Thr branch. Ile synthesis is strongly controlled by the allosteric inhibition of TD by Ile (with an elasticity of –2.3). In addition, demand of Thr for protein synthesis is nearly saturated, so any increase in the synthesis of Thr has almost no effect on the flux, but a large increase in [Thr]ss (Figure 7A).
Importance of reaction reversibility
Reaction reversibility can be an important consideration in metabolic models (Hofmeyr and Cornish‐Bowden, 1997; Cornish‐Bowden and Hofmeyr, 2005), although in reality it is not reversibility as such that must be taken into account, but sensitivity to inhibition by reaction products or, better, downstream metabolites (Cornish‐Bowden and Cárdenas, 2001).
In this model, reversibility is experimentally detectable for only three enzymes (AK, ASADH and HSDH, Supplementary Figures S1, S3 and S6). The activity of HSDH in the non‐physiological direction was ignored because it is quantitatively negligible under normal conditions (see Supplementary Figure S6). To appreciate the importance of the reversibility of AK and ASADH under normal conditions, reactions in the reverse directions were suppressed from the model, producing only minor changes to the dynamics of the system, the steady‐state variables and the regulatory responses. It follows that the strong feedback controls in the model are sufficiently effective in transferring information to the early steps to render the effects of reversibility insignificant.
This study represents the first detailed kinetic model of a real branched metabolic pathway in which all the allosteric controls and isoforms are taken into account. Kinetic characterization in conditions that mimic the in vivo metabolic context, together with the use of purified recombinant enzymes, explains the predictive power of the model. Specifically, it shows that proximal controls of DHDPS by Lys and of TS1 by AdoMet are involved in the coupling between demand and supply. Distal controls (synergistic inhibition of AK1 and inhibition of AK2 by Lys) might also contribute to this function but are much less efficient. Their main function is rather to attenuate flux changes in the branches where demand remains constant, and they thus participate in the independence between pathways. The existence of isoforms with different regulatory patterns also contributes to this independence. The presence of an AK sensitive to Lys alone, together with an AK sensitive to Lys and AdoMet, implies that part of the Asp flux escapes control by AdoMet. This helps to increase the independence between the Lys and the Met–AdoMet branches. In the same manner, the existence of an isoform of HSDH nearly insensitive to physiological [Thr] (HSDH II), together with one that is sensitive (HSDH I), contributes to the independence between the Met–AdoMet and Thr branches, because part of the flux escapes control by Thr.
The explanation of the specific role played by AdoMet is an important contribution of the model. The physiological significance of the allosteric activation of TS1 has long been unclear. Kocsis et al (2003) described the metabolic perturbations in an Arabidopsis plant with an insertion mutation in the gene of methionine S‐methyltransferase (an enzyme involved in the S‐methylmethionine cycle). A 1.15‐fold increase in [AdoMet] and a two‐fold decrease in [Thr] were measured. In mto3, another mutant deficient in AdoMet synthetase (Shen et al, 2002), [AdoMet] decreased 1.5‐fold compared with the wild‐type plant and [Thr] increased 1.4‐fold. This ‘apparent lack of an in vivo relationship between AdoMet level and TS activation’ surprised the authors (Kocsis et al, 2003), who expected a positive correlation. However, for [AdoMet] ≈20 μM (physiological concentration) the changes in [Thr]ss are predicted by the model with a good accuracy: when [AdoMet] was increased 1.15‐fold [Thr]ss decreased 1.1‐fold, and when [AdoMet] was decreased 1.5‐fold [Thr]ss increased 1.3‐fold. The agreement between the model and the plant mutants constitutes an additional and independent strong validation of the model.
The efficiency of the controls by Thr is less than those by Lys and AdoMet, even though the inhibition of AK and HSDH by Thr attenuates changes in [Thr]ss when demand is modified. [Thr]ss is thus the most sensitive variable in the system to perturbations. This striking property of the model is also seen, in less pronounced form, in the natural system: in the Arabidopsis DHDPS2 knockout plant (Craciun et al, 2000), [Lys] decreased 1.5‐fold, whereas [Thr] increased 4.6‐fold; when DHDPS2 is suppressed in the model [Lys]ss decreases 2‐fold and [Thr]ss increases 23‐fold. Similarly, in an Arabidopsis plant mutant in which Lys‐sensitive AK was desensitized with respect to Lys inhibition (Heremans and Jacobs, 1995), [Lys] increased 1.5‐fold but [Thr] increased 6‐fold; in the model, desensitization of AK1 to Lys (replacing [Lys] by a constant of zero in the AK1 equation), leads to a 1.5‐fold increase in [Lys]ss but a 12‐fold increase in [Thr]ss. The sensitivity of [Thr] to perturbations suggests that it might be a good candidate for higher‐level control processes, i.e. processes outside the limits of the model. Moreover, its position at the start of a major biosynthetic pathway makes it very suitable for such a role of integrating all of the signals coming from the different amino‐acid products. Experimental support for such a role is provided by antisense potato plants with reduced levels of TS and of [Thr], which showed higher [Asp] than normal (Zeh et al, 2001). Conversely, tobacco cells transformed with a TS from E. coli had 5.7‐fold higher [Thr], and showed a 10‐fold decrease in the rate of incorporation of 14C‐aspartate (Muhitch, 1997). If Thr does indeed directly or indirectly regulate [Asp], this would be an efficient way to adjust demand and supply because AK activities are proportional to [Asp] in the physiological range (see Supplementary Figure S2). The mechanism of such a putative control remains to be determined.
At the level of integration examined here some questions remain. One of the puzzles concerns the need for isoforms such as DHDPS1 and the AKI domain of AKI–HSDH I, which both make apparently negligible contributions to the flux, and removing them from the model makes very little difference to the variables in the reference state (data not shown). The normal‐looking aspect of the DHDPS2 knockout plant mentioned above (Craciun et al, 2000) suggests that isoforms contribute to the robustness of the system. Elimination of DHSDP2 from the model results in a slight decrease in [Lys]ss and a slight increase in [ASA]ss, and the resulting increase in DHDPS1 activity compensates to a large extent for DHDPS2 deficiency (data not shown).
A second puzzle concerns the activation of Thr‐sensitive AKs by Ala, Cys, Ile, Ser and Val, all of which are efficient activators. Why do Thr‐sensitive AKs need activators and why are there so many different ones? Intuitively, the activation of Thr‐sensitive AK by various amino acids, among which some, Val and Ile, are products of the pathways, seems as a way to increase sensitivity of [Thr] to changes in the demand for protein synthesis. However, preliminary results (data not shown), using more complex equations for AKs than those used in the model, show that Thr efficiently displaces the activators bound to Thr‐sensitive AKs, and so any amplifying effect is quantitatively limited. The function of Thr‐sensitive AK activation thus remains an open question that requires integration at a higher level of organization.
A third puzzle concerns the role played by the antagonism between Val and Ile in their control of TD activity. As the control of TD by Val is distal, it might be involved in stabilizing [Ile] when demand in the Leu/Val branches is specifically modified. Extension of the present model to incorporate branched‐chain amino‐acid metabolism in plants should help to clarify this point.
Finally, the reason why AK and HSDH activities are present on the same polypeptide is an old question that has not had any answer. Although bifunctionality might be related to metabolic channeling, in vitro reconstitution (see section S2 in Supplementary data) as well as in vivo data (Giovanelli et al, 1989) exclude a direct transfer of ASA between ASADH and HSDH, nor is bifunctionality justified by a structural coupling at the enzyme or system level. Indeed, kinetic analysis (Curien et al, 2005) indicated a high degree of independence between the AK and the HSDH domains. In the model, the AK domains of AK–HSDHs contribute only 12% to the total flux, whereas the HSDH domains contribute about 70% (Figure 3). The advantage of bifunctionality probably needs to be explained in terms of coordination of changes in AK and HSDH protein concentrations, minimizing effects on the concentrations of the branch‐point metabolite ASA.
As many amino‐acid concentrations are modified by perturbations to plant metabolism, our model will be useful for interpreting metabolic effects of changes to the genome or the environment. This work also provides the mathematical basis for a more complete quantification of flux distribution during the life of a plant once the magnitude of enzyme concentration changes becomes known in detail and the putative control of Asp supply by Thr deciphered. The results for Arabidopsis are transposable to other plants, as the allosteric control pattern and the presence of isoforms are not specific to Arabidopsis. For example, Lys‐sensitive and Lys/AdoMet‐sensitive AKs are found in other plants (Rognes et al, 1980; Azevedo et al, 1992), and all plants characterized to date have AdoMet‐activated TS and Thr‐sensitive AK–HSDH. This might contrast with the case of hexokinase isoforms in vertebrates, for which wide variations in distribution are found even between quite closely related animals (Ureta, 1982), but it might simply reflect a less exhaustive study of isoforms in plants. Finally, the instability of [Thr] is also observed in other plants, such as Lemna paucicostata (Giovanelli et al, 1988) and maize mutants (Muehlbauer et al, 1994). The model should therefore be useful for understanding Asp‐derived amino‐acid metabolism in other plants and predicting the effects of modifications, especially for optimizing the production of compounds of economic interest in crop plants without affecting essential cellular processes.
Analysis of the integrated behavior of a biological system is a powerful tool for revealing the specialized roles of elements that seem redundant in superficial examination of a map, such as the one in Figure 1B, and thus explaining why such complexity has been conserved in the course of evolution.
Materials and methods
Kinetic characterization of isolated enzymes and enzyme quantification
Recombinant enzymes were purified to homogeneity as described previously (Curien et al, 1998, 2005, 2007; Ravanel et al, 1998; Wessel et al, 2000; Paris et al, 2002a, 2002b, 2003) or according to the procedures summarized at the end of Supplementary data (purification of DHDPS1,2 and HSK). Pure proteins were quantified as described by Scopes (1974). Enzyme concentration is expressed on a monomer basis. Kinetic experiments were carried out with purified recombinant enzymes. The kinetic characterizations were focused on the response of the enzymes to the substrates specific to the pathway (others, such as ATP, being maintained at physiological concentrations) and to the various allosteric effectors. Assays were carried out at 30°C in the presence of 50 mM HEPES pH 8.0, 50 mM NaCl, 10 mM MgOAc, 100 mM KCl and 50 mM KOAc (buffer A). Under these conditions, pH and salt concentrations were close to physiological: [Cl−]=96±11 mM, [K+]=194±6 mM, [Na+]=72±8 mM (Robinson and Downton, 1984); [Mg2+]=10 mM (Portis and Heldt, 1976), pH=8.0 (Heldt et al, 1973). Slight quantitative differences with our kinetic data published earlier result from differences in the experimental conditions, closer here to the physiological operating conditions of the enzymes. All the enzyme ligands were present at concentrations corresponding to a physiological metabolic state of reference (photosynthesizing leaf chloroplasts, Table IV). The concentrations of the ligands of interest (internal metabolites, allosteric effectors) were then varied. We assumed that the general metabolites aspartate, ATP, ADP, NADPH, NADP, Pi and pyruvate have their own homeostatic regulatory mechanisms, and treated their concentrations as external quantities with fixed values (1.5 mM aspartate, 2 mM ATP, 0.5 mM ADP, 0.2 mM NADPH, 0.2 mM NADP, 10 mM Pi and 1 mM pyruvate). The modeling of the more complex enzymes was thereby considerably simplified without loss of biological relevance. Detailed information on each enzyme activity and the derivation of simplified rate laws are in Supplementary data. Enzyme concentrations in the stroma of isolated chloroplasts from Arabidopsis (Ferro et al, 2003) were determined by ELISA (Table V) with antibodies raised against the purified recombinant proteins. Purified recombinant proteins were used as standards. Owing to cross reactions (AK2 antibodies with AK1; DHDPS1 antibodies with DHDPS2; AKII–HSDH II antibodies with AKI–HSDH I), the ELISA gave an estimation of the total amount of the monofunctional AKs, of the bifunctional AKs and of the DHDPSs, but did not give access to the relative amount of each isoform. Western blots (Supplementary Figure S14A) indicated that AK1 and AK2 are co‐expressed in leaf chloroplasts and are present at similar abundance. DHDPS1 and DHDPS2 are also co‐expressed in Arabidopsis chloroplast stroma at similar levels (Supplementary Figure S4B). Molecular data for the two bifunctional AK–HSDHs indicate that the two genes are simultaneously expressed in leaves (Rognes et al, 2003), and in the absence of more detailed information we assumed identical concentrations for these proteins.
Reconstitution of enzyme systems
To check whether the enzyme rate equations adequately accounted for the behavior of the enzymes under approximately physiological conditions (i.e. enzymes working in concert, high enzyme concentration, low substrate concentrations), systems of enzymes were reconstituted and the observed behavior was compared with the prediction of models corresponding to the experimental conditions. Owing to the good agreement between the experiments and model predictions the experiments were confirmatory and are not detailed in the main text (see Supplementary Figures S11 and S12).
Model assembly and simplifications
Data on the nine enzymes operating in the upper part of the Asp‐pathway (AK1, AK2, AKI–HSDH I, AKII–HSDH II, ASADH, DHDPS1, DHDPS2) were compiled with the model of the Met–Thr branch‐point validated earlier (Curien et al, 2003), including HSK, TS1 and CGS. The model was further expended with the enzyme TD, specifically characterized for the present study (Supplementary Figure S10). The demand in the amino acids for protein synthesis is determined by the rate at which aminoacyl‐tRNA synthetases (AaRS) form aminoacyl‐tRNA by condensation of an uncharged tRNA with its cognate adenylated amino acid. We did not distinguish between cytosolic or chloroplastic AaRS, as this is not relevant to the internal properties of the Asp‐derived amino‐acid pathway examined here. A simple Michaelis–Menten equation was used for each lumped (cytosolic plus stromal) AaRS (equations (S15)–(S18)) with amino‐acid concentration as sole variable. Km values for the amino acids were from the literature on plant enzymes when available, or from yeast (Ile‐tRNA synthetases). VAaRS is an apparent limiting velocity in equations (S14)–(S16), and corresponds to the rate at which tRNAs are set free by translating ribosomes. It is used as an adjustable parameter for simulating modification in the demand for the different amino acids for protein synthesis. Lys, Thr and Ile are present in similar proportions in plant proteins and a single value of VAaRS was used for all three amino acids.
Lys and Thr can enter into carbon metabolism on account of the activity of lysine ketoglutarate reductase (LKR) and threonine aldolase (THA), respectively. LKR condenses lysine and α‐ketoglutarate to form saccharopine, which is subsequently transformed in several steps into acetyl CoA (Galili et al, 2001). Kinetic parameters of Arabidopsis recombinant LKR (kcat=3.7 s−1 and KmLys =13 mM) were from Miron et al (2000) (equation S17). THA cleaves Thr into acetaldehyde and Gly, with two isoforms in Arabidopsis that have similar values of KmThy, in the range 4–7 mM (Joshi et al, 2006). The catalytic constants of the plant enzymes are unknown, and the value of 1.7 s−1 determined for the enzyme from yeast was used (Liu et al, 1997) (equation S18). LKR and THA protein concentrations are adjustable parameters in the model. The six enzyme steps involved in the transformation of dihydrodipicolinate into lysine were omitted. The omission is justified because the DHDPS reaction is irreversible and insensitive to its product, so the immediate downstream part of the pathway is isolated from the upstream part.
The model of the Ile branch was reduced to the irreversible deamination of Thr to 2‐ketobutyrate by TD (Figure 1), which is inhibited by Ile, with the inhibition antagonized by Val. [Ile] is an internal variable in the model and [Val] is external (fixed), set at a physiological value of 100 μM (Table IV) to ensure that the TD activation level is physiological. The model does not allow quantification of the interaction between Leu plus Val synthesis and Ile synthesis when the demand in one or the other end‐product is modified. It will be possible, however, as soon as a model of the branched‐chain amino‐acid pathway in plants becomes available.
The model of the Met branch comprises cystathionine γ‐synthase (CGS) enzyme, [Cys] and [AdoMet] fixed, and the branch‐point metabolite [PHser] as internal variable. As Met does not control any enzyme activity in the system shown in Figure 1 (the sensor of Met–AdoMet metabolism in plants being AdoMet), it was not necessary to include the demand of Met for protein synthesis at this step of integration. Modulation of the demand for Met–AdoMet was simulated by a change in [AdoMet] (set at a physiological value, 20 μM, or varied as indicated in the figures).
Domain of validity of the model
Strictly speaking, the domain of validity of the model concerns the leaf mesophyll cell of an Arabidopsis plant in the daytime. As protein concentrations are external parameters in the model, the time scale of its validity is about an hour (i.e. before protein concentration is modified). The domain of validity can, however, be extended, as the model reproduces the phenotypes of several plant mutants (see Discussion). In addition, experimental measurements carried out with reconstituted enzyme systems (data not shown) indicated a low sensitivity to [ATP]/[ADP], [NADPH]/[NADP], [Pi] and [pyruvate] when varied within physiological limits. Together the data suggest that the domain of validity of the model can be extended to various metabolic situations. As we do not know whether [Asp] changes in vivo, this concentration was not modified in the model. It cannot therefore be used to simulate changes in [Asp], which would require much more complex equations for the AKs.
Mathematical model and simulations
The rate laws (Section S1) and the set of differential equations (Section S4) are available in Supplementary data. Simulations, both steady‐state and time‐dependent, were carried out with COPASI (Hoops et al, 2006), and in most cases checked with Berkeley Madonna (http://www.berkeleymadonna.com/), with identical numerical results from both. The SBML file (Hucka et al, 2003) is included in Supplementary data. For time dependence studies the time increment was set at 0.1 s and initial internal variable values were set at zero. The integration method was either LSODA (COPASI) or 4th order Runge–Kutta (Berkeley Madonna).
AK1: At5g13280; AK2: At5g14060; AK3: At3g02020; AKI–HSDH I:At1g31230; AKII–HSDH II: At4g19710; ASADH: At1g14810; DHDPS1: At3g60880; DHDPS2: At2g45440; HSK: At2g17265; CGS: At3g01120; TS1: At4g29840; TS2: At1g72810; TD: At3g10050.
We thank Norbert Rolland for kindly providing us with Arabidopsis chloroplast stroma.
Conflict of interest
The authors declare that they have no conflict of interest.
Supplementary data, Supplementary figures S1‐14, Supplementary tables S1‐2The original version of the supplementary information contained errors in equations S3, S4 and S5. The expression for d[Lys]/dt contained ‐vTHA instead of ‐vLKR. These errors are corrected in this file. The corrections do not affect either the SBML or the Berkeley Madonna listings. We have also made a trivial correction to the last sentence in Section 4.1.The Supplementary Materials file was corrected on 16 February 2010. [msb200929-sup-0001.pdf]
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2009 EMBO and Nature Publishing Group