Simulating Plant Metabolic Pathways with Enzyme-Kinetic Models^1,[C]

THE MATHEMATICS BEHIND METABOLIC NETWORKS

A metabolic network can be translated in mathematical terms by relatively easy means: The concentration of a metabolite is described by a variable S_i. As a result of mass balances (no matter can appear or disappear), the change of this variable over time (dS_i/dt) is given by the sum of the rates of the enzymes synthesizing the metabolite minus the sum of the rates of the enzymes utilizing the metabolite. The rate of each enzymatic step can be described by enzyme-kinetic rate laws, such as the Michaelis-Menten equation, as a function depending on metabolite concentrations and parameters such as the maximal velocity of a reaction, or binding constants. This process yields a system of ordinary differential equations (ODEs) in which dS_i/dt is on one side and the metabolite-dependent rate laws are on the other side of the equations (see Eqs. 4 and 5 below). With this system of differential equation, the metabolic network can be simulated, and by solving the system of ODEs the steady state can be calculated, in which all reaction rates and metabolite concentrations are constant.

The process described above is simplified by considering the cell as a homogenous volume. While for many scenarios this is a valid simplification, for other instances, ODE-based kinetic models will not be the method of choice. For example, if spatial gradients are of importance to address a certain question, partial differential equations would be used instead, which increases the complexity of a model dramatically. Furthermore, if metabolites with very low concentrations are considered (e.g. hormones), stochastic effects are coming into play, so that the model would include stochastic components rather than being deterministic. More details on different methods for metabolic modeling are given in a recent comprehensive overview of computational models of metabolism (Steuer and Junker, 2009).

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 that has been thoroughly analyzed over time and that 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:

(1)

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:

(2)

With the stoichiometric data and the derived stoichiometric matrix at hand structural problems like inactive reactions or nonbalanced metabolites can be identified by structural modeling techniques such as elementary mode analysis (Schuster et al., 1999). An elementary (flux) mode is a stoichiometrically balanced route defining an input-output character of the metabolic network. In our example there is a single elementary flux mode of the system in which all reactions operate at the same rate. Thus, there are no structural problems such as deadlocks. The next step is the detailed description of all reactions involved in the model with rate laws. Again there are different levels of model accuracy. The rate laws for the example model are kept simple:

(3)

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:

(4)

(5)

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).

(6)

WHERE TO GET THE DATA FROM?

For establishing a detailed kinetic model, a variety of data sources have to be consulted 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 I. 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 that 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 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.

SIMULATION: THE MODEL AT WORK

Once the stoichiometry, the rate laws, and all parameters are collected, the equation system of the kinetic model is fully specified. Using specialized software (Table I), the model can thus be deterministically simulated in a stepwise manner. Starting from a defined point, the differential equations are solved to obtain the expected changes in metabolite concentrations 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 that results in different behavior. In Figure 2A the system converges toward 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 Figure 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 I).

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₀) 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₀), the net flux has to be positive to reach the steady state. Accordingly, if the concentration is higher after perturbation than before (S₀), 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 Suc phosphorylase in potato (Solanum tuberosum) 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 (Kacser and Burns, 1973; Heinrich and Rapoport, 1974). 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 overexpress one enzyme 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).

The main tools for the evaluation of a global or systemic effect of a perturbation are the control coefficients. The flux control coefficient is defined as:

(7)

And the concentration control coefficient is defined as:

(8)

These coefficients are the quantitative connections between a change in the maximal catalytic activity (A) of an enzyme and the resulting change of flux (J) or metabolite concentration (S), respectively. So these relations compare the named factors between two different steady states. Control coefficients can be computed given a kinetic metabolic model, but they can also be measured by observing the changes in fluxes or metabolite concentrations obtained upon genetic modification of the activity of an enzyme. A flux control coefficient close to 1 indicates an enzyme with rate-limiting capabilities. However, such a value is rather unusual (Fell, 1992, 1997), 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 and 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 (BST), 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 nonlinear rate equations. All metabolic rate equations are represented via power-law functions

(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 (linear-logarithmic kinetics). In this approach, all rate equations are defined by their deviation from a reference state (v⁰, S⁰, E_T⁰ for the reference reaction rate, metabolite concentration, and total enzyme activity, respectively) depending on a logarithmic change in concentration:

(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 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 built for the Calvin cycle, arguably the most important metabolic pathway in plants (for review, see Poolman et al., 2001). 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 I) such as Java Web Simulation 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₂ 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₃ leaves might be suboptimal for maximizing the light-saturated rate of photosynthesis.

The Suc-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 Suc-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 (Saccharum officinarum) operates in the direction of Suc synthesis. In an earlier model of Suc 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 Fru or Glc transporter or the vacuolar Suc 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 sugar cane parenchymal tissue matures in different internodes of the stem.

One of the latest established kinetic models simulates the Asp-derived amino acid pathway in Arabidopsis (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 Asp 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 nonintuitive prediction of the simulation.

Another complex kinetic model in Arabidopsis simulates the synthesis of aliphatic glucosinolates (Knoke et al., 2009), which are important defensive metabolites in plants. The model compared the wild-type Arabidopsis leaves with knockout 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 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 knockout 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 subcellular fractionation, rapid sampling, and noninvasive online measurements. A topic that has to be addressed in the 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 to estimate future developments in plant kinetic modeling. Specifically the metabolic networks in bacteria are well characterized in comparison to plant networks. In an approach for modeling bacterial growth, Fang et al. (2009) 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 that 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 I) and assembly of high-throughput data into new metabolic models, thereby increasing the comparability of the derived models. First steps in the standardization of biological systems from an engineering point of view were recently undertaken for both synthetic biology (Canton et al., 2008) and systems biology (Le Novère et al., 2009).

Ultimately, the success of metabolic modeling (and all other systems biology areas) depends on the interdisciplinary education of students and young scientists in all contributing fields such as mathematics, systems science, statistics, computer science—and last but not least biology.