Simulating plant metabolic pathways with enzyme-kinetic models

knock-out mutants of MAM1 and MAM3 (methylthiolalkylmalate). These enzymes define the length of glucosinolate chains through control of the number of chain elongation cycles. The profiles of the glucosinolates in

Deterministic, ODE-based enzyme-kinetic models have the longest history in the area of metabolic pathway modeling. The local complexity of kinetic models leads to the fact that for reasons of feasibility usually single pathways or even just some reactions are modeled. For decades kinetic models have been the most frequently used mathematical approach to metabolism, starting with the work of B. Chance who published the first numerical simulation of a biochemical system in 1943, solving the equations for the behavior of a simple enzymatic system using a mechanical differential analyzer (Chance 1943).

A model needs a defined purpose
Before starting the construction of a kinetic model it is necessary to define the scope of the model. This definition must be as exact as possible, because it will have a profound impact on all further steps of the modeling process. Once this purpose has been defined, the borders of the model can be set accordingly. Due to the connection between model size and processing complexity, a kinetic model will normally describe just some reactions up to a few metabolic pathways. We will illustrate the generation of a kinetic model by the example of the Higgins- Sel'kov oscillator (Higgins, 1964;Sel'kov, 1968), a simple kinetic model which has been thoroughly analyzed over time and which has been used as a teaching example earlier (Klipp et al., 2005). The modeling process is summarized in Figure 1.

Putting together the parts
The first step for the assembly of a kinetic model is the examination of the stoichiometric relations, thus obtaining information about the network structure. The size of the conceived model is of high importance, because the inclusion of too many reactions and metabolites would make the processing computationally inefficient. The reactions of our example model are: That means, S is imported into the modeled system by R1, converted to P by R2, and finally taken out of the system by R3. The stoichiometric matrix N of the reaction system is the following: With the stoichiometric data and the derived stoichiometric matrix at hand structural problems like inactive reactions or non-balanced metabolites can be identified by structural modeling techniques such as elementary mode analysis (Schuster et al., 1999). An elementary ( in which v i specifies the flux through the i th reaction, k i are parameters indicating the velocities of the reactions, and S and P are the concentrations of the two metabolites. R1 carries a constant flux, while R2 and R3 follow mass action kinetics, in the case of R2 the product P is acting as an activator of the reaction. At this point it is important to account for reversibility as well as inhibition or activation of an enzyme, since omitting these effects is a common cause of unrealistic behavior (Cornish-Bowden and Cárdenas, 2001). Using the rate laws and the stoichiometric matrix, the time dependent differential equations for the metabolites can be established: Provided that the parameters are known, the steady state of the system can be determined by calculating values of the metabolite concentrations that make the right-hand side of each equation equal to zero. While the rate laws in this toy model are very simple, other more realistic models will contain details about enzyme saturation and inhibition, which for example can be included using Michaelis-Menten kinetics (Eq. 6). Eq. 6) Where to get the data from?
For establishing a detailed kinetic model, a variety of data sources have to be consulted in order to obtain detailed and reliable information. The stoichiometry of a reaction is usually well-established and can be taken from standard biochemistry textbooks or from online databases such as the examples listed in Table 1. The structure of the rate law depends on the stoichiometry of the reaction and on the enzyme's mechanism and may be derived from enzyme kinetics textbooks (Segel 1975). The maximal velocity of the reaction (V max ) is dependent on the concentration of the enzyme in the cell or the extract, but rather a result of regulated gene expression, and thus has to be measured for each situation again. We note here that the V max used in kinetic models is not the specific activity which is measured with the purified enzyme. Thus, V max is dependent on the enzyme's concentration and needs to be measured in vivo or in crude extracts. For reversible reactions the V max is different for the forward and reverse reaction, respectively. The binding constant (or Michaelis constant) K m results from the structure of the enzyme, and thus is independent of the enzyme's concentration. We therefore argue that in many cases it is a legitimate simplification to take the K m values from the literature, also from related species if the data situation requires.
The state variables of the system, v(S) and S, i.e. the fluxes and metabolite concentrations, respectively, can be measured for a specific state of the system, usually at (quasi) steady state, but sometimes also as time-course. If enough parameters of the model are known, it is theoretically possible to estimate a limited number of missing parameters by fitting the system variables (metabolite concentrations) to data sets from steady-state or time-course measurements. Care has to be taken as the data used for the parameter estimation is not to be used in the model validation. over a short time period, the concentrations are then updated with these changes, and the process is repeated. Depending on the model and its parameters, different behaviors might be observed. Figure 2 shows the example model characterized by different sets of parameters which results in different behavior. In Fig. 2a the system converges towards a steady-state, which is the case for most biological systems if they are kept under constant conditions. At the steady-state, all fluxes in this linear pathway are the same, so that the concentrations of the metabolites do not vary. In Fig. 2b, a divergent behavior is depicted: one metabolite concentration is constantly increasing. A similar scenario is obtained if one metabolite is depleted with such a high rate that it eventually reaches zero. Consequently, the metabolic network might come to a complete hold. Changing the same parameter in the example model to yet another value results in oscillatory behavior of the system (see Fig. 2c). It has been shown that the probability of a model to exhibit oscillatory behavior or other instability actually increases with its size (May 1973). Consequently, for large models it is an exception rather than common to reach a stable steady state. The reason why large biological systems nevertheless tend to steady states might be that they are the product of evolution.
Modeling tools or other mathematical frameworks are indispensable for building and interrogating a model. To simplify using different tools to analyze one model, data standards and exchange languages have been defined. Finally, models can be deposited in online databases to make them accessible for the scientific community (Table 1).

