Enzyme characterisation and kinetic modelling of the pentose phosphate pathway in yeast

We present the quantiﬁcation and kinetic characterisation of the enzymes of the pentose phosphate pathway in Saccharomyces cerevisiae . The data are combined into a mathematical model that describes the dynamics of this system and allows for predicting changes in metabolite concentrations and ﬂuxes in response to perturbations. We use the model to study the response of yeast to a glucose pulse. We then combine the model with an existing glycolysis model to study the eﬀect of oxidative stress on carbohydrate metabolism. The combination of these two models was made possible by the standardised enzyme kinetic experiments carried out in both studies. This work demonstrates the feasibility of constructing larger network models by merging smaller pathway models.


Introduction
The pentose phosphate pathway (PPP) is a central and widely conserved metabolic pathway of carbohydrate metabolism which, in eukaryotic cells, is located in the cytoplasm (see Figure 1). This pathway serves two major functions: production of precursors for biosynthesis of macromolecules and production of reducing equivalents in the form of NADPH. Accordingly, these two roles are reflected in the two major phases of the PPP: in the "oxidative phase", glucose 6-phosphate (G6P) is converted into ribulose 5-phosphate (Ru5P) through the sequential action of glucose-6-phosphate dehydrogenase and 6-phosphogluconate dehydrogenase, with lactonase enhancing the already significant non-enzymic hydrolysis of the 6-phosphogluconolactone (G6L) product of the latter. The "non-oxidative phase" carries out the isomerisation of Ru5P to ribose 5-phosphate (R5P), the epimerisation of Ru5P to xylulose 5-phosphate (X5P) and, through the actions of transketolase and transaldolase, a series of carbon skeleton transfers that can interconvert pentose phosphate into the glycolytic intermediates fructose 6-phosphate (F6P) and glyceraldehyde 3-phosphate (GAP) and erythrose 4-phosphate (E4P). The oxidative branch is considered to be largely irreversible under normal cellular conditions, whilst the non-oxidative branch is reversible [Saggerson, 2009]. The PPP is not a simple linear pathway (see Figure 1) since several carbon atoms are recycled back into glycolysis. Furthermore, the enzyme transketolase catalyses two different reactions in the pathway, resulting in the two substrates of these reactions being competitive inhibitors of each other.
The PPP has three main products: reduced equivalents in the form of NADPH, produced in the oxidative phase, needed in biosynthetic pathways and for maintenance of the oxidative level of cells; R5P, for the biosynthesis of all nucleic acids; and E4P, for biosynthesis of the three aromatic amino acids. Different physiological states require operation of this biochemical network in different modes: in actively growing cells, such as during culture growth in reactors, the pathway must produce a sufficient amount of all three products, since all are required in the construction of new cells. Under stress conditions growth slows down and the only product in considerable demand is NADPH.
Oxidative stress causes damage to all living organisms. They have developed a number of defence and repair mechanisms that are conserved from unicellular to multicellular organisms. Cells typically respond with post-translational modification of a number of proteins, affecting both their localisation and functionality [Godon et al., 1998, Ishii et al., 2007. In particular, oxidative stress in yeast leads to repression of glycolysis and induction of the PPP; this is crucial for maintaining the NADPH/NADP + ratio, which provides the redox power for antioxidant systems [Ralser et al., 2007].
Since the seminal work of [Glock & McLean, 1953], the pentose phosphate pathway has been subjected to a number of quantitative studies, including in yeast [Bruinenberg et al., 1983]. Mathematical models of the pathway have been created in yeast [Vaseghi et al., 1999, Ralser et al., 2007, trypanosome [Kerkhoven et al., 2013], rat [Haut et al., 1974, Sabate et al., 1995 and human [Joshi & Palsson, 1989, Mulquiney & Kuchel, 1999. However, such studies have neglected, or over-simplified, the non-oxidative branch of the pathway. Here we kinetically quantify and characterise various enzymes in the pathway, combine these properties into a mathematical model that describes the dynamic behaviour of this system, and compare the model's predictions to experimental observations of transient metabolite concentrations following a glucose pulse. We go on to examine the response of a combined glycolysis-PPP model to oxidative stress, and compare this to measured metabolite levels.

Kinetics
To determine the kinetic parameters of individual enzymatic reactions of the pentose phosphate pathway, isoenzymes were purified as described previously [Malys et al., 2011]. Spectrophotometric assays were then performed for most of the isoenzymes, following a similar strategy to [Smallbone et al., 2013]. Enzymes were assayed spectrophotometrically through detection of NADPH or NADH, by using coupling reactions where needed, with the exception of ribulose-5-phosphate-3-epimerase (RPE1) and ribose-5-phosphate ketol isomerase (RKI1) which where assayed using circular dichroism (CD, [Kochetov et al., 1978]). Assays were coupled with enzyme(s) in which NAD(P) or NAD(P)H is a product or substrate so that its formation or consumption could be followed spectrophotometrically at 340 nm using an extinction coefficient of 6.62 mM −1 cm −1 , unless the reaction of a particular enzyme consumes or produces NADH or NADPH, in which case no coupling enzymes were needed.
Absorbance measurements were carried out with a BMG Labtech NOVOstar plate reader (Offenburg, Germany) in 384-well format plates in a 60 l reaction volume. All assays were performed in a standardised reaction buffer (100 mM MES, pH 6.5, 100 mM KCl, and 5 mM free magnesium in the form of MgCl 2 ) at 30 • C and were automated so that all reagents in the reaction buffer (including any coupling enzymes) are in 45 µl, the enzyme (to be assayed) in 5 µl and the substrate in 10 µl volumes as described in [Messiha et al., 2011]. For each individual enzyme, the forward and the reverse reactions were assayed whenever possible.
Assays for each individual enzyme were either developed or modified from previously published methodology to be compatible with the conditions of the assay reactions (e.g. pH compatibility or unavailability of commercial substrates). The assay conditions used for each enzyme were as follows: 6-phosphogluconate dehydrogenase GND1 and GND2 were assayed in the reaction buffer in the forward reaction by direct measurement of the production of NADPH as in [He et al., 2007]. The kinetic parameters for each isoenzyme were determined by varying the concentration of each substrate (6-phosphogluconate and NADP) at fixed saturated concentration of the other.
Transaldolase TAL1 and NQM1 were assayed in the reaction buffer in the forward and reverse directions according to [Tsolas & Joris, 1964, Wood, 1972. Since sedoheptulose 7-phosphate was not available commercially, its barium salt was synthesised by Chemos GmbH and converted to the sodium salt just prior to assay according to [Charmantray et al., 2009].
Transketolase TKL1 and TKL2 were assayed for both of their participatory reactions in the reaction buffer in the forward and reverse directions according to [Datta & Racker, 1961, Kochetov, 1982. The kinetic parameters were determined by varying the concentration of each substrate at a fixed saturated concentration of the other for the forward and reverse reactions.
Glucose-6-phosphate dehydrogenase ZWF1 was assayed in the reaction buffer in the forward reaction by direct measurement of the production of NADPH according to [Gould & Goheer, 1976].
Ribose-5-phosphate ketol-isomerase RKI1 was assayed for the forward and reverse reaction by circular dichroism (CD) measurements. The assay was developed based on the fact that ribulose-5-phosphate has a maximum absorbance at 278 nm, with a measured coefficient of -2.88 m • mM −1 mm −1 , and ribose-5-phosphate has an absorbance at 278 nm with a measured coefficient of -0.131 m • mM −1 mm −1 . The data were collected in 400 µl in a 1 mm path length cuvette. In both directions, the CD angle θ can be used to calculate reactant concentrations, from which we can calculate the rate of reaction.
D-ribulose-5-phosphate 3-epimerase RPE1 was assayed for the forward and reverse reaction by CD measurements. The assay was developed and modified from [Karmali et al., 1983]. Xylulose-5-phosphate has an absorbance at 278 nm with a measured coefficient of 0.846 m • mM −1 mm −1 . The signal of CD θ can again be followed to infer the rate of reaction in both directions.
All measurements are based on at least duplicate determination of the reaction rates at each substrate concentration. For each isoenzyme, the initial rates at various substrate concentrations were determined and the data obtained were analysed by the KineticsWizard [Swainston et al., 2010] and COPASI [Hoops et al., 2006] and fitted to Michaelis-Menten kinetics. Data were obtained for all PPP isoenzymes, with the exception of SOL3 and SOL4 which showed no activity (see Table 1). Any missing kinetic parameters were taken from previous models [Vaseghi et al., 1999, Ralser et al., 2007, or given initial estimates using typical values (kcat = 10 s −1 , Km = 0.1 mM, [Bar-Even et al., 2011).

Proteomics
We attempted to measure the absolute quantities of all isoenzymes in this pathway through the QConCAT technology [Benyon et al., 2005]. Total cell protein was extracted from turbidostat yeast cultures as described earlier . Data analyses were performed using the PrideWizard software [Swainston et al., 2011] (see Table 2). Concentrations were then calculated from copy number using a typical cytoplasmic volume of 5 fL [Smallbone et al., 2013].
Only four of the isoenzymes (Gnd1, Sol3, Tal1 and Tkl1) were detected using this approach. In the case of Gnd1/Gnd2 and Tal1/Nqm1, only the most abundant isoenzyme was detected in each case, and it is likely that the expression level the less abundant isoenzyme was not necessarily zero but at least it was below the detection limit. The remaining three enzymes (Rki1, Rpe1 and Zwf1) were detected at higher levels in a previous study ([Ghaemmaghami et al., 2003], detailed in Table 2). Moreover, these are soluble cytoplasmic proteins, so we can assume they were likely present in the extracted protein preparations (rather than sequestered to membranes, and subsequently lost as insoluble material). There are two possible explanations for the failure to detect these proteins: poor or incomplete proteolysis (trypsin miscleavage) or unexpected post-translational modifications, either naturally occuring or inadvertently introduced during the experimental protocol.
We note that copy numbers of PPP enzymes measured here were twenty-fold higher, on average, than in a previous study [Ghaemmaghami et al., 2003]. Thus using the data from Ghaemmaghami et al. directly to fill in for the missing measurements is not appropriate here. Rather, in cases where one of two isoenzymes was not quantified (Gnd2, Nqm1), the same ratio was maintained as as in [Ghaemmaghami et al., 2003] (i.e. we use the same proportions of the two isoenzymes). For the remaining three undetected enzymes (Rki1, Rpe1, Zwf1) the value reported in that study was multiplied by twenty to provide an initial estimate.

Model construction
From a modelling perspective, the enzyme kinetic constants and protein concentrations represent the parameters of the system, while the metabolite concentrations (Table 3) represent the variables. Combining the protein concentration data with those for the enzyme kinetic parameters, together with the measured steady state metabolite levels, allows a mathematical model to be produced for this system (Table 4) in ordinary differential equation format. This model considers, in the first instance, the PPP in isolation. Thus we consider three boundary metabolites to be fixed: F6P, G6P and GAP.
To consider oxidative stress, however, we expanded the model to combine it with our recently published model of glycolysis (that includes trehalose and glycerol metabolism) [Smallbone et al., 2013], where the enzymatic parameters were determined in the same condiditions as described here. This combined glycolysis:PPP model contains 34 reactions, and allows calculation of the concentration of 32 metabolites (variables). Importantly, it allows us to compare the joint response of both pathways to environmental perturbations.

Glucose pulse
We used data published earlier for a perturbation in the G6P concentration following a glucose pulse [Vaseghi et al., 1999] with the following format G6P = 0.9 + 44.1 t 48.0 + t + 0.45 t 2 and compared our model's predictions to the experimental observations in [Vaseghi et al., 1999] with respect to NADPH and P6G concentrations (see Figure 2). Whilst the present model contains many parameters that we measured under standardised conditions, a few parameters were not possible to determine experimentally and were therefore obtained from the literature. The relative contribution of each of these parameters to the quality of fit to data may be ranked using sensitivity analysis. If we were unable to closely match the data using only the literature value of the most important parameter, we tried using two parameters, and continued adding parameters until a satisfactory match was obtained. Only five parameters had to be obtained from the literature in this way (see Table 5) to provide the good match seen in Figure 2. Of these five, three were initial guesses, one (ZWF:Kg6l) was measured under other conditions, and only one ([Gnd1]) had been measured by us (see above), but nonetheless fitted to the data.

Oxidative stress
One of the proteins that responds to oxidative stress is the glycolytic enzyme glyceraldehyde-3phosphate dehydrogenase (TDH). In response to high oxidant levels this enzyme is inactivated and accumulates in the nucleus of the cell in several organisms and cell types [Chuang et al., 2005, Shenton & Grant, 2003]. Thus, we simulate in silico oxidative stress through reduction of TDH activity in the combined glycolysis:PPP model to 25% of its wild-type value, following the approach of [Ralser et al., 2007]. Cells also respond to the presence of oxidative agents through slower growth, which we translate in our model as reducing the requirement for E4P and R5P (the biomass precursors); we thus reduce the rate of consumption of these by two orders of magnitude from their reference values. The defence against the oxidant agent requires reductive power which is ultimately supplied by NADPH (e.g. through glutathione); we thus also increase the rate of NADPH consumption by two orders of magnitude. We may then compare predicted changes in metabolite concentrations to those measured in response to H 2 O 2 treatment [Ralser et al., 2007], a typical oxidative stress agent [Godon et al., 1998].
The results of these simulations are presented in Table 6. They show that seven of the eight qualitative changes in metabolite concentrations are correctly predicted by the model. A difference between the experimental data and the predictions was only observed for the metabolite glycerol 3-phosphate (G3P), where the simulation predicts a small increase, but experimentally we observe a small decrease.
As the qualitative predictions reasonably matched the experimental data set, we moved on to calculate the influence of oxidative stress on carbon flux. We found that the ratio of fluxes into glycolysis (via PGI) and into PPP (via ZWF) was 18:1 under reference conditions, but this value reduced to 9:1 under the oxidative stress conditions, corroborating the hypothesis that oxidative stress leads to a redirection of the carbohydrate flux [Ralser et al., 2007]. The ratios are consistent with experimental measurement showing that, in aerobic growth conditions on glucose minimal medium, PPP activity accounts for some 10% of the total consumption of glucose [Blank et al., 2005].

Control analysis
Metabolic control analysis (MCA) is a biochemical formalism, defining how variables, such as fluxes and concentrations, depend on network parameters. It stems from the work of [Kacser & Burns, 1973] and, independently, [Heinrich & Rapoport, 1974]. In Table 7 (a), we present the scaled flux control coefficients for the (fitted) PPP model. These are measures of how a relative change in enzyme activity leads to a change in steady state flux through the system. For example, from the third row of the table, we predict that a 1% increase in GND levels would lead to a 0.153% decrease in RPE flux.
The table shows us that flux into the pathway (via ZWF) is entirely controlled by ZWF, SOL, GND and NADPH oxidase (the latter representing all processes that oxidise NADPH). Returning to Figure 1, we see that these correspond to the first three steps of the pathway plus NADPH recycling -the oxidative phase. The table also shows the overall control of each step of the pathway, taken in COPASI [Hoops et al., 2006] to be the norm of the control coefficients. We see that little control is exerted by the RPE and TKL (R5P:S7P) steps. The three sinks have high overall control, and as such we would expect fluxes through the pathway to be highly dependent on growth rate and stress levels.
In the oxidative stress simulation the control distribution changes, as presented in Table 7 (b). The main observation from these data is that the control of the pathway input flux is now much lower by the NADPH oxidase -this is somewhat expected since the rate of this step increased 100× and thus became less limiting. Less intuitive is the reduction of overall control of the network by RKI (the reaction that produces ribose 5-phosphate). However this result implies that under oxidative stress the PPP is essentially insensitive to the "pull" from ribose use for nucleic acid synthesis, which agrees with the observation that growth is arrested under these conditions.

Discussion
The pentose phosphate pathway, depicted in Figure 1, is a central pathway in yeast and in most organisms and serves two main functions: maintenance of the NADPH:NADP + ratio, and production of several precursors for biosynthesis of macromolecules. These two roles of the pathway are mirrored in its structure and it consists of two semi-independent parts; the oxidative branch reduces NADP + , whilst the non-oxidative branch creates R5P, a precursor for nucleic acid biosynthesis, or E4P, a precursor for aromatic amino acids and some vitamins. The PPP is intimately connected with glycolysis as it diverts some of its flux away from energy production. Furthermore, the two pathways have three metabolites in common: G6P, F6P and GAP.
In order to describe a biological system such as PPP quantitatively, the kinetic properties of all its components need to be established in conditions close to the physiological [van Eunen et al., 2010. Where possible, they should represent a system in steady state, where all measurements, even if carried out at different times, are performed under identical conditions. Following the methodology previously applied to glycolysis [Smallbone et al., 2013], robust and standardised enzyme kinetics and quantitative proteomics measurements were applied to the enzymes of the pentose phosphate pathway in the S. cerevisiae strain YDL227C. The resulting data are integrated in a kinetic model of the pathway. This is in contrast to previous studies [Vaseghi et al., 1999, Ralser et al., 2007, where kinetic parameters were taken from various sources the literature: "The kinetic constants were determined using enzymes from five different species (human, cow, rabbit, yeast, E. coli) in different laboratories over a period of more than three decades." [Ralser et al., 2007] We use our standardised model to study the response of the pentose phosphate pathway to a glucose pulse, finding a good agreement between measured and predicted NADPH and P6G profiles ( Figure 2). We go on to use model to study the combined response of glycolysis and PPP to oxidative stress, and find that a considerable amount of flux is rerouted through the PPP.
Our modelling approach also reveals a discrepancy between the observed change in G3P levels following stress cannot be predicted by current understanding of glycolysis and PPP; following the "cycle of knowledge" [Kell, 2006], it is of interest to direct future focus towards glycerol metabolism in order to improve the accuracy of this model.
It is important to highlight that we were not able here to quantify the concentration of all enzymes in the pathway, thus having to rely on crude estimates. It is quite possible that the physiological conditions under which the cells were measured by [Ghaemmaghami et al., 2003] were very different than those used here, which could result in inaccurate estimates for the concentration of several enzymes. However the fact that we have measured k cat parameters for those enzymes will allow easy correction of the model if accurate enzyme concentrations are determined later. Indeed, these data will allow to account for changes in enzyme concentrations resulting from a longer term response of the cells, through protein degradation or increased protein synthesis rate due to changes at the level of transcription and translation.
The combined PPP and glycolysis model demonstrates the value of standardised enzyme kinetic measurements -models thus parameterised can be combined to expand their scope, eventually forming large-scale models of metabolism [Snoep, 2005, Snoep et al., 2006. Indeed the combined glycolysis:PPP model could be expanded to consider enzyme concentrations as variables (through accounting for their synthesis and degradation, reflecting gene expression and signalling) which would improve its utility in predicting a broader array of conditions. Such an expansion of models to cover wider areas of metabolism and cellular biochemistry will lead to digital organisms, as shown in a recent proof of principle for the simple bacterium Mycoplasma genitalium [Karr et al., 2012].
The "bottom-up" strategy used here is to combine compatible kinetic models (PPP and glycolysis) exapnding them towards a larger metabolic model. An alternative ("top-down") strategy is to start with a large structural yeast network [Herrgård et al., 2008, Dobson et al., 2010, Heavner et al., 2012, Heavner et al., 2013, Aung et al., 2013, then add estimated kinetic parameters and, through successive rounds of improvement, incorporate measured parameters [Smallbone et al., 2010, Stanford et al., 2013, in an automated manner where possible , Büchel et al., 2013. Can these two strategies be combined into a more robust and scalable approach?
In summary, we present here a model of the yeast pentose phosphate pathway that we believe is the most realistic so far, including experimentally determined kinetic parameters for its enzymes and (some) physiological enzyme concentrations. A more complex model resulting from the combination of this PPP model with a previous glycolytic model [Smallbone et al., 2013] was possible due to the standardized way in which the kinetic parameters were measured. This opens up the prospect of expanding models to eventually cover the entire metabolism of a cell in a way that makes them compatible with a further improvement, by including the effects of changes in gene expression.