VU Research Portal

The evolution of metabolic strategies

Wortel, M.T.

2015

document version

Publisher's PDF, also known as Version of record

Link to publication in VU Research Portal

citation for published version (APA)

Wortel, M. T. (2015). The evolution of metabolic strategies.

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal ? Take down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

E-mail address:

7

|

Lost in transition

startup of glycolysis yields subpopulations of non-growing cells

In collaboration with: Johan H. van Heerden, Frank J. Bruggeman, Joseph J. Heijnen, Yves J.M. Bollen, Robert Planqué, Josephus Hulshof, Tom G. O’Toole, S. Aljoscha Wahl and Bas Teusink

Abstract

Cells need to adapt to dynamic environments. Yeast that fail to cope with dynamic changes in the abundance of glucose can undergo growth arrest. We show that this failure is caused by imbalanced reactions in glycolysis, the essential pathway in energy metabolism in most organ-isms. The imbalance arises largely from the fundamental design of glycolysis, making this state of glycolysis a generic risk. Cells with unbalanced glycolysis co-existed with vital cells. Sponta-neous, non-genetic metabolic variability among individual cells determines which state is reached and consequently which cells survive. Transient ATP hydrolysis through futile cycling reduces the probability of reaching the imbalanced state. Our results reveal dynamic behavior of glycolysis and indicate that cell fate can be determined by heterogeneity purely at the metabolic level.

Main text

hexoki-nases (HXK1 and HXK2), with respect to glucose (Hohmann et al, 1996). This negative feedback loop was hypothesized to slow down the upper part of glycolysis and restore the bal-ance in wild-type cells (Teusink et al, 1998), but T6P-insensitive hexokinase mutants do grow on glucose (Bonini et al, 2003). An alternative hypothesis assumes a reduced activity of glycer-aldehyde 3-phosphate dehydrogenase (gapdh) because of the low amounts of its substrate, cytosolic phosphate ( P_i), in Tps1∆mutants. Accordingly, P_irelease should enhance gapdh ac-tivity and restore balance. In line with this hypothesis, enhanced glycerol production, which releases P_i, restores growth of tps1∆mutants on glucose (Luyten et al, 1995). Trehalose pro-duction from glucose 6-phosphate (G6P) also releases P_i (Fig. 1A); however, the capacity of trehalose synthesis was believed to be too low to provide enough Pi for the high flux of

glycolysis (François and Parrou, 2001).

We used a computational approach to better understand the complex phenotype of tps1∆

mutants. We adapted an existing kinetic model of glycolysis (Teusink et al, 2000) by (i) introduc-ing P_i as an explicit variable in the model (rather than as a fixed “commodity” metabolite); and (ii) allowing for the mobilization of Pifrom vacuolar stores, based on in vivo NMR data that

de-scribe this behavior (Van Aelst et al, 1993). The latter was necessary, as the net accumulation of phosphate-containing glycolytic intermediates that is observed experimentally in the imbal-anced state (Fig. 7.1B), is not possible without the import of P_i(details on these and other small adjustments to the original model are provided in supplementary online text).

Two glycolytic states co-exist

Simulations representing the tps1∆-mutant resulted in dynamic metabolite profiles that were qualitatively similar to experimentally observed profiles (Fig. 7.1B), i.e. all metabolites were bal-anced, except for the intermediates between the upper and lower parts of glycolysis (Fig. 7.1B, S3). Known experimental rescue mechanisms for the tps1∆-mutant, such as reduced activity of hexokinase (Hohmann et al, 1996) or enhanced glycerol production (Luyten et al, 1995) could be reproduced in silico (Fig. S11C, S14C and D).

Our tps1∆-mutant model could reach another state, which resembled the wild-type steady state with proper flux, high ATP and Pilevels, and normal FBP levels. Whether this state was

reached depended on the initial concentrations of the metabolites, as shown for FBP and P_i (Fig. 7.1B). Hence two stable outcomes (states) co-existed in the model, a global steady state and the imbalanced state: the latter is a non-typical stable state as some variables are not con-stant in time but rather accumulate. Other systems with more than one stable state, as found in sporulation (Veening et al, 2005) or differentiation (Wang et al, 2009), often result in phenotypi-cally different subpopulations in an isogenic population. We assessed experimentally whether we could find evidence for such subpopulations as well, by asserting that only the functional glycolytic state would support growth. Based on serial-dilution plating of wild-type and tps1∆

FBP F6P transport hk glk ⁄ ATP tpi glt poly(Pi ) Vacuole Pi exchange Membrane transport Enzymatic reaction Allosteric activation Allosteric inhibition Glc_ext Glc_int 2PGA PEP PYR ACE Succinate Ethanol NAD NADH ADP NADH NAD Pi gapdh ADP ATPase Pi CO₂ 2 pyk GAP BPG 3PGA DHAP G3P Glycerol gpd Pi P_i P_i NADH 3 NADH 3 NAD NAD Upp er Gl yc ol ysis Lo w er Gl yc ol ysis pgk pgm eno adh pdc ATP ATP ald G6P pgi pfk ADP ATP ADP ATP ADP TREH T6P 2 AMP tps2 tps1 Glycogen A C D E tps1Δ

2% Gal 2% Glu 2% Gal 2% Glu WT Dilution 100 10-1 10-2 10-3 10-4 10-5 0 5 10 15 20 25 30 35 40 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 time (hrs) OD600 1 2 3 4 5 0 10 20 FBP (mM) P i (mM) 0 2 4 6 10 20 30 time (min) mM 0 2 4 6 10 time (min) mM FBP ATP Pi FBP ATP Pi Imbalanced state Steady state B (i) (ii) Glycerol stock Glu Gal Glu Gal 1/10

dils glucose 1/10dils tolerant colony single colony wash CBSgal Buff er Buffer grow to midlog 2% 2% Gal 2% tps1Δ + 2% Gal tps1Δ + 2% Glu WT + 2% Glu WT + 2% Gal Dilution 100 10-1 10-2 10-3 10-4 10-5 100 10-1 10-2 10-3 10-4 10-5 Dilution G al Glu Gal Glu tps1Δ WT G al Glu Gal Glu tps1Δ WT 0 20 40 60 0 10 20 30 0 1 2 3 tmin mmol L Cytosol 1 20 40 60 0 20 40 FBP ATP Pi 0 1 2 3 FBP (mM) ATP , P i (mM) time (min)