Steady state: constant fluxes and concentrations
A steady state or stationary state of a metabolic system is defined by the fact that all of the state variables (i.e. metabolite concentrations and fluxes) stay constant in value over time, which means that the system will not leave steady state unless a change of these state variables, a perturbation occurs. A perturbation can be a change in metabolite concentration or enzyme activity. The behavior of a system after a perturbation differs if its steady state is stable (i.e. the perturbed system returns to a steady state), unstable (i.e. the perturbed system leaves steady state) or metastable (i.e. the perturbed system oscillates).
A very useful visualization of steady state conditions is the impact of metabolic fluxes on metabolite concentrations. A metabolite concentration (S) is defined to be changed by two fluxes, synthesis (v syn ) and consumption (v con ). If the steady state at a given metabolite concentration (S 0 ) of this system is stable, the metabolite concentration and fluxes should return to their original status after an infinitesimal perturbation. In case of product inhibition, the synthesis will diminish with increasing product concentration, while the consumption is usually increased (see Fig. 3).
After a perturbation the metabolite concentration will be either higher or lower than the concentration at steady state. This dictates the necessary net flux to reach steady state conditions: if the concentration after perturbation is lower than before (S 0 ), the net flux has to be positive to reach the steady state. Accordingly, if the concentration is higher after perturbation than before (S 0 ), the net flux must be negative. In a system in which the net flux shows a more complicated course with the increase of metabolite concentration, more than one steady state can occur (e.g. Poolman et al., 2000). But, as mentioned before, not all steady states are stable after a perturbation. If the net flux is negative and the concentration at the same time-point is lower than at the steady state (or vice versa), the system will further shift away from this steady state.
A stable steady state represents a local or global minimum; therefore the gradients of all points close to the steady state must be negative. These steady states are comparable to the deepest point of a curved surface, so the system will return to this point after perturbation.
Unstable steady states would correspond for example to saddle points on this surface. This information on the gradients in the vicinity of a steady state can be derived from the Jacobian matrix, which describes the flux changes occurring upon small perturbations of the metabolic system (Klipp et al., 2005;Steuer and Junker, 2009).
After identifying a steady state, the response to a system perturbation can be simulated, e.g. by inducing a change of a metabolite concentration and observing the following return to a steady state (metabolite pulse). Further applications are in silico overexpression or silencing experiments. Effectively, the overexpression or silencing of an enzyme can be simulated by increasing or decreasing, respectively, the maximal activity (v max ) of the enzyme in question.
The expression of an additional enzyme working in the metabolic network is simulated in the same way, for example the simulated expression of sucrose phosphorylase in potato tubers (Junker, 2004).