Figure 7.1. The co-existence of two glycolytic states underlies glucose tolerant tps1∆subpopulations. (A)

glucose plates and subjected to another transition via galactose to glucose, the original subpop-ulation structure with less than 1 in 103 glucose-tolerant colonies was restored (Fig. 1E and Fig. S5C). This argues against a genetic basis. We therefore conclude that there is a small sub-population of glucose-positive tps1∆ cells that arises from spontaneous phenotypic—as opposed to genetic—variability.

We re-examined the reported inhibitory effect of low concentrations of glucose on tps1∆

mutant growth in the presence of excess galactose (Neves et al, 1995). Plating experiments again showed that the growth inhibition observed at the population level is in fact caused by a glucose-dependent increase in the size of a subpopulation that was unable to grow (Fig. S6A). A simple mathematical model of population growth dynamics with different growing and non-growing subpopulation sizes (supplementary online text) reproduced the experimental data (Neves et al, 1995) very well (Fig. S6B).

Intracellular pH reveals two metabolic subpopulations

To visualize the two subpopulations, we made use of the observation that tps1∆mutants, when exposed to glucose, are unable to maintain pH homeostasis because they produce too little ATP (Van Aelst et al, 1993). Hence, after glucose addition, the intracellular pH (pH_i) of tps1∆

populations decreased by more than 1 pH unit compared to that of wild-type cells (Fig. 7.2A). After confirming that different subpopulations could be distinguished on the basis of pHi signals (Fig. 7.2B) we used flow cytometry as a high-throughput approach to study the structure of both tps1∆and wild-type populations exposed to galactose and glucose. We found two subpopula-tions of tps1∆ cells of sizes that agreed with the plating assays and growth lag phases (Fig. 7.2C). We also tested wild-type cells, because analysis of the wild-type version of the model also showed two stable states (Fig. 7.9). Indeed, a subpopulation of wild-type cells with low pH appeared in cultures exposed to glucose, but a similar response was observed for galactose (Fig. 7.2C). The size of the subpopulation, about 7%, was larger than the model sug-gested (but not unexpectedly, Fig. 7.9). These results indicate that the imbalanced state is a general property of glycolysis that cannot be fully prevented by regulatory mechanisms operative in wild-type cells.

Metabolic variability determines state of glycolysis

Phe-A C B 0 1 Ratio (405nm/488nm) Lower pH 0.25 0.5 0.75 Higher pH 2% Gal 2% Glu 0 5 10 15 20 25 5.5 6 6.5 7 7.5 8 time (min) pH i tps1Δ 2% Glu 2% Gal WT 2% Glu 2% Gal 0 0.5 1 1.5 2 2.5 0 5 10 15 Gal Glu Pre-perturb 7.3 % tps1Δ _{Gal + 2.5mM Glu} 0 0.01 0.02 0.03 0.04 0.05 0.4 % 1 1.5 2 2.5 0 0.050.1 0.150.2 0.25 3.4 % Normalized percentage Ratio (405/488 nm) 0 2 4 6 8 Gal Glu Pre-perturb 7.0 % 8.4 % WT 2% Gal+ 2 mM Glu

Figure 7.2. Intracellular pH reveals distinct metabolic subpopulations. (A) Population-level pHiresponses show

the disruption of pH homeostasis of tps1∆cells in response to 2% glucose. Following a 2% glucose or galactose pulse, the pHi of WT and tps1∆populations are shown in time. (B) Fluorescent microscopy shows that distinct

metabolic states can be visualized using pHi readouts. (C) Flow cytometery measurements based on pHluorin

signals, reveal the presence of distinct subpopulations in both WT (top graph) and tps1∆populations (bottom graph) following perturbations with glucose and galactose (black, pre-perturbation sample in wash buffer; blue, 2% Gal; red, 2% Glu; orange, 2% Gal + 2.5mM Glu).

notypic variation by non-genetic variability is usually studied in genetic circuits that naturally operate at a stochastic regime with low copy numbers for key components (Raj and Van Oude-naarden , 2008). It is relevant to compare the above observations with frameworks such as fluctuation-induced bistable switching (Balaban et al, 2004; Levy et al, 2012) which has previ-ously been linked to the emergence of phenotypic heterogeneity. In contrast to such stochastic-switching phenomena, the emergence of the two distinct phenotypes (viable and non-viable) described here does not depend on the co-existence of two qualitatively different physiological states prior to a glucose perturbation. Our interpretation of the above data is that spontaneous, non-genetic variation between cells creates a continuous probability distribution for metabolite concentrations and metabolic fluxes. Within the space of these initial physiologi-cal states, a subspace exists that characterizes the cells that survive a sudden glucose excess exposure (Fig. S7) (also see supplementary online text, section 3.4).

To assess which parameters and initial conditions most affect the probability to reach an imbalanced or functional steady state, we performed a linear discriminant analysis over all our model simulations (supplementary online text). A single discriminant accounted for 99% of the differences in initial conditions that lead to either the balanced or imbalanced states (Fig. S11A). This discriminant identified parameters and initial metabolite concentrations that either tend to reduce the flux through the upper part of glycolysis or enhance the flux through lower glycolysis (Fig. 7.3A). The parameters related to the primary mechanisms known to rescue the tps1∆

Upper Glycolysis Lower Glycolysis Gly cer ol Br anch LD2 LD1 1.0 0.5 0.0 -0.5 -1.0 1.0 0.5 0.0 -0.5 -1.0 Vmax glt V_max hk V_max pfk Vmax g3pdh Vmax gapdhfwd Pi time (hrs) OD 600 0 4.3 10 21.5 30 43 50 75 0 1 2 3 4 5 6 7 8 Ethanol (mM) Viable subpopulation %

Random Sampling Results Inferred from growth curves

B C A 5 10 15 20 25 30 35 40 0 0.2 0.4 0.6 0.8 1 2% Glu only 2% Glu + 4.3 mM EtOH 2% Glu + 21.5 mM EtOH 2% Glu + 43 mM EtOH

Figure 7.3. Metabolic subpopulations are caused by small variation in metabolic variables, and their sizes can be manipulated. (A) Linear discriminant analysis of randomly sampled initial conditions (metabolites and Vmaxes)

Size of subpopulations can be manipulated

We looked for ways to experimentally influence the size of the two subpopulations. Respiratory inhibitors, in particular antimycin A (Blázquez and Gancedo, 1995), have been shown to improve growth in the presence of glucose. However, ethanol, the solvent of antimycin A, pro-duced pronounced decreases in the lag phase (Fig. 7.3B), completely dominating the antimycin A effect (Fig. S5D). These results were confirmed by plating assays (Fig. S6C), and flow cytom-etry measurements of pHi (Fig. S9) (supplementary online text). We repeated the in silico random-sampling approach at different ethanol concentrations and reproduced the positive effect of ethanol (Fig. 7.3C). The model showed that the increase in ethanol concentration in-creased Pirelease through glycerol formation (Fig. S11B) driven by an increased NADH/NAD

ratio (Overkamp et al, 2002). This model prediction was experimentally tested in tps1∆mutant cultures by addition of formate. Although formate cannot be used as a carbon source by yeast, its conversion to CO2 by formate dehydrogenase (Geertman et al, 2006) enhances NADH for-mation. Indeed, formate additions similarly decreased the lag phase, implying a NADH-driven increase in glycerol production that subsequently increases the proportion of tps1∆cells in the population that could grow on glucose (Fig. S5E).

Core model explains and generalizes dynamics

To generalize our findings and to provide a deeper understanding of the observed co-existence of two stable states, we captured the essential features of the large model in a reduced model. This generalized core model (supplementary online text) only considers the concentrations of FBP, ATP and P_iand four reactions: a lumped upper glycolysis reaction (vupper), a lumped lower part

of glycolysis (v_lower), an ATP-demand reaction (v_ATP), and a P_i import or export reaction (v_Pi) (Fig. 7.4). Detailed mathematical analysis showed that such a generalized glycolytic pathway has two stable states representing a functional steady state and an imbalanced state. Figure 4 shows the system dynamics that lead to these two states: the difference between the left and right panel is only the initial P_i level (10.4 and 9.4 mM, respectively). How can the different outcomes in these simulations be explained? At the start of the simulation (when glucose is added), vupper > vlower, and this difference causes the concentration of the intermediate FBP to

increase (asdFBP_dt = v_upperv_lower). For a balanced steady state, v_lowerneeds to accelerate to equal the rate of vupper (note that this challenge becomes bigger if the activity of vupper is higher). Piis

the phosphate source for FBP and therefore FBP accumulation results in a drop in P_i(Fig. 7.4). As both FBP and P_i are substrates of v_lower the accumulation of FBP stimulates v_lower, while the drop in Pi tends to slow it down: Which effect is dominant may determine the fate of the

system. If P_iis high initially, a drop in P_iwill not affect v_lower and the FBP increase will dominate, resulting in the functional steady state being reached (Fig. 7.4, left panel). Similarly, if Pi is

liberated quickly enough directly by storage/uptake (v_Pi) or indirectly by ATP hydrolysis (v_ATP), P_i will drop less quickly and the balanced steady state can also be reached. If, however, P_iis low at the onset of glucose addition, or Pimobilization is too slow, or both, the decrease in Piwill

system will collapse to the imbalanced state (right panel Fig. 7.4). Once in this imbalanced state, continuous Pimobilization from uptake or storage paradoxically maintains the imbalance. In this

low P_i, low ATP state, imported P_ienhances the rate of v_lower, but the concomitant production of ATP will increase v_upper two times more (due to the stoichiometric coupling of ATP in glycolysis, Fig. 4B). Hence, the imbalance and thus FBP accumulation will only get bigger with faster Pi

import, as observed in the core model and the detailed model (Fig. 7.6).

Trehalose metabolism constitutes transient futile cycling

The core model predicted that enhanced Pimobilization through ATP hydrolysis by vATP could

enlarge the set of initial conditions that leads to a functional steady state, or could even result in the disappearance of the imbalanced state altogether (Fig. 7.12). This was confirmed in the de-tailed model. We realized that, in yeast, the full trehalose cycle (Fig. 7.1A) could act as a mechanism for ATP hydrolysis through futile cycling. This cycling of trehalose should be able to remove the existence of the imbalanced state or at least reduce the probability to reach it, pro-viding a rationale why trehalose metabolism affects glycolytic function. To estimate the dynamic fluxes through the trehalose network upon a transition to glucose excess, we used a dynamic [13C]-labeling approach (supplementary online text). Wild-type cells were grown in a glucose-limited chemostat and treated with either 110 mM [12C]-glucose or uniformly labeled [U-13C6]-glucose pulses. The time course of the concentrations and the average carbon label-ing enrichments for key intermediates are shown in Fig. 7.5A. For most metabolites up to full enrichment was achieved very rapidly; in contrast, the large trehalose pool was enriched to only 14% at the end of the experiment. From these data the flux of glucose through the different gly-colytic enzymes was estimated based on a hybrid modeling approach (Abate et al, 2012). The flux profiles indicated that: (i) fluxes changed rapidly after a glucose pulse, at similar time scales as key metabolites, such as ATP (Fig. 7.5B), (ii) fluxes through TPS1 and TPS2 increased and subsequently decreased between 0 and 5 min and (iii) at its maximum, as much as 28% of the glucose taken up was branched into trehalose (Fig. 5C). The transient nature of the flux through the trehalose pathway is consistent with the existence of two stable states in glycolysis, i.e. once the system has reached the viable steady state, the need for excessive ATP hydrolysis through futile trehalose cycling has disappeared. Thus, we suggest that trehalose cycling con-stitutes a transient futile cycle, large enough to push the system’s dynamics into the functional steady state.

Different mechanisms contribute to robustness

Finally, we examined the contribution of various aspects of trehalose metabolism to establish proper glycolytic functioning, through random sampling of initial conditions in the full kinetic model (Fig. 5D). Whereas the combined trehalose cycling and negative (trehalose 6-phosphate mediated) feedback on hexokinase resulted in 100% viability in the wild-type ver-sion of the model, apparent futile cycling of trehalose alone resulted in a regular steady state in 76% of the sampled cases. Removal of G6P without Pirelease (only possible in silico) resulted

tre-A

B

C

D

E

F

Steady State

Imbalanced State