No bottlenecks
One of the long term applications of metabolic analysis of plants is the change of metabolite levels for product improvement (i.e. increase of biosynthesis, reduction of toxic compounds, etc). For decades there has been the quest for "rate limiting steps" or "bottlenecks" in metabolic pathways. It took a relatively long time to appreciate the theory of Metabolic Control Analysis (MCA) for quantifying the impact of a change in enzymatic activity on a flux and/or a metabolite, which has been already developed in the 1970s (Heinrich and Rapoport, 1974;Kacser and Burns, 1973). The main point of the theory is that the control of metabolic fluxes and metabolite concentrations is distributed across many enzymes of a metabolic pathway, and that this control can be quantified for each enzymatic step. Nowadays it has been internalized that the strategy to over-express one enzyme in order to increase the flux through an associated pathway will most likely fail (Morandini and Salamini, 2003). This is supported by a recent study in which was shown for a certain pathway that maximal catalytic activities (i.e., limiting rates) of most enzymes were in large surplus of the fluxes through these reaction steps (Junker et al., 2007). indicates an enzyme with rate limiting capabilities. However, such a value is rather unusual (Fell, 1992 and1997), meaning that the classic idea of a true rate limiting enzyme is mostly not valid, as stated above.

Related and current approaches
The generation of kinetic models using explicit kinetics is arguably the most detailed depiction of a metabolic system. Nevertheless it also requires values numerous parameters that are usually not easy to determine. Thus, the process kinetic modeling suffers from a notorious lack of information. Several other approaches have been proposed to overcome this problem, usually by simplifying the underlying kinetic formulas without losing the physiological closeness to the in vivo situation.
An earlier approach for a simplified mathematical description of biochemical networks is Biochemical Systems Theory, mainly put forward by M. Savageau and E. Voit (Savageau, 1969;Voit, 2000). In contrast to simple linear mass-action kinetics that fail to accurately account for several biochemical properties of enzymes, BST mainly focuses on the description of a unified mathematical treatment of non-linear rate equations. All metabolic rate equations are represented via power law functions Eq.9) in which α i are apparent rate constants, S j are the system variables, such as metabolite concentrations, while γ ij are the apparent kinetic orders. Using these power-law functions and the mass balances of all metabolites a synergistic system of differential equations ("S-System") can be established. The mentioned kinetic orders correspond with the scaled elasticities used in MCA, which both do not depend on substrate concentrations and will therefore not saturate with increasing substrate concentration. However, a limitation of BST is the necessity of a fully specified reference state, the definition of which is not always possible.
An approach that is recently gaining more and more attention is the lin-log framework (linearlogarithmic kinetics). In this approach, all rate equations are defined by their deviation from a reference state (v 0 , S 0 , E T 0 for the reference reaction rate, metabolite concentration, and total enzyme activity, respectively) depending on a logarithmic change in concentration: Eq. 10) in which ε denotes the reference elasticity (Steuer and Junker, 2009).
The performance of the lin-log kinetics is in certain cases better than Power-Law functions or linear approximation (Heijnen, 2005), because the elasticities and therefore the kinetic orders are not constant in the lin-log framework, but dependent on metabolite concentrations. This approach is not restricted to an analytical solution, but maintains the interpretability of the involved parameters in the metabolic context. However, the lin-log function shares the problem of a required reference state with the BST as well as inaccuracies at small concentration values. However, the latter is sometimes considered to be a lesser problem due to the limited concentration ranges in metabolic pathways (Heijnen, 2005).
Structural Kinetic Modeling (SKM) is a recent approach to close a gap between pure structural modeling derived from stoichiometric data and kinetic modeling, depicting the dynamic values of a metabolic system (Steuer et al., 2006;Steuer and Junker, 2009). The main difference to other approaches is the parametric representation of the Jacobian matrix for all possible parameters. This gives the possibility to access certain dynamic behaviors (e.g. oscillation) without knowing the exact functional form of the rate laws or the values of the parameters.