GLUCOSE vupper vlower vATP 2ATP polyP 2ETOH FBP vP 2P_i 4ATP GLUCOSE vupper vlower vATP 2ATP polyP 2ETOH vP 2Pi 4ATP FBP x₂ 0 10 20 30 5 10 15 20 time (min) mM vupper - vlower = dFBP/dt = 0 mM.min -1 time (min) 0 5 10 15 20 25 30 -1 0 1 2 3 vP = 0 0 5 10 15 20 25 30 -1 0 1 2 3 mM.min -1 time (min) vupper - vlower = dFBP/dt > 0 vP > 0 FBP ATP Pi FBP ATP Pi time (min) mM 0 10 20 30 5 10 15 20

Figure 7.4. Generalized core model of glycolysis can reach two stable, co-existing, states. The left panel shows

the global steady state, the right panel the imbalanced state. The difference between panels is the initial Pilevel (10.4

and 9.4, respectively). (A and B) Stoichiometry of the core model, with red arrows emphasizing the vacuolar flow of Pi

from polyphosphates (polyP). The coupling between the upper and lower part of glycolysis through ATP is emphasized by the red dashed line (B). (C and D) Metabolite levels for a simulation of the core model, resulting in steady state (metabolite levels constant in time, C) or imbalance (FBP accumulation at very low Piand ATP levels, D). (E and F)

characteristic rates that specify the states: the red dashed lines indicate the difference in rate between upper and lower glycolysis (vuppervlower), which is zero at steady state (E) and is positive at the imbalanced state (F). The dashed

blue lines represent the vacuolar import rate of phos (vPi), which should be zero at steady state. In Fig. 4F, the

constant positive vPiindicates mobilization of Pi, which sustains accumulation of FBP (red dashed line) through the

100 200 300 0 1 2 3 4 5 ATP (µmol .gDW -1 )

time (s)

C

(i)

100 %

ATP Pi Trehalose Glci G6P

(ii)

76 %

ATP Pi Trehalose Glci G6P

(iii)

4 %

Trehalose Glci G6P

(iv)

0.06 %

Glci G6P

TREHALOSE CYCLE

UDPG Glucose v5 v6 0 100 200 300 0 0.05 0.1 0.15 0.2 0 tps2 0 100 200 300 0 20 10 0 00.2 0.4 0.6 0.8 1 0 100 200 300 0 0.05 0.1 0.15 0.2 0 tps1 0 100 200 300 0 200 100 0 00.2 0.4 0.6 0.8 1 0 100 200 300 0 0.05 0.1 0 Trehalase 0 100 200 300 0 0.8 0.4 0 0 100 200 300 0 4 2 0 00.2 0.4 0.6 0.8 1 T6P TREH 20 10 0 0 100 200 300 0 0 0.2 0.4 0.6 0.8 1 G6P F6P 0 100 200 300 0 0.5 1 Glu uptake/hxk 1.5 2 G1P v1 v7 v8 v9 v12 v2 pgi pgm upg v3 _pfk g6pdh A 100 200 300 0 0 0.1 0.2 0.3 0.4 _tps1/hxk (v7/v1)

time (s)

Flux r

atio

B D 100 200 300 0 0 0.05 0.1 0.15 0.2 tps2/hxk (v8/v1)

Figure 7.5. 13C tracer enrichment reveals highly dynamic flux distributions through the trehalose cycle. (A)

A selection of data [complete dataset in supplementary online material in Science Online], superimposed on the tre-halose cycle. White boxes contain metabolite data, with concentrations (µmol·gDW−1) in blue and tracer enrichment (fraction) in red (symbols represent measurements, lines represent model fits); x-axes show time in seconds. Grey boxes show flux profiles (µmol·gDW−1 ·s−1) of the indicated reactions in time (s). (B) ATP concentration profile show a response time similar to that of flux channeling toward the trehalose pool via the tps1 reaction. (C) Dynamic flux ratios of tps1 (v7) and tps2 (v8) relative to hxk (v1) show that up to 28% of glucose is dynamically routed into

halose synthesis by P_i (Londesborough and Vuorio, 1993) reinforces this picture as low P_i concentrations will relieve inhibition and ensure Pimobilization. These results provide a strong

basis for the interpretation of population-level phenotypes of various mutants in trehalose metabolism and hexokinase (supplementary online text).

Startup of glycolysis requires dynamic regulation

Intriguingly, PKM2 expression in tumors coincides with an alternative phosphoglycerate mutase activity that instead of ATP produces Pifrom PEP (Vander Heiden et al, 2010). Such an activity

would fit with the role of P_i release to ensure robust functioning of glycolysis. Interfering with protective mechanisms against the metabolic imbalance state in tumor cells, or perhaps, counter-intuitively, enhancing glucose uptake rather than inhibiting it, may provide rationales for sophisticated treatment strategies.

Materials and Methods

Kinetic modelling

All models were implemented and analyzed using Mathematica 9.0 (Wolfram Research). The time simulations were performed with the NDSolve function. A steady state is defined as a state, characterized by the metabolite concentrations, where all time derivatives of internal metabolite concentrations are equal to zero. Steady states were calculated by solving these equalities with the FindRoot function. Metabolite concentrations from the time simulations, after 250 simulation minutes, were used as initial values for the steady state estimations with the FindRoot function.

The detailed kinetic model and the core model described in this paper are available as sup-plementary files in SBML format. In addition, an interactive web application is provided at: http://www.ibi.vu.nl/sysbio/tpsmodel/, and demonstrates the effects of several parameters on the detailed kinetic model.

Modelling metabolic heterogeneity

Model variations and random sampling descriptions

The kinetic model, as detailed in supporting online, was used to explore the relationship be-tween initial conditions and the probability of obtaining a regular steady state with a model representing the tps1∆state (see supporting materials, section 2). We took enzyme expression levels and metabolite concentrations as two sources of inter-individual metabolic variation in a population. In order to simulate this heterogeneity, we assumed Gaussian distributions for both enzyme expression levels (represented by vmax values) and initial metabolite concentrations.

Variance was set such that the probability of a sampled value deviating more than 20% from the reference value is less than 1 in 103; this equals a coefficient of variation (i.e. a mean-normalized standard deviation) of 6.1% (3.29 standard deviations = mean ± 20%, for Gaussian distributions). To additionally evaluate the effect of ethanol on the probability of reach-ing a steady state, we performed samplreach-ings at varyreach-ing ethanol concentrations: 0, 4.3, 10, 21.5, 30, 43, 50, 62.5, 75, 100, 200 and 500 mM.

We randomly drew more than 106 unique Vmax and initial metabolite concentration sets, at

would arise when no free phosphate is available to drive FBP accumulation. Based on the evaluation outcome each randomly drawn data set was categorized and saved.

Random sampling, time simulations and steady state evaluations were performed, as above, using Mathematica 9.0 (Wolfram Research).

Discriminant Analysis

Using the output from the random sampling evaluations, 5000 samples were drawn (randomly) from the saved data sets, for both imbalanced and regular steady state solutions, at each ethanol concentration, yielding a data set with 28 variables (14 initial metabolite and 14 v_max values) and 24 independent classes (2 groups: imbalanced vs. steady state, at 12 different ethanol concentrations, see above). The discriminant analysis was performed with the lda function of the MASS package in the R (version 2.14.2) statistical environment (Ihaka and Gentleman, 1996).

Supplementary material section 2:

Kinetic model of Glycolysis

Detailed description of changes to the original model

The detailed kinetic glycolysis model from Teusink et al. (Teusink et al, 2000) was used and updated by introducing minor adjustments relevant for the current study. Firstly, we introduced a proxy for the inhibition of hexokinase by T6P. Since most of the study deals with the mutant case where this inhibition is actually absent, we only incorporated this feedback loop in a coarse-grained fashion, to have a representation of wild type to benchmark against. Since T6P is not explicitly defined in the model and the magnitude of change in the G6P pool is similar to that of T6P dynamics (Hohmann et al, 1996), we used the G6P levels as a proxy to simulate the inhibition of hk by T6P (Eqn (7.1)). Several different values are reported for the T6P inhibition constant of hk (Hohmann et al, 1996). We decided to use the average of 0.04 mM and 0.1 mM (see (Hohmann et al, 1996) for details), i.e. 0.07 mM (Tabel 7.1). The rate equation for hk becomes: v_hk = vmax ,hk  ATP· GLCi kATP· k_GLCi − ADP· G6P keq,hk· kATP· k_GLCi

  ADP kADP+ ATP kATP+ 1   G6P ki,G6P + G6P kG6P + GLCi k_GLCi + 1  (7.1)

After the introduction of hk inhibition, the vmax of the transporter (Tabel 7.1) was increased in

order to restore the flux to the originally published level. This is in agreement with the observation made by Teusink and colleagues, that the glucose transport rate was most likely underestimated by the zero-trans influx assay employed (Teusink et al, 2000).

Secondly, we improved the trehalose and glycogen branch by making the flux dependent on the G6P concentration (k * G6P) (Table 7.1).

Thirdly, since phosphate dynamics is a key feature of the tps1∆-mutant’s phenotype, as we report in this work, we included phosphate ( Pi) as a free variable. Piexchange with the vacuole,

Table 7.1. Parameter adjustments of the detailed kinetic model of Teusink et al (2000)

Parameters New value Original Value unit

vmax ,glt 198 97.264 mM min−1 k_i,G6P,hk 0.07 mM k_trehalose 2.32 min−1 kglycogen 5.8 min−1 npyk 4 L_0,pyk 60000 k_PEP,pyk 0.19 0.14 mM kADP,pyk 0.3 0.53 mM k_FBP,pyk 0.2 mM k_ATP,pyk 9.3 1.5 mM vmax ,pdc 1062.58 174.19 mM min−1 k_phosex 0.1 min−1 P_T 10 mM AXP 3.1 4.1 mM EtOH 0 50 mM GLC_o 110 50 mM 50 100 150 200 250 300 0 20 40 60 80 100 120 140 10 20 30 40 50 60 70 0 10 20 30 40 time (min) FBP (mM) A B kPHO Sex = 0.2 kPHOSex = 0.1 kPHOSex = 0.01 kPHO Sex = 0.2 kPHOSex = 0.1 kPHOSex = 0.01

Figure 7.6. FBP dynamics are dependent on the rate of phosphate exchange Long (A) and short (B) time

identical to the value originally fixed in the model (Eqn (7.3)). This ensures that P_iis a dynamic variable, but reaches the original 10 mM at any global steady state. Since the phosphate, stored as polyphosphates, in the vacuole is so much higher than in the cytosol, we did not model the vacuolar phosphate but rather we kept it as a fixed external metabolite. The rate constant

( kphosex) for this equation directly influences the rate of FBP accumulation in the TPS mutant

(Fig. 7.6). We chose the rate constant such that the FBP accumulation was in agreement with data from Hohmann et al. (Hohmann et al, 1996) , in the regime where phosphate does not seem to be limiting (Tabel 7.1). In addition, we adjusted the rate equation for gapdh to include free phosphate and assumed a Km of phosphate of 1 mM (Eqn (7.3)).

vphosex = kphosex·(PT−Pi) (7.2)

v_gapdh =

C_gapdhvmax ,gapdh,f· GAP· NAD· Pi

kGAP· kNAD −

vmax ,gapdh,r· BPG· NADH

kBPG· kNADH  ( P_i+ 1)  BPG kBPG + GAP kGAP + 1   NAD kNAD + NADH kNADH + 1  (7.3)

Fourthly, the activation of pyk by FBP has been shown to be an important regulatory feature (46) and since FBP potentially accumulates in the tps1∆-mutant, we have adjusted the pyk rate equation accordingly (Eqn (7.4)). We retained the originally measured v_max values from (Teusink et al, 2000) and implemented the other parameters as reported by (Van Eunen et al, 2012) (Tabel 7.1).

v_pyk =

v_{max ,pyk}· PEP

kPEP  PEP kPEP + 1 n−1 kPEP  L0, PYK  ATP kATP+1 FBP kFBP+1 n +  PEP kPEP + 1 n ADP ADP + k_ADP (7.4)