Examples of kinetic models for plant metabolic pathways
It is not astonishing that the first models of plant metabolic pathways were build for the Calvin cycle, arguably the most important metabolic pathway in plants (reviewed in Poolman et al., 2000). Comprehensive lists of kinetic models of any plant metabolic pathway have been published in recent reviews (Morgan and Rhodes, 2002;Rios-Estepa and Lange, 2007).
Furthermore, various models can be found in model repositories (see Table 1) such as JWS Online (Olivier and Snoep, 2004), which also allows amateur users to explore the modeling and simulation process.
The first model of parts of photosynthesis linked equations for RubisCO kinetics with stoichiometric equations of the photosynthetic carbon reduction cycle and the photorespiratory carbon oxidation cycle (Farquhar et al., 1980). In a recent approach, Zhu and coworkers constructed a large compartmented model of leaf carbon metabolism comprising more than 40 reactions and transporters (Zhu et al., 2007). To simulate the effect of evolutionary adaptation to higher CO 2 concentrations, they allowed the model to adapt by redistributing the amount of protein-nitrogen between the enzymes by means of an evolutionary algorithm. They found a substantial increase in photosynthesis, suggesting that the "typical" partitioning in C 3 leaves might be suboptimal for maximizing the light-saturated rate of photosynthesis.
The sucrose-to-starch transition has been the focus of two kinetic models of potato tuber metabolism: one to study the cold sweetening effect observed when tubers are detached of the mother plant and stored (Assmus, 2005), and the other one to simulate the effect of increased activities in sucrose cleaving enzymes on metabolic fluxes (Junker, 2004). In both cases the models correctly predict the trends observed in experimental data.
In contrast to potato tubers, sugar metabolism in sugar cane operates in the direction of sucrose synthesis. In an earlier model of sucrose accumulation in developing sugar cane (Rohwer and Botha, 2001), each reaction step was modeled with accurate kinetics by using detailed experimental and literature data. The differences between experiment and simulation never exceeded a factor of two. The simulations suggested that overexpression of the fructose or glucose transporter or the vacuolar sucrose import protein, as well as reduction of cytosolic neutral invertase levels appear to be the most promising targets for genetic manipulation. The extended version of the model (Uys et al., 2007) also contains information about growth; so that the model describes the metabolic behavior as sugarcane parenchymal tissue matures in different internodes of the stem.
One of the latest established kinetic models simulates the aspartate-derived amino-acid pathway in Arabidopsis thaliana (Curien et al., 2009). This model was based on in vitro kinetic measurements and successfully reproduces in vivo data like metabolite concentrations and fluxes. The aspartate-pathway includes several branch-points and many enzyme isoforms.
According to the simulation, these isoforms contribute differently to the regulation of the respective fluxes and are therefore not redundant. This is a non-intuitive prediction of the the simulation matched with the measured chain length distributions. Furthermore, the model predicted a connection between the regulations of the two enzyme forms in so far, as the overexpression of MAM3 is critical to achieve the measured chain length profile in MAM1 knock-out plants.

Final thoughts
Models of plant metabolic pathways were already built long before the term systems biology was coined, and will continue to play an important role in systems-wide approaches to plant metabolism. The fast increase in biochemical data quantity and quality will support these approaches in the future. However, the resolution of the data will certainly increase through technologies such as sub-cellular fractionation, rapid sampling, and non-invasive online measurements. A topic that has to be addressed in future is the standardization of input data (Jenkins et al., 2004). Furthermore, a standard for modeling steps, specifically for decomposing complex systems into modules and assembling modules together, needs to be established (Oberhardt et al., 2009). The connection of metabolic models with other layers of regulation is another example for future improvements.
It is generally useful to have a look at kinetic models of other systems, such as animal cells and bacteria in order to estimate future developments in plant kinetic modelling. Specifically the metabolic networks in bacteria are well characterized in comparison to plant networks. In an approach for modeling bacterial growth, Fang et al. started by modeling the inhibitory effect on metabolic fluxes, followed by the translation of this information into the growth rate of a single bacterium. Finally the growth of a whole population of bacteria was simulated to estimate the effective bacteria concentration (Fang et al., 2009). This methodology might be adapted to multicellular organisms such as plants by simulating a population of cells rather than a single "average" cell which is currently usually done.
But despite the aforementioned better knowledge base in bacterial systems and the generation of metabolic networks from genome data it is still necessary to manually construct a kinetic model from a given metabolic network (Durot et al., 2009). This process is further complicated by the patchy character of the experimental data. Therefore, the main task for the future in all the fields connected with systems biology is the reliable and automated acquisition (either experimental or from established databases, see Table 1) and assembly of high-throughput data into new metabolic models, thereby increasing the comparability of the  C: Rate laws defining the reaction rate of each enzymatic step. D: Combining the stoichiometric data with the rate laws into a system of ordinary differential equations describing the change in metabolite concentrations over time.