Lastly, we addressed the observation that pyruvate accumulation under certain conditions is caused by an insufficient vmax of pdc (as reported by (Teusink et al, 2000), which should be

increased to reflect experimental measurements. We therefore increased the vmax of pdc by 6.1

fold (Tabel 7.1), as originally suggested, which did indeed resolve this behavior.

The original model was used for predictions of steady state profiles. In this study we are interested in the transient dynamics of the model and therefore we had to make some additional adjustments. Prior to perturbation with a glucose pulse ethanol will not be present or will be very low, and consequently the ethanol concentration is set to 0. This hardly influenced the steady-state fluxes. For the simulations we used initial metabolite values from a study by (Fendt et al, 2010). Specifically, we used metabolite data for S. cerevisiae cells, grown in batch cultures with excess galactose as carbon source and sampled during mid-logarithmic growth (Tabel 7.2); this condition closely resembles our own experimental setup.

Metabolite values were converted to mM (L cell volume−1) concentrations, using a volume of 2 ml per gram dry weight (gDW). This conversion factor is very similar to the value arrived at by (33), and is calculated by assuming a single cell volume of 3 x 10−14L (as in (Van Eunen et al,

1

Table 7.2. Initial values used in the kinetic model based on Fendt et al (2010)

Metabolite Initial value (mM) GLC_i 0.087 G6P 3.085 F6P 0.752 FBP 0.836 BPG 0.111 P3G 0.825 P2G 0.138 PEP 0.140 PYR 0.884 ACE 0.047 TRIO 0.518 NADH 0.044 P_i 101 ATP 2.06 ADP 0.870 AMP 0.165

2010)) and a single cell weight of 1.5 x 10−11gDW·cell−1(47). Missing metabolite values were calculated using the published steady-state flux distribution map (Fendt et al, 2010) and the rate equations from our kinetic model. The sum of AMP, ADP and ATP has also been adjusted to the measurements from (Fendt et al, 2010) (see Table 7.1).

The final model was compared to the original model, and the steady-state metabolite levels and fluxes slightly changed, mostly reflecting the specific changes we made to the model (Fig. 7.7). To mimic the tps1∆mutant in our model, we set the rate constant of the trehalose branch to 0 and removed the feedback inhibition of G6P on hk. Simulations of the tps1∆mutant with our model show canonical metabolite profiles (Fig. 7.8A). While the fluxes within the higher and lower glycolysis modules become equal, the fluxes around FBP and the triose phosphates are unbalanced (Fig. 7.8B), which cause the accumulation of FBP.

Initial condition-dependent outcomes and basins of attraction

GLCi G6P F6P FBP BPG P3G P2G PEP PYR ACE TRIO NADH AMP ADP ATP 1 2 3 7 8 2000 model New model Concentration (mM )

Figure 7.7. Steady-state metabolite levels of the metabolites for the original Teusink et al. model (Teusink

et al, 2000) and our adapted model with extracellular concentrations of 110 mM Glucose The fluxes through

upper glycolysis are 88 (2000 model) and 86 mM min−1in the model presented here, with the fluxes through lower glycolysis at 138 and 133 mM min−1respectively (see Fig. 7.1 for glycolysis scheme). The difference in intracellular glucose is due to the increase of the transporter Vmaxand the compensation by G6P inhibition of hk, the decrease in

PYR concentration is due to the increased Vmaxof pdc and the difference in the total of AMP, ADP and ATP is due to

the adaptation to the Fendt and Sauer data (Fendt et al, 2010). The levels of BPG are very low, but similar.

0.2 0.4 0.6 0.8 1 0 50 100 150 4.5 3.5 mM min -1 Vpfk 1 2(Vg3pdh + Vgapdh) 0.5 1 0 50 100 150 time (min) A B 0 _{2 4 6} 10 0 3 6 5 10 FBP TRIO ATP Pi mM

Figure 7.8. In-silico phenotype of the tps1∆mutant Insets are profiles of wild-type model simulations. (A)

Metabo-lite profiles of simulations of the tps1∆model. The tps1∆mutant phenotype is simulated by setting the flux towards trehalose to 0, and eliminating the feedback on hexokinase. In the tps1∆mutant, FBP accumulates, while ATP and Piare almost depleted. (B) Flux profiles of simulations of the tps1∆model around FBP and the triose phosphates

1 2 3 4 5 0 5 10 15 20 FBP (mM) P i (mM) Steady state time (min) time (min) mM 3 6 0 5 10 Imbalanced state 3 6 0 5 10 mM FBP ATP Pi FBP ATP Pi

Figure 7.9. Wild-type model also shows a co-existence of two stable states Basins of attraction for the model

with the simple feedback of G6P on hk are shown for FBP and Pi, with all other initial conditions kept the same.

The basin is much smaller than in the tps1∆model (see Fig. 7.1). Experimental data (see Fig. 7.2C) suggests an underestimation of the in vivo size of the basin of attraction in WT cells; a discrepancy that is not unexpected given the simplified approximation of the trehalose pathway in our wild-type model.

were taken from Tabel 7.2.

The co-existence of the imbalanced and steady state shown in the main text for the tps1∆

mutant is also present in the model with G6P feedback on hxk - with the effect of the TPS branch represented in a simplified way (Fig. 7.9); this design is comparable to other glycolysis systems, such as in human muscle cells. The co-existence of the two states in the wild-type model shows that this could be a much more general feature of glycolysis, but it will in general be difficult to observe a small nonviable fraction in a population of viable cells. Our FACS data (Fig. 7.2C) on wild-type yeast and tps1∆mutants on galactose do show that such a subpopulation exists; their relatively large size (approximately 7%), suggests that the basin of attraction of the imbalanced state is larger than is modelled here. This discrepancy is not unexpected, as the size of the basin of attraction is dependent on the complex kinetic characteristics of the trehalose pathway, which we have implemented as a simplified approximation.

Supplementary material section 5:

Core model of glycolysis

Model description

GLUCOSE v1 v2 v3 2ATP polyP 2ETOH FBP 2GLYC v4 2Pi 2Pi 4ATP v5

Figure 7.10. Depiction of the core model of glycolysis, representing tps1∆mutants Underlined metabolites are

held fixed. Rates of reactions are denoted by a “v” with the reaction number as an index. polyP denotes phosphate storage. The glycerol branch, shown in grey, is omitted from the analysis detailed in the main text and Fig. 7.11.

core model is depicted in Fig. 7.10. The corresponding differential equations are:

˙

FBP = v₁( ATP)−v₂( ATP, FBP, P_i)−v₅( FBP) ˙

ATP =−2v₁ATP + 4v₂( ATP, FBP, P_i)−v₃( ATP) ˙

P_i=−2v2( ATP, FBP, Pi) + 2v5( FBP) + v3( ATP) + v4( Pi)

(7.5)

We made the dependencies of the reaction rates on the involved substrates and products explicit, e.g. v₁( ATP) indicates that reaction 1 is sensitive to ATP (we will usually exclude the dependencies below). Note that the synthesis reaction of FBP, implemented as a simplified kinetic description of the glucose transporter and hexokinase, is insensitive to FBP, which mimics the situation of a tps1∆mutant strain where regulation of hk is absent and, moreover, pfk is kinetically insensitive to its product, FBP, in yeast (Teusink and Westerhoff, 2000).

Stable and attracting states

We are interested in showing that this system, given suitable kinetics, can reach two states; one of the stable solutions is a normal steady state and the other corresponds to an imbalanced state, in which FBP continuously accumulates while P_iand ATP remain constant. The definition of a steady-state or the equilibrium of a dynamical system is where all derivatives are equal to 0 (in this case FBP = 0,˙ ATP = 0 and ˙˙ P_i= 0). Equilibria can be attracting or repelling, which means that the system evolves towards the equilibrium or away from it (in the figures attracting equilibria are shown as closed dots while repelling equilibria are shown as open dots). The imbalanced state is technically not a steady state or equilibrium, since FBP is continuously increasing and therefore FBP˙ > 0. However, the other metabolites are in steady state ( ˙ATP = 0 and ˙Pi = 0)

and some of these imbalanced states are attracting. 2 The system can evolve towards these attracting imbalanced states and therefore they are biologically relevant.

2

The normal steady-state

A steady state of the entire system, i.e. when FBP = 0,˙ ATP = 0 and ˙˙ P_i = 0, is only possible if v4 = 0, because 2 ˙FBP + ATP + ˙˙ Pi = v4(the sum 2 FBP + ATP + Pidepicts the total phosphate

content of the system). The rate of phosphate exchange (v₄) is 0, because the whole system does not consume or produce phosphate, whereas in the imbalanced state phosphate is contin-uously mobilized with (v₄ 6= 0). Since v₄ depends on P_i only, v₄ = 0 defines the steady-state concentration of P_i, which we denote by P_T. Since v₄= 0 implies a linear combination of differ-ential equations is 0, the steady state concentrations of FBP and ATP can be calculated from two more nullclines:

˙

ATP =−2v1ATP + 4v2( ATP, FBP, PT)−v3( ATP) = 0

˙

Pi=−2v2( ATP, FBP, PT) + 2v5( FBP) + v3( ATP) + v4( PT) = 0

(7.6)

These two equations allow for a graphical exploration of the steady states of the system in the (ATP,FBP)-plane, which we shall illustrate below.

The imbalanced state

One of the curious features of the imbalanced state is that while ATP and Piattain a steady state

value, FBP accumulates. We will show with the core model how this is possible and that the phosphate that accumulates as a component of FBP, leads to a relation between FBP (the FBP˙ accumulation rate) and v4 in the imbalanced state. FBP accumulates in the imbalanced state

(denoted by FBP→∞), therefore, since FBP is a substrate for reaction 2 and 5, the corresponding

enzymes will quickly become saturated with FBP and consequently become insensitive to its concentration. If we write the resulting rate equations as ˆv2( ATP, Pi) = v2( ATP, FBP→∞, Pi)

and V_5,max = v₅( FBP→∞), we obtain the following conditions for the imbalanced state:

˙

FBP = v₁( ATP)−vˆ₂( ATP, P_i)−V_5,max >0 ˙

ATP =−2v1ATP + 4 ˆv2( ATP, Pi)−v3( ATP) = 0

˙

Pi=−2 ˆv2( ATP, Pi) + 2V5,max + v3( ATP) + v4( Pi) = 0

(7.7)

This analysis demonstrates why the insensitivity of v1 to FBP is a key feature of the

imbal-anced state; if v1 would be sensitive to FBP this inhibition would prevent the accumulation of

FBP ( ˙FBP > 0 would not hold), and a regular steady state would result. We can use these equations for a graphical exploration of the imbalanced states in the (ATP, Pi)-plane, which we

shall illustrate below.

Furthermore, in the imbalanced state, we have:−1

2ATP˙ −

1

2P˙i= v1−vˆ2−V5,max−

1

2v4 = 0

(since ATP = ˙˙ P_i= 0), and therefore FBP =˙ 1₂v4. It is logical that the FBP accumulation occurs

at half the rate at which phosphate is released from the vacuole, because FBP contains two phosphate groups per molecule.

We can conclude from this section that in the imbalanced state ATP and P_i can reach a steady state while FBP accumulates with a rate equal to half the vacuolar export rate of phos-phate, and this can only occur when feedback inhibition on v1 is lacking (as in yeast glycolysis).

on hk serves to slow down substrate supply for pfk, making the imbalanced state in the wild-type less likely to occur (see Fig. 7.9).

The importance of phosphate homeostasis

The above insights highlight an important feature of phosphate homeostasis which is central to the imbalanced state. At the normal steady state, there is no net exchange of Pibetween cytosol

and the vacuole, whereas a key feature of the imbalanced state is the continued import of P_iinto the cytosol from the vacuolar storage pool; a process which sustains the accumulation of FBP. Most models of glycolysis do not include Pias a free variable, but rather completely ignore it or

keep it fixed at a measured value. In these cases, there is no moiety conservation of phosphate, because either reactions are not phosphate-balanced or there is an infinite source of phosphate. Such models are unlikely to show an imbalanced state, as the gapdh reaction will have unlimited P_i. Recent models with P_ias a free variable, such as in two lactic acid bacteria (Levering et al, 2012) or in muscle (Voit, 2003), do have a conserved moiety of total phosphate; consequently, glycolytic imbalance can potentially occur, but only if exchange of phosphate with either vacuolar storage or external environment is possible. Without such exchanges, continuous accumulation of FBP cannot occur; imbalance in glycolysis would result in a complete depletion of ATP and P_i, and a stop in FBP accumulation.

A kinetic parameterization of a core model that displays imbalanced states or

steady states depending on initial conditions

In the previous section, we carried out an initial analysis of the core model in both a regular steady state and an imbalanced state. In this section, we will present a kinetic parameterization of the core model that will either end up in the imbalanced state or steady state depending on the initial conditions. In Tabel 7.3, the rate equations and the kinetic parameters are summarized.

In Fig. 7.11, the dynamics of the core model are shown for two different sets of initial condi-tions (specified in Tabel 7.3). The phase plane plots (Fig. 7.11, top panels) show the balanced and imbalanced state as intersections of the graphs d Pi

dt = 0 and dATP

dt = 0. Calculations for the

phase plane plots were facilitated by the use of a variable transformation,ϕ= _{FBP+ k}FBP

M,LG,FBP. This

variable transformation ensured that all time derivatives of internal metabolite concentrations were equal to zero in the balanced and imbalanced state â ˘A ¸S to simplify numerical analysis â ˘A ¸S and that saturating FBP concentrations are defined byϕ= 1. The phase plane plots were constructed by solving the d Pi

dt = 0 and dATP

dt = 0 forϕ(using Mathematica with the Solve

func-tion) for different concentrations of ATP. The flow diagrams were created by calculating the time derivatives for Pi and ATP for different concentrations of Pi and ATP and subsequently using

the Mathematica function ListStreamPlot.

In one case (Fig. 7.11B), the system evolves towards a regular steady state where all the concentrations become constant. Note that phosphate becomes equal to the concentration in the vacuole, which was set to 10 mM (Tabel 7.3). In the second case (Fig. 7.11A), only ATP and Pi reach a steady state while FBP accumulates. Phosphate reaches a level below 10

mM, resulting in the net export of phosphate from the vacuolar compartment, i.e. v4 > 0. In

A B 1 0 2 5 10 ATP (mM) P i (mM) Pi = 0 . . ATP = 0 (i) (ii) (iii) GLUCOSE vupper vlower vATP 2ATP polyP 2ETOH vP 2Pi 4ATP FBP 0 10 20 30 40 5 10 time (min) mM FBP ATP Pi 0.5 1 1.5 0 2 4 6 ATP (mM) FBP (mM) FBP = 0. . ATP = 0 (i) (ii) (iii) 0 _{10 20 30 40} 5 10 GLUCOSE vupper vlower vATP 2ATP polyP 2ETOH FBP vP 2Pi 4ATP

Figure 7.11. Nullcline plots indicate states where the concentrations of FBP, ATP and/or Piare balanced (A)

In the imbalanced state FBP is saturatingly high and therefore we can find the imbalanced states in the plane FBP =∞. (i) The intersection of the nullclines ATP = 0 and ˙˙ Pi = 0 show the imbalanced state where ATP and Piare

balanced. (ii) In this state Pimobilized from the vacuolar stores sustains the continued accumulation of FBP (iii). (B)

Table 7.3. Overview of the kinetic parameterization of the core model of yeast tps mutants.∗Even though this is a core model and not all parameters can be identified with real constants, all affinity constants and maximal rate constants have been assigned realistic values. Units: concentrations are in mM (per liter cytosol) and time in minutes.∗∗In the main text and Fig. S13 the model without the glycerol branch (i.e. , Vmax,5= 0) is used.

Rate Equations

“Transporter +hk +pfk” v1 = Vmax,1 ATP ATP

 1+_K1,i,aATP

 + K1,a

“Lower glycolysis” v2 = Vmax,2_{FBP+ K}FBP_2,f

αT − ATP K2,ADP Pi K2,p  1+ ATP K2,ADP  1+ Pi K2,p  “ATPase” v₃ = k₃· ATP

“Vacuole phosphate export” v4 = k4( PT− Pi)

“Glycerol branch” v5= Vmax,5_{FBP+ K}FBP_5,FBP

Kinetic parameters∗

“Transporter +hk +pfk” V_max,1= 10, K_1,i,a = 3, K_1,a= 0.1

“Lower glycolysis” Vmax,2= 10, K2,f= 1, K2,ADP=0.1, K2,p= 2,αT = 5

“ATPase” k3= 10

“Vacuole phosphate export” k₄ = 0.3, P_T= 10 “Glycerol branch” V_max,5= 1∗∗, K_5,f= 1

Initial conditions

Regular steady state FBP = 2, ATP = 1, P_i= 10 Imbalanced state FBP = 0.01, ATP = 0, Pi= 0.01

in agreement with the experimental data. The initial condition-dependent behavior, shown in Fig. 7.11, illustrates the co-existence of two stable states in the core model, as was found in the detailed model (main text).

Graphical exploration of the steady states

In the steady state we have the condition 2 ˙FBP + ATP + ˙˙ P_i= v₄( P_i) = 0 which implies P_i= P_T. Therefore we can show the isoclines ˙Pi= 0 and ATP = 0 in the plane P˙ i = PT to determine the

concentrations of ATP and FBP in the steady state (Fig. 7.11B). A similar analysis can be done for the imbalanced state (Fig. 7.11B). In the latter case, the concentration of FBP is set to infinity such that the rate equations become saturated and are not dependent on FBP, as explained above. The steady state concentrations for P_iand ATP can then be determined from the ˙P_i= 0 and ATP = 0 nullclines (Fig. S13A).˙

1 2 0 5 10 Pi. .

= 0

ATP

= 0

ATP (mM) P i (mM) 0 _{10 20 30 40} 5 10

FBP

ATP

Pi

time (min) mM C D mM A B 0 10 20 30 40 5 10

FBP

ATP

Pi

0 0.1 0.2 0.3 5 10 0 0.05 0.1 1 Pi. .

= 0

ATP

= 0

P i (mM) time (min) ATP (mM)

Figure 7.12. Enhancement of ATPase or glycerol branch activity leads to disappearance of the metabolic imbalance (A) When ATP hydrolysis is increased (6 fold) the intersection of the nullclines at saturating FBP

dis-appears, consequently an imbalanced state no longer exists. Although not visible in (A), a steady state exists with ˙

ATP = ˙Pi = FBP = 0. This steady state is shown as solid lines in (B), for initial values that would lead to an im-˙

balanced state when kATPaseis low (dashed lines). (C) and (D) Increasing the maximal activity of the glycerol branch

Acknowledgements

We thank G. Smits for pHluorin plasmids; J. Thevelein for W303-1A strains and HXK2 overex-pression plasmids; M. Walsh for fruitful discussions; D. Molenaar and H. Hoefsloot for valuable assistance with statistical methods; J. Pronk for experimental suggestions; K. Hellingwerf, R. Leurs, and I. Stulemeijer for critical reading of the manuscript; and L. da Cruz, Z. Zhao, and A. Deshmukh for assistance with tracer experiments. Funding: This work was supported by fund-ing from AIMMS, FP7 UNICELLSYS, and Kluyver Centre for Genomics of Industrial Fermentation and NCSB, funded by the Netherlands Genomics Initiatives