Forward design of a complex enzyme cascade reaction

Enzymatic reaction networks are unique in that one can operate a large number of reactions under the same set of conditions concomitantly in one pot, but the nonlinear kinetics of the enzymes and the resulting system complexity have so far defeated rational design processes for the construction of such complex cascade reactions. Here we demonstrate the forward design of an in vitro 10-membered system using enzymes from highly regulated biological processes such as glycolysis. For this, we adapt the characterization of the biochemical system to the needs of classical engineering systems theory: we combine online mass spectrometry and continuous system operation to apply standard system theory input functions and to use the detailed dynamic system responses to parameterize a model of sufficient quality for forward design. This allows the facile optimization of a 10-enzyme cascade reaction for fine chemical production purposes.

: Measured concentration time series and s imulation with the best identified parameter set. The start of reaction by enzyme injection is indicated by the stippled red line. The only observed concentration is GLC (blue line), the concentration of other compounds is a result of the simulation. Stippled lines in blue or green indicate the simulated feed profile for GLC (blue) or ATP (green). Column [c] Error distribution achieved by the optimizer in 100 runs. The number in green at the left of each picture indicates the number of runs with an error which was not more than 1% larger than the smallest achieved error. The 1% criterion is indicated by the red line. The number in red at the right of each picture indicates the number of runs with an er ror more than 1% larger. Column [d]: Mean and standard deviations of the parameters that led to a qual ity of fit within 1% of the minimal error. C: Experimental confirmation for an inhibition of Glk by ADP. Recorded concentration time series for GLC in the presence or absence of ADP or G6P as possible inhibitors already in the feed. Mentioned concentrations refer to the initial compound concentrations in the reactor at the time point at which the experiment was started (addition of 1 U of Glk after 3 min) and to the concentrations in the feed. Note that the addition of only G6P does not produce a difference in the concentration time series as compared to the control, suggesting that G6P is not an inhibitor.  Table  8), followed by the concentration time series (measured: blue, best-fit simulation: green). The parameters for the best-fit simulation were obtained from Supplementary Table 9. A: C1; B: C2; C: C3; D: C4. Note that for experiment C3, simulation with the parameter set of Supplementary Table 9 showed a lower fit quality than the other experiments. This can be explained by impurities of the enzymes as can be seen for the injection of 1 U of Gdh after 15 min: theoretically only BPG should be formed but instead one can see an increase in both XPG and A TP, both products of Pgk. This indicates that the applied preparation of the enzyme Gdh might have been contaminated with Pgk.  (see also Supplementary Table  10), followed by the concentration time series (measured: blue, best-fit simulation: green).
The parameters for the best-fit simulation were obtained from Supplementary Table 11 (see also Supplementary Table  12), followed by the concentration time series (measured: blue, best-fit simulation: green). The parameters for the best-fit simulation were obtained from Supplementary Table 13. A: E1; B: E2; C: E3; D: E4; E: E5. Please note: In experiment E5, the XAP dynamics at 60 min showed an increase in DHAP in response to Eno injection. Until then we never had observed such an e ffect and al so in the experiments in the following sections this type of dynamics never occurred again. Hence we considered this a measurement artefact which we did not follow up.  (see also Supplementary  Table 14), followed by the concentration time series (measured: blue, best-fit simulation: green). The parameters for the best-fit simulation were obtained from Supplementary Table  15 Table 17) of subsystem G was performed by re-using experiments B1 to B6 from subsystem B and including them in the estimation. Shown are the data for B1 to B6 (experimental data as in Supplementary    Figure 13: Illustration of the data processing after calculation of concentrations from the specific quantification method. Idealized behavior is shown for two compounds, one of which is present in the feed (upper panel), and one of which is not present (equivalent to a f eed concentration of 0, lower panel). Data set c1 contains the compound concentrations (black lines) calculated with the help of the different standards after correcting for ADP/ATP crosstalk. Data set c2 contains the compound concentrations after correction of offsets (blue lines). Next, a t ime dependent correction factor was calculated from the ideal and real mass balances (symbolized by the lines deviating from the constant balance after the enzyme injection in the middle panel). To obtain data set c3, this correction factor was applied to compounds with initial concentrations of zero, resulting in the pink curve (only in the lower panel). In data set c4, the correction factor was additionally applied to the remaining compounds (initial concentration larger than zero), resulting in the pink curve of the upper panel.  feed   c1  c2  c3  c4  balance  and  correction  factor  Supplementary Tables   Compound  Enzyme  abbreviation Full name  abbreviation Full name  GLC  glucose  Glk Table 3: Commercial sources of used enz ymes and key elements of the definition of activity applied by the supplier. One unit of activity corresponds to the formation of 1 µmol min -1 of product at the given pH and temperature for the reaction for which the enzyme is responsible in the glycolysis. According to the product description the rate was either determined in the direction of standard glycolysis ("forward") or in the opposite direction ("reverse"). Indicated are the type of experiments (first column, "Type of experiment"); an i nternal identifier of the specific experiment (second column); the composition of feed 1 (provided from the first HPLC pump; please note that the composition of feed 1 is equal to the initial reactor concentration, third column); the composition of feed 2 (required to produce gradients, all of which were linear, fourth column); the amount of enzyme injected at t=0 to start the experiment (fifth column); time points of events (sixth column); and the nature of the events (seventh column). MgCl 2 was provided in all cases in excess of the maximal ATP concentration in the feed. Indicated are experimental type, number, feed concentrations of compounds (equals initial concentration in the reactor), injected amount of enzyme at t= 0 and t iming of pulses. For compound pulses the (theoretical) reactor concentration and for enzyme pulses the absolute activity added to the reactor are given.   19 for either the specific host ("Host range") or a set of hosts ("General range"). Parameters without indicated boarders were not estimated. Not considered in the range analysis were parameters measured for mutated or immobilized enzymes as well as singular extreme values. Green indicates K M values within host range and bl ue K M values inside the general range. Pink indicates that the estimated K M -value was identical with either the upper or the lower set boundary. Orange indicates that the estimated K M -value was out of range of the general range of K M -values obtained from BRENDA.  Table 21: Optimized full reaction system. Activity distributions of the G3P production system as used for equi-activity reference and optimized experiment of 20 U and optimized 40 U experiment. The design included an ex perimental time of 60 min, feed concentrations of 2 mM GLC, 2 mM ATP, 1 mM NAD and 2 mM P i , and initial pulses of the required enzymes to the reactor.

Enzyme
[PEP]  [NADH] [NAD] Literature research for E. coli glycolytic enzymes was undertaken and rate equations were derived considering known mechanism, sequence of metabolite binding, and effectors. In cases where enzymes were obtained from organisms other than E. coli, the literature was consulted for deviations in mechanism and effectors. In general the enzymes from the different host organisms had comparable properties to those from E. coli. This indicates that the metabolic regulation of the very basic pathway of glycolysis, mainly by Pfk and Pyk, seems to be conserved within our scope of hosts and application.
Reactions were modeled as irreversible if the calculated thermodynamic equilibria were more than 1000-fold on the product side (in glycolytic direction). For the majority of reactions, rapid equilibrium Michaelis-Menten (MM) type kinetics was applied. In the case of the enzyme Eno (catalyzing a dehy dration reaction) we did not consider water as reactant as in aqueous solutions it is ubiquitous. Enzymes G3d and Ldh, which were shown to use a ping-pong mechanism, were modeled as Theorell-Chance (TC) kinetics. For the two cooperative enzymes with allosteric interactions Pfk and Pyk the Monod-Wyman-Changeux (MWC) model was utilized. As in these cases cooperativity is only towards one of the two substrates, the non-cooperative substrate was modeled as MM kinetics and the rate equation was completed by multiplying the MWC and MM terms. Special care was taken to include regulatory properties of the pathway members. Complex regulatory mechanisms for the allosteric enzymes Pfk and Pyk are described in the literature. These include for Pfk an allosteric activation by all nucleoside diphosphates (due to the model scope only ADP was considered in the model), competitive inhibition of ADP by ATP, and an allosteric inhibition by PEP 4 . All of these effects were included in the model. For Pyk, these include allosteric inhibition by ATP, GTP, and succinyl-coenzyme A as well as allosteric activation by FBP 14 . Due to the model scope only the interactions of ATP and FBP were integrated into the model. Kinases use the Mg 2+ -nucleotide complex as enzyme substrate, meaning that MgATP is the real substrate for Glk, Pfk, Pgk and Pyk. Our model does not directly account for this fact, but we prevent the Mg 2+ concentration from influencing the experimental data by supplementing Mg 2+ in excess. An overview of the regulation pattern is provided in Fig. 2B.
The reaction mechanisms including sources and the formulated rates laws are shown in Supplementary Table 24.

Abbreviations
Abbreviations of compounds and enzymes are summarized in Supplementary Table 1.

Source of chemicals
Chemicals were obtained from Sigma-Aldrich (Buchs, Switzerland) in the highest available purity unless mentioned otherwise. Labelled GLC [U- 13 C 6 ] and ATP [U-15 N 5 ] were obtained from Cambridge Isotope Laboratories (Tewksbury, MA, USA).

Setup of experimental system
Buffers Operation of the reaction system required specific buffers that were compatible with enzyme activity and the subsequent analysis by MS. The basic reaction buffer for the experiments in the enzyme membrane reactor consisted of 50 mM NH 4 HCO 3 , 2.5 mM MgCl 2 and 0.4 mM MOPS at pH 7.8. For enzyme reactions, the reaction buffer was modified where required in order to ensure enzymatic reaction: addition of 1 mM dithiothreitol (DTT) for the stabilization of Gdh and Pgm (as recommended by the suppliers), 1 mM KCl (Pyk requires potassium as a cofactor), and 2 μM 23BPG (the Pgm variant that we used required 23BPG for activation 11 ). Where indicated MgCl 2 concentration was adapted to ATP concentration to ensure formation of MgATP as substrate for kinases. Depending on the quantification methods (see below) standards were added to this buffer. After the reaction mix had left the reactor, an MS matrix buffer was used to dilute the reactor effluent in order to condition it for MS-analysis (see Fig. 1A of the main text). This buffer consisted of 5 mM NH 4 HCO 3 in a m ix of 75% MeOH/25% water (vol/vol). The buffer was adjusted to pH 4.8 by the addition of pure formic acid (typically 5 mL L -1 ) and degassed before usage. Again depending on which quantification method was applied, standards were added before usage.
Hardware The experimental system consisted of a f unction generation part, a part in which the enzymatic reaction system was operated, a sample conditioning part, and the final measurement part (Fig. 1A). Function generation was achieved by combining a ( multichannel) HPLC pump and an injection loop. The HPLC pump (Agilent, Santa Clara, CA, USA) was connected through polyether ether ketone (PEEK) HPLC capillaries to the main reactor and operated at a rate of 0.25 mL min -1 with reaction buffer and reaction compounds as detailed for the individual experiments. Between pump and reactor we placed an injector with loops of different but defined volumes (varying from 20 to 1,000 µL), thus allowing redirection of the flux from the pump through the loop and f lushing the loop content, e.g. enzymes or compounds, into the reactor. Supplying defined compound concentrations in the reaction buffer established a c ontinuous reaction feed; switching abruptly from one to a second reaction buffer allowed implementing step functions. Switching gradually over time allowed implementing a ramp function and even more complex input functions such as a sine can be r ealized by programming the pump. Using the injector for compounds allowed approximating an impulse by producing a pulse from an injector loop volume small in relation to the flux provided by the pump. Finally, using it for enzyme supply produced a step function in terms of reactor enzyme concentration. Reactions were carried out in a continuously operated enzyme membrane reactor (Bioengineering AG, Wald, Switzerland, nominal internal volume of 10 mL) with a c ustommanufactured inlay that reduced the volume to approximately half and t herefore increased the system dynamics. The vessel was constantly stirred by a magnetic stirrer at 600 rpm and the temperature was maintained at 30°C by flushing the integrated cooling jacket with water from a t hermostat. The reactor contained a c ellulose triacetate membrane with a c utoff at 20 kDa to retain the enzymes but to let the compounds pass. Bovine serum albumin was injected to a reactor concentration of 1 mg mL -1 prior to injecting enzymes for improving enzymatic stability (presumably by passivating critical surfaces, e.g. on the membrane). The reactor effluent was conditioned for subsequent online MS detection by diluting it between 30-and 100-fold with the MS matrix buffer, supplied via an HPLC PEEK capillary system. Afterwards, the flow was split in such a way that approximately 0.03 mL min -1 were directed into the MS. The exact fluxes were calculated from standard compounds as discussed below.

Mass spectrometry
In order to measure compound concentrations in the reactor effluent, mass spectrometry was used. Specifically, a triple quadrupole mass spectrometer (MS) (MDS Sciex 4000 Q-TRAP, Applied Biosystems, CA, USA) with electrospray ionization (ESI) was employed in multiple reaction monitoring (MRM) setting and in negative ion mode. The ESI-MS was generally operated at an electron spray temperature of 200°C, an ion spray voltage of -4,200 V, a channel electron multiplier (CEM) of 2,100 V, and the following gas pressures: curtain gas 30 psi, collision gas 6 psi, ion source gas one 30 psi and ion source gas two 40 psi. The identification of compounds by mass required appropriate settings for the three quadrupoles Q1 to Q3, which were responsible for electron spray ionization, analysis and separation (Q1), collision for further fragmentation (Q2), and fragment analysis (Q3), respectively. The corresponding operating parameters required optimization with respect to signal intensity. For this, the compounds to be analyzed were dissolved to a concentration of 20 µM in 50% methanol/50% water (vol/vol) and injected into the instrument using a syringe pump with a volume flow rate of 30 µL min -1 . First, to initiate the analytic MS-routines, the Q1 spectrum was checked for the presence of molecular ion of the compound to be identified, which was always found to be the strongest signal. On the basis of this ion the instrument manufacturer's characterization optimization routine was started, and the optimization routine scanned through the instrument settings (entrance potential EP, declustering potential DP, collision energy CE, and collision cell exit potential CXP) to identify the optimal parameter set. Results are summarized in Supplementary Table 2. These settings (adapted to the differences in masses between compound and isotopologue) were also applied to labeled compounds. The routines described above did not allow resolving structural isomers, specifically the pairs G6P/F6P and 2P G/3PG. Instead, we measured and r eport the sum of concentrations of these pairs in two pools, termed X6P (for G6P/F6P) and XPG (for 2PG/3PG). In an earlier publication, we had reported separate measurements of DHAP and GAP based on a unique molecular ion for DHAP and a persistent negligible concentration of GAP in a phosphate-ion rich background, which tempered the effect of ion suppression 20 . In the present system with a much reduced ion-background, measurements of this DHAP-specific ion did not comply 53 with our accuracy requirements, and we discontinued this practice. Instead, we determined concentrations of DHAP and GAP as a third pool termed XAP. Finally, there was no standard available for 13BPG, therefore we relied on 23BPG.
In the final measurement protocol the dwell time for measuring a single compound in MRM mode was set to 500 ms resulting in a measurement frequency of 1 Hz for a compound and its standard. Depending on the size of the enzymatic system and the number of compounds the measurement protocol was reduced in order to lower the time for a measurement cycle . The EP can be set only once for one protocol and it was set to -10 V.

General microbiological and molecular biology methods
General microbiology and molecular biology work followed established standard procedures unless otherwise mentioned.

Production of the E. coli glucokinase Glk
The E. coli glk gene was expressed from plasmid pET30-Glk, a pET-30b(+)-derivative (Novagen) that carries the gene encoding an N-terminally 6xHis-tagged Glk under the control of the T7 RNA polymerase-dependent promoter. Cells were prepared after induction from 1 L of LB medium by centrifugation, resuspended in 10 mL of lysis buffer (50 mM Tris/HCl pH8, 300 mM NaCl, 10 mM imidazol) and lysed in a high pressure homogenizer (Microfluidics M-110P, Westwood, MA) at 4°C and a 1,000 bar pressure drop across the orifice. The clarified supernatant was applied to two self-packed Ni-Sepharose gravity flow columns (2 mL bed volume, GE Healthcare, UK) and eluted with 400 mM imidazol. The eluted protein was dialysed twice against a 200-fold excess of Tris/HCl pH 7.5 containing 1 mM MgCl 2 . The protein was 95% pure (SDS-PAGE, visual inspection). Aliquots were stored with 10% glycerol at -80˚C. The final enzyme concentration was determined as 4.8 mg mL -1 by photometric measurement. The enzyme activity was determined spectrophotometrically by oxidation of glucose 6-phosphate with concomitant reduction of NADP. The specific activity of the Glk preparation was determined as 176 U per mg of purified protein.

Production of cell-free extract (CFX)
For CFX production (used for the synthesis of isotopologues, see section 10 of Supplementary Methods), we constructed the strain E. coli BW25113 ΔfbaB pfkB pykA::kan tpiA::Cm by phage transductions based on strains from the KEIO-collection 21 . The resulting strain was then used for the production of CFX as described previously 20 . CFX was stored in aliquots at -80°C.

Identification of experimental system
In order to separate the identification of the parameters of the system setup, e.g. volume of the reaction vessel, from the identification of the parameters of the enzymatic system, we determined the transfer function of the experimental system first. The experimental system consisted of the HPLC pump, the injector, the reactor and the dilution device. For this, the system was decomposed into subunits which were identified one after the other (Fig. 1C). A constant flow through the subunit that was to be identified was provided by a syringe pump (for the identification of the injector) or by one channel of the HPLC pump (for all other identification experiments). We applied 1 m M GLC in reaction buffer to generate a perturbation. A step input function was impressed on the system through either switching the HPLC pump to a second channel or to an injection loop with a volume of 1 mL. The response of the system to the step function was recorded by the MS. Neglecting the setup-specific dead times, the injector response was close to an ideal step. For the subsequent identification steps the injector function was taken as an i deal step function. The other system elements generally behaved as first order lag elements (PT 1 ) with dead time. Again disregarding the dead t ime as a s etup-specific property, their responses were fitted to a continuous stirred-tank reactor (CSTR) model:

Compound quantification
The signal generated by the mass spectrometer (detected ions over time) (see section Mass spectrometry) had to be linked to the concentration of the specific compound in the reactor. There were two potential main error sources in this process: First ion suppression 22 due to the presence of a variety of different compounds in varying concentrations over time in the reactor and therefore also in the ESI-chamber; and secondly flux instabilities in the passive mixing and s plitting steps in the devices after the reactor, where the conditioning of the reactor efflux for MS-analysis relied on a s ystem of thin capillaries and the generated flux pattern depended on pr essure drops defined by diameters and l engths of the applied capillaries. In consequence, the flux pattern was sensitive to influences changing the involved pressure drops, such as partial blockage of capillaries by solid particles introduced from the reaction buffer or accumulation of drained enzyme (fragments) or other compounds on the capillary walls. We therefore implemented a s et of 3 pr ocedures allowing us to compensate for these effects and to quantify the concentrations of compounds in the reactor based on the MS signals. First, flux monitoring based on the addition of two different isotopologues of taurine (which is not a par t of the reaction system) to the reaction and t he MS matrix buffer enabled the determination of the flux ratios post-reactor and in consequence also the dilution over time of the reactor effluent with MS matrix buffer (Supplementary Fig. 2A and B). Secondly, as ion suppression effects should apply in the same way to isotopologues of the same compound, referencing to isotopologues should lead to accurate measurements of concentrations 23 . Therefore, once the flux ratio was known, the concentrations of a compound in the reactor effluent could be calculated by comparing them to the concentration of its corresponding isotopologue in the MS matrix buffer (Supplementary Fig. 2C). Thirdly, where an isotopologue was not available, another standard (HEPES) that was not part of the reaction system was added to the MS matrix buffer and metabolite concentrations were calculated from predetermined calibration functions of compounds to the used standard ( Supplementary Fig. 2D). The methods will be described in more detail in the following sections.

Determination of the flux ratio of reaction effluent and MS matrix buffer
The flux ratio was calculated as follows (see Supplementary The mole flux of compound i is the product of the flow rate F and its concentration i c , and so Supplementary Equation (4) Inserting the balance for the flow rates into Supplementary Equation (5) Using the dilution factor R D for the dilution of the reactor efflux reactor R reactor dilution (8) and D D of the MS matrix buffer flux Supplementary Equation (7) can be simplified to TAU with the standard 14 The parameter TAU y is defined as the ratio of the MS-determined signals x for 14 N-and 15 N- Assuming that both labeled and unlabeled taurine ionize the same way, TAU 57 Using the Supplementary Equations (11) and (12) this can be written as: 15 15 (15) and Supplementary Equation (15) Two assumptions were made in the deriving of this equation: [i] The ionization properties of TAU are not influenced by the heavy N-isotope; and [ii] the only source for unlabeled TAU at the mixing tee is the reactor effluent and the only source for labelled TAU is the MS matrix buffer. Regarding the former assumption, injection of TAU and 15 TAU at identical concentrations into the MS via the syringe pump resulted in the same signal intensity (data not shown). Regarding the latter assumption, the certificate of analysis by the supplier for 15 TAU stated a pur ity of 99.3% 15 N (corresponding to 0.7% 14 N). The natural isotope distribution of TAU is known (99.636% 14 N and 0.364% 15 N). At the targeted dilution rate of 1:100, the TAU mole flux from the reactor is essentially identical to the 15 TAU mole flux from the MS matrix buffer, as the TAU flux from the reactor is diluted 100-fold. Hence the isotopological impurities partially compensate each other restricting the error to the difference in purity of 0.336%. In sum it can be said that both assumptions are reasonable.

Concentration determination by reference to isotopologues
As ion suppression effects should apply in the same way to isotopologues of the same compound, referencing to isotopologues should lead to accurate measurements of concentrations 23 . To calculate the reactor concentration of a c ompound using its isotopologues as a reference ( Supplementary Fig. 2C), Supplementary Equation (16) can be re-formulated for each compound i whose non-labeled version is present in the reactor effluent and whose heavy isotopologue is present in the MS matrix buffer: 58 Now the combination of Supplementary Equations (17) and (18) We further analyzed the quantification range for different compounds (GLC, FBP, AMP, ADP, ATP, and t he measurement pools of X6P, XAP and XPG), and found the relations to be linear under all relevant conditions over more than 3 orders of magnitude (data not shown).

Concentration determination in reference to a standard in the MS matrix buffer
In order to calculate compound concentrations in the reactor in reference to the HEPES standard ( Supplementary Fig. 2D), the ratio i z of the number of ions per second i x of compound i detected at the MS and HEPES x of HEPES was taken: The relation between metabolite concentration at the instrument i,MS c and the ratio i z was taken from a c alibration curve that had been obtained as follows: A mix of all relevant compounds at a defined concentration was prepared, distributed over aliquots, and stored at -80°C. Reaction buffer was prepared freshly and di luted 1:100 with MS matrix buffer containing 2 µM HEPES to give the calibration solution. Then an aliquot containing metabolite mix was thawed and diluted 1:100 with calibration solution. This initial metabolite mix was then repeatedly diluted two-fold with calibration solution, so that the metabolite concentration decreased quadratically while the concentrations of the background compounds in the reaction buffer (e.g. MgCl 2 ) and in the MS matrix buffer (HEPES) remained constant. In this, samples corresponding to reactor concentrations between 0.61 nM and 4 mM were prepared. These samples were injected directly into the MS with a syringe pump at a rate of 30 µL min -1 which was in the range of typical fluxes of the dilution system to the MS. For each metabolite equation (20) was evaluated at different metabolite concentrations and these data points were used for the identification the parameters a and b of a quadratic These calibration curves describe the difference in ionization properties between HEPES and the metabolite which, as already suggested by the equation above, was a non-linear relation.  (22) and, with the dilution term defined as in Supplementary Equation (8) with the flux ratio from the TAU flux monitoring as given in Supplementary Equation (17).

Consideration of errors
The 15 TAU/TAU standard was used in general to observe the flux pattern. While the flux from the MS matrix buffer was constantly maintained at 1 mL min -1 by an HPLC pump, the flux from the reactor was prone to variation. During experiments we typically observed dilution rates in the range of 1:25 to 1:100. This supports clearly the idea of measuring relative to a standard that is provided by the MS matrix buffer as its corresponding dilution factor was close to 1 (  (9)). Therefore the maximal difference in concentration at the MS of any compound added to the MS matrix buffer is 3%, which for all practical purposes means the signal is constant. Hence the concentration of a standard, such as HEPES, added to the MS matrix buffer can be used as a reliable reference signal in the measurement system. By comparing the concentration measurements via the isotopologue standard with those that would have resulted if these concentrations had been evaluated using the HEPES standard for those components for which both methods could be appl ied, we concluded that the remaining ion-suppression-derived error in the HEPES-method is smaller than 20% for the conditions applied in the experiments described here. The relation between ionization properties of the standard and the compound of interest are dependent on the instrument. As an example the position and also the wear level of the electrode used in the ESI have an influence on ion production. We considered these effects by renewing calibration curves in a biweekly manner or whenever changes to the instrument were made. In view of the remaining limitations of the HEPES-based quantification, we limited the use of this method to those situations where no suitable isotopologues were available. This was the case for some of the compounds that were unique to the lower glycolysis, and therefore those experiments that involve sections of the lower glycolysis were evaluated with the HEPES-standard.

Production of isotopologues
As mentioned, we had only heavy isotopologues for ATP and G LC commercially available and therefore proceeded to produce further isotopologues ourselves using CFX in the enzyme membrane reactor setup. The system was provided with reaction buffer containing 2 mM 13 GLC, 2 mM 15 ATP, 4 mM MgCl 2 , 2 mM sodium phosphate (Na 2 HPO 4 / NaH 2 PO 4 , adjusted to pH 7.2), 1 mM NAD and 400 µM 15 TAU. NAD was provided in unlabeled form.
The MS matrix buffer that was used for dilution after the reactor contained 10 µM unlabeled glycolytic compounds and 4 µM TAU for monitoring and quant ification purposes. The reaction was started by injection of 7 mg of total protein derived from CFX via the injection loop resulting in a total protein concentration in the reactor of 1 mg mL -1 . The reactor was continuously fed with reaction buffer at 0.25 mL min -1 for 2.5 h. Once the reactor had reached steady-state, the reactor effluent was directly (without further dilution by MS matrix buffer) collected first on i ce and then stored in aliquots at -80°C. In order to account for possible metabolite degradation through freeze and thaw processes in subsequent experiments, a thawed aliquot was always measured against a standard of unlabeled glycolysis compounds and concentrations were determined.  15 ATP, 0.955, 1.910. NAD, NADH: n.a. c) ; n.a. c) . 15 TAU; 0.395; 0.790. n.d. not determined; n.a. not applicable a) Concentration in self-prepared set of labelled compounds remained below or at limit of detection b) The 13 LAC signal was strongly perturbed by an unidentified ion c) No non-radioactive isotopologues of NAD and NADH were available.

Enzyme stability
The stability of enzymes is a prerequisite for successful experiments and for the validity of the kinetic model which was applied, which is time-invariant in the amount of enzymes. However, enzymes are prone to degradation over time or could be lost through leakage over the membrane in the reactor. We addressed the evaluation of a pos sible loss in enzyme activity repeatedly in the course of this work. However we did so indirectly. Through whatever underlying effect, a loss of activity would result in the present reaction system in a change of metabolite concentrations over time under steady state conditions. However, we did not notice such system behavior in any of our experiments, which typically did not exceed 120 min after charging the reactor with enzymes. As an example for the long time stability of reactions we refer to the stability in metabolite concentrations in three different reaction systems as shown in Supplementary Fig. 3: first Glk only (panel A), secondly upper glycolysis (including Glk, Pgi, Pfk and Ald, panel B) and thirdly full glycolysis (panel C). We conclude that no s ignificant loss of activity resulted in the experimental system for at least 120 min.

Calculation of thermodynamic equilibria
Reaction equilibria are a function of a v ariety of influences, such as pH, temperature and ionic strength. In order to decide on the reversibility of enzymatic rate laws and to compute a starting point for the estimation of selected equilibrium constants of the model, we used the biochemical thermodynamics calculator eQuilibrator 1.0 (equilibrator.weizmann.ac.il) 24 . For our reactor conditions we estimated the equilibria for the minimal and maximal value of ionic strength that was deemed applicable. The results are summarized in Supplementary Table  23.

Kinetic model
The kinetic model was constructed by balancing the compounds, formulating rate laws for the enzymes, and e mbedding them into the identified 3 v essel system of the experimental setup. The kinetic model describes the (essentially glycolytic) reaction network displayed in    (25) Optimization In order to search the parameter space for parameter estimation and optimizing enzyme concentrations a Gaussian adaptation algorithm 25 served as optimizer. The used algorithm is a probabilistic method resulting in different outcomes for different optimization runs and so by default multiple runs were performed to ensure optimal results. For the estimation of parameters p , the mean squared error function between measurement and simulation was minimized as a quality criterion:

Concentration balances
From a practical perspective we wanted to optimize the distribution of a given amount of total enzyme 0 E (in units) resulting in a constrained optimization problem with as the constraint. We realized this by optimizing the proportions i k of enzymes instead of directly optimizing the enzyme amounts i E to fulfill the conservation of Supplementary Equation (28). For example, for a two enzyme system the relation between the first and the second enzyme would be formulated as and the conservation from Supplementary Equation (28) would be In the optimization only the proportion 1 k would be adapted while inserting Supplementary Equation (29) into (30) allows the calculation of the last, here second, enzyme amount: This process can be expanded for n enzymes in n-1 calculations.

Data processing and error model
In the section "Compound quantification" we presented the strategies and m ethods for the conversion of MS signals into concentrations in the reactor. The practical implementation of these strategies required a work-flow that was able to deal with the imperfection of real data. Different types of errors like crosstalk between compounds, offsets in the measurements and errors in quantification were addressed in the form of a pipeline -a fixed sequence of steps to treat data consistently. The general workflow of the data processing pipeline is illustrated in Supplementary Fig. 13 and included the following steps:

Elimination of crosstalk in raw data
First of all we corrected the raw data from crosstalk between adenosine nucleotides. Injection of ATP of the highest purity available commercially (≥ 99% according to the certificate of analysis) into the MS resulted in signals for ADP and AMP that could not be solely explained by the remaining impurities in the ATP preparation. We therefore assume that this phenomenon is a consequence of the electron spray ionization process, during which one or two phosphate groups could be l ost from ATP, leading to the in situ creation of two other compounds, one of which is part of the experimental measurement. During experiments, the ATP-signal from the feed-only phase (before the injection of the enzyme(s)) was used to calibrate this percentage anew for each experiment. In the following experiment, the ADP-signal was then corrected as a function of the ATP signal according to this calibration. We did not observe similar behavior for any other metabolite.

Calculation of metabolite concentrations
The metabolite concentrations in the reactor were calculated according to the different quantification methods as described earlier. These raw concentration time series were stored as data set "c1" for further processing (see also Supplementary Fig. 12 at the bottom of this section).

Matching the observed and set concentrations
Where these calculations led to non-closing mass balances during those times where the system was under complete control (e.g. before the addition of enzymes to start the reaction), these differences were used for corrections, leading to data set "c2". Specifically, signals of compounds that were sure to be absent in the feed at the start of the reaction were taken as offsets and s ubtracted from calculated concentrations. These typically low errors can be t racked to different types of contamination: remaining quantities within the experimental system (in particular the post-reactor tubing), pollution of the MS and impurities of used compounds. For compounds with initial experimental concentrations larger than zero, the ratio between the set and the measured concentration  c t kc t = . (33) Possible reasons for these differences for compounds with initial concentrations larger than 0 have to be distinguished according to the quantification method: For the HEPES standard based quantification method, ion suppression by different ion backgrounds is the most likely source of error. In agreement with this, the observed error in the compound balance was typically much lower for the isotopologue standard. This presumable error source suggested that a scaling factor was a more suitable correction than a constant addition or subtraction as applied above. However, even for the isotopologue standard the error was not always negligible. Possible reasons included instrument tuning-specific limitations with the exact distinction between TAU and 15 TAU at low dilution rates in the post-reactor tubing section (i.e., when dilution deviated substantially from the 1:100 target value, the two signals might overlap in the MRM analysis) and pos sible errors in the quantification of the labelled standard itself, which needed to be repeated regularly and showed distinct variations as a function of the continuous use of the MS. Furthermore, for economic reasons we used the isotopologues in the corresponding standard in the lowest possible concentration, which allowed a good and reliable quantification of the reactor-concentrations in most scenarios, but was susceptible to contaminations in some situations, the reasons for which were not always possible to discern.
For experiments involving only the enzymes of upper glycolysis with Ald as the last enzyme of the reaction system, the concentrations of DHAP and GAP were obtained by dividing the measured concentration pool XAP into two equal parts as both compounds are produced with identical stoichiometry by Ald.

Calculation of compound balances
In order to further improve data consistency of the observed concentrations we took into account the conservation of species, which had to be computed. For the most general casethe full reaction system -four conserved moieties were selected, reflected in the carbonbalance, the NAD(H)-balance, the AXP-balance and the P i -balance. The carbon balance can be written down in different ways, e.g. taking either compounds with 3 carbon atoms (C3) or with 6 carbon atoms (C6) as basis, but these can be transformed into each other by a factor. Written for C6, the carbon balance can be formulated as Although these 4 formulated balances hold true generally, not all balances were helpful or meaningful in each scenario. We practically operated with 3 di fferent scenarios, namely upper, lower and full glycolysis. For upper glycolysis, the NAD(H) balance does not play a role as neither NAD nor NADH is part of the system, but the C, the AXP, and the P i -balance could be used. For the lower and the full glycolysis, NAD(H), C, and AXP-balance could be used while the P i -balance was obsolete due to not measureable P i . For experiments with a constant feed and without perturbations, the balances are constant over time. However, this is not true for other perturbation scenarios, for which the total amount of compounds in the system could change over time. For dynamic feed profiles and i njections, the feed profiles were taken into account to compute the theoretical total compound concentrations in the system in the reactor for each time point of the experiment allowing calculating the balances over time.

Matching of ideal and real mass balances
In the previous section we described the calculation of ideal This correction factor was then used to correct the c2 concentration series This means essentially that the error in the mass balance was distributed equally over the considered concentrations in order to force the concentration balance to close. Certainly the assumption that the observed error is caused by all compounds to the same extent is not necessarily true, but it seemed the least biased way to handle unknown error sources.
In choosing which concentrations were to be corrected according to equation (39) we followed two paths: in the first scenario (leading to the data set "c3") we considered the already performed corrections described in "Matching the observed and set concentrations" for compounds with an initial concentration larger than zero as sufficient and refrained from correcting them again. Hence the remaining error was distributed only over those compounds with initial values equal to zero. In the second scenario (data set "c4") all 67 compounds of the balance were corrected regardless of previous corrections. In consequence the difference lies only in the few substrates provided with the initial feed, e.g. GLC and ATP. Both ways are also shown in the Supplementary Fig. 13. While the usage of the correction factor is straight forward, one aspect had to be considered: the presented balances are only partially interdependent as they share compounds through the phosphate balance. For example FBP is part of the carbon (see Supplementary Equation (34)) as well as of the phosphate balance (Supplementary Equation (37)) and ATP is part of the adenosine (Supplementary Equation (36)) and of the phosphate balance ( Supplementary  Equation (37)). In consequence, if a balance was used to correct certain concentrations this influenced the remaining balances and therefore a certain order in applying the balances had to be chosen. This was only necessary for the upper glycolysis as only there the phosphate balance was used. We applied the balances in the sequence of carbon, phosphate, and finally AXP balance, which in our hands delivered the best results. For the lower and full glycolysis the phosphate balance could not be ev aluated as the phosphate was not part of the measured compounds and t he remaining balances did not overlap in compounds, so no ranking was required.
Injecting the diluted compound mix from the section "Concentration determination in reference to a s tandard in the MS matrix buffer" into the MS showed for the compounds 13BPG and A TP a r etained signal in time reaching steady-state later than the other compounds. At the same time, when the injection was stopped and the system was flushed by pure buffer again a time delayed response giving a signal trail was observed. This effect, reminiscent of chromatography, is most likely the result of binding/unbinding interactions between the PEEK tubing and t hese compounds. While 13BPG concentrations were experimentally very small and t he compound was not used in the feed (commercially unavailable), consequences of the time delay for the widely used ADP/ATP cofactors were of concern. Together with the data processing, which uses initial concentrations as reference where possible, a delay in signal establishment at the beginning could result in overcorrecting this compound, which also could influence the correction of the concentrations of other compounds via the closing of mass balances. At the start of an experiment only ATP was present but no A DP, hence the only compound responsible for a possibly wrongly determined initial concentration was ATP which could be c orrected. Later, when due t o reaction ATP and ADP are present, both compounds may contribute to a measurement error. As we did not know which compound contributed how strongly to the error we decided to distribute the error to both compounds relative to their actual values. If this measure is wrong we would have corrected over time with an increasing error and this could explain why the deviations in ADP/ATP concentrations were only observed in parts of the data. We did not observe this for the other compounds. This is also the reason why the ADP/ATP cofactors pair was not fitted and therefore not used for parameter estimation.

Boundary conditions for parameterization
For the initial parameterization of the reaction model and the definition of parameter boundaries for the parameter estimation we used different sources. V max values were obtained from the known amount of added enzymes, equilibrium constants were calculated from thermodynamic data and intrinsic enzymatic properties such as K M were obtained from databases such as BRENDA 19 and literature research.
Hereby the following considerations applied to the process of parameterization of the reaction system: Generally the upper and l ower boundaries in which parameters were allowed to vary during estimation were set following fixed rules: providing 1 U of activity in the reactor would translate into a V max 0.143 mM min -1 . The activity definition of most enzymes of lower glycolysis by the supplier is valid for 25°C (see Supplementary Table 3) while the reactor temperature was set to 30°C. Using the logarithmic Arrhenius relationship for the temperature dependency of reactions as a first estimate for the temperaturedependent change of enzymatic activity suggested an activity of 0.230 mM min -1 . Therefore, we fixed boundaries for V max values between 60 and 250% of this temperature-adapted reference activity, from 0.143 to 0.572 mM min -1 in terms of activity in the reactor.
Because of frequently broadly varying literature values, we attributed the lowest certainty to parameters describing enzyme-compound interactions, and consequently the estimation of affinities or K M values was constrained only by what were considered reasonable physiological ranges, e.g. between 1 µM and 2 mM. This also allowed testing the ability of our procedures to perform effectively de novo parameter estimation, as in many similar scenarios which do not rely on well-known enzymes such as those derived from glycolysis, little enzyme data will be available from the literature. In contrast, V max -values (as the product of k cat and enzyme concentration), thermodynamic equlibria, and H ill coefficients, were considered to be o f a more certain character and for them the estimation process was constrained to more narrow ranges, e.g. between 0.5 and 2 t imes the initial value for V max and thermodynamic equilibria while Hill coefficients were constrained around the number of reported enzyme subunits, e.g. between 3 and 4 for a tetrameric enzyme. As the substrate FBP and the XAP pool were distinguishable and the stoichiometry of the reaction was pre-defined, the problem of not sufficiently different MS-signals for GAP and D HAP did not limit the parameterization. However, during preparatory estimation runs we observed that the measured concentrations of DHAP and GAP could not be reproduced with the calculated Ald equilibrium even if widely open parameter boundaries were selected for the affinities. One simple explanation for such 69 a systematic error is an error in the determination of the equilibrium constant of the Aldcatalyzed reaction. The value obtained from the formalism used to calculate the equilibrium constants was extremely low, which seemed to be i nconsistent with the naturally required reaction rates through this enzyme, pointing indeed to a systematic error in the computation of the equilibrium constant. This constant is derived from literature adapted to user-specified modifications in the experimental conditions (see section "Calculation of thermodynamic equilibria") and t he measurement of enthalpy changes. However, the determination of the equilibrium constant of Ald is known to be error prone 26 . Therefore, we used the experimentally observed concentrations of FBP and t he XAPs to calculate the XAP/FBP equilibrium as 0.036 mM (instead of 1.2e-3 to 1.7e-4 mM as calculated above). This value was consistent with values theoretically calculated and used in other models (e.g. 27,28 ) and we used it as a constant in the model together with the affinities from literature 29 19 and were able to fit the experimental data within these boundaries.

Parameterization of subsystems in general
The model of the enzymatic reactions has 60 parameters whose values have to be identified by computationally fitting experimental time series. Given the available computation power and the parameter space to be searched (which grows exponentially with each parameter), the estimation of all parameters at the same time was not feasible. Hence the overall problem was split into manageable sub-problems of maximally 4 enz ymes of around 20 parameters as shown in Supplementary Fig. 5. The breakdown of the system was structured according to the availability of substrates, the observability of compounds, and the requirement that the thermodynamic equilibria of reactions should allow a sufficient extent of reaction.
Subsystem A The 4 parameter system was elaborated as an ex ample for structural identification. It contained only Glk, whose few parameters allowed manually investigating results and t est procedures, in the course of which we identified and corrected an inconsistency between the rate equation reported in the literature and our results.

Subsystem B
The system size was increased to upper glycolysis from Glk to Ald, already producing DHAP (22 parameters).

Subsystem C
The 23 parameters choice was influenced by two considerations: firstly GAP was not a good starting point as it was not commercially available in suitable purity, implying that FBP and Ald were needed as part of C to produce GAP. This implied also the production of DHAP which then could be consumed by G3d. Secondly, thermodynamic calculations suggested that the equilibrium of the Gdh reaction is strongly on t he GAP side and so to observe appreciable reaction the product 13BPG must be removed by the action of Pgk.

Subsystem D
The composition of subsystem D (8 parameters) was motivated by the fact that 2PG was again unavailable at suitable purity. Further, the Pgm-catalzyed reaction could not be observed as substrate and product and thus could not be di stinguished by MS. However, 2PG could be transformed by Eno into PEP, which could be observed.

Subsystem E
Both, Pyk-and Ldh-catalyzed reactions are strongly on t he side of the products and were identified together in subsystem E (12 parameters).

Subsystem F
Reconciliation of the various subsystems into one single model required re-estimation of some parameters. The two NADH consuming reactions Ldh and G 3d had to be es timated again together as subsystem F (16 parameters). Pyk was included to provide consistency.
Subsystem G Furthermore, Pfk and Ald were re-estimated as subsystem G (11 parameters) integrating data from reduced and full reaction system to maintain consistency.

Parameterization of subsystem A (Glk only)
We first applied different input functions ( Fig. 2A and Supplementary Fig. 4A, column a, also Supplementary Table 4) to obtain a parameter estimate for purified E. coli Glk using a rate equation derived from literature 1 . As this was the first subsystem we investigated, we looked for each type of experiment at the optical fit quality ( Supplementary Fig. 4A, column b), the distribution of the achieved errors ( Supplementary Fig. 4A, column c), and the distribution of parameters allowing the representation of the experimental data within 1% of the minimal error ( Supplementary Fig. 4A, column d). Clearly, it was impossible to obtain with any of the different input functions a narrow distribution of all the parameters for runs that fulfilled the 1% criterion ( Fig. 2A, upper panel and Supplementary Fig. 4A, column d). This indicated a problem of practical non-identifiability for this combination of experiment and model. Next, we incorporated a product inhibition by ADP into the model rate equation for the Glk reaction (in the form of an apparent affinity ) (see also Supplementary Note 1, Supplementary equations (I) and (II)) and repeated the analysis (Fig. S4B). The experiments had again very good optical qualities of fit, similar to the situation for the enzyme mechanism without inhibition (data not shown). The quality of fit distributions for all types of experiments also remained comparable (Supplementary Fig. 4A and B, column c). The best qualities of fit for type I (constant) and II (step) experiment stayed almost constant irrespective of which mechanism was applied for Glk and similar problems with respect to the affinities remained, suggesting that the non-identifiability problem persisted for these experiments. However, the quality of the fittings for the type III (ramp) experiment was markedly improved and delivered a unique minimum which could be assigned to a single parameter set with very good accuracy ( Fig. 2A, lower panel, Supplementary Fig. 4B, column d, and Supplementary Table  5). Finally, as a control, we conducted an experimental confirmation for an inhibition of Glk by ADP ( Supplementary Fig. 4C) by running constant feed experiments, with GLC or G6P and with GLC and ADP in the feed. We placed 1 U of Glk in the reactor and tracked the concentration of GLC.

Parameterization of subsystem B
Experiments were planned to rationally target different aspects (i.e., parameters) of subsystem B with different types of experiments, specifically constant feed experiments and pulse experiments (Supplementary Table 6). Pulse experiments were expected to allow improved estimation of affinities: Tailored injection of compounds allowed the setting of individual compound concentrations independently from the reaction system by deflecting it from the steady-state of substrate and p roduct concentrations. The dynamic response of (constant) enzyme to changing compound concentrations would then reflect the interactions between enzyme and m etabolites (affinities). In the first pulse experiment (B4 in Supplementary Table 6) we sequentially pulsed metabolites of upper glycolysis after an initial steady-state had established (Supplementary Fig. 6). In experiment B5 we exposed the cooperative Pfk to pulses in known allosteric effectors, namely the inhibitor PEP and the activator ADP with a second role as competitive inhibitor to ATP, and in B6 we screened for effectors of upper glycolysis by sequentially pulsing metabolites of lower glycolysis. However, we did not observe such effects besides the already known inhibition of Pfk by PEP (Supplementary Fig. 6F). The resulting set of parameters for the enzymes from Glk to Ald is summarized in Supplementary Table 7.
Parameterization of subsystem C Please note that for the estimation of parameters for subsystems C to E, we had one experiment (E5) that was used -together with other experiments -for the estimation of all three subsystems. In this experiment, the lower glycolysis was built up s tepwise by subsequently injecting single enzymes once the injection of the previous enzyme had resulted again in a steady state with respect to the metabolite concentrations. This generated an experiment that contained all reactions of lower glycolysis. From this experiment, three individual data sets C3, D3 and E5 were obtained by considering different end times. The subsystem C contains a l arge number of parameters for 4 reactions and is interdependent as the G3d reaction requires NADH formed by the Gdh reaction. Furthermore the predicted equilibria and also the experimental observations showed that the equilibrium of the Gdh reaction could be e xpected to be on t he substrate side and t hus efficient flux through this enzyme required removal of the reaction product 13BPG by Pgk. Thus we performed a two-step approach in experiments and in parameter estimation. First, we used experiments C1 to C3 (Supplementary Table 8, Supplementary Fig. 7) to estimate the parameters for Ald, Gdh and P gk. Next, we used experiments C1 and C4 to estimate the parameters of G3d and performed simultaneously a re-estimation of the parameters of Gdh while keeping the parameters of Ald and Pgk constant. As both enzymes Gdh and G3d act on the NAD/NADH cofactors this allowed (correctly) balancing both reactions and at the same time reducing the numbers of parameters to be es timated. This delivered the parameter set summarized in Supplementary Table 9 which was used for the simulations in Supplementary Fig. 7.

Parameterization of subsystem D
Compounds 3PG and 2PG were only measured as a pool (XPG), impeding the experimental observability of the Pgm reaction. All three compounds were commercially available but 2PG only in low purity of around 70%. As we suspected 3PG to be the most important impurity in the commercial preparation, we abstained in the experiments from using 2PG as we did not want to work with undefined feeds or pulses. Hence the experiments D1 and D 2 (Supplementary Table 10, Supplementary Fig. S8) started with 3PG as feed substrate in a concentration typically observed by us in glycolysis experiments. During the fitting process we observed that the experiments could be fitted exceptionally well within the broad de novo boundaries, such as K M values between 1 µM and 2 m M, but 3 o f the 4 K M 's were at the lower boundary. This outcome seemed to be rather unrealistic as when compared with the literature 19 . Therefore, we selected boundaries more in the physiological range (Supplementary Table 11). Please note that only one of three concentrations in subsystem D could be obs erved directly limiting the identifiability of parameters as e.g. high 2PG concentrations can be compensated by low affinity of an enzyme for 2PG.

Parameterization of subsystem E
Subsystem E contained 12 parameters and was estimated using the experiments E1 to E5 (Supplementary Table 12, Supplementary Fig. 9) to obtain the parameter set of Supplementary Table 13. While the fitted concentrations are generally in good agr eement, the cofactors ADP, ATP, NAD and N ADH qualitatively reproduce the measurements but show slight offsets.

Parameterization of subsystem F
The previous parameter estimations for the subsystems C to E had led to a fully parameterized model of the enzymatic reaction system of lower glycolysis which we had verified for the production of LAC using Ldh (Supplementary Fig. 9). As an additional control we also evaluated it against experimental data of the all lower glycolysis enzymes including G3d for the production of G3P. In this control we observed a mismatch of the model (data not shown), suggesting that the parameters for both enzymes needed t o be estimated in the same experiment. To assure consistency with previous results and data bases we reestimated subsystem E together with G3d resulting in the subsystem F. We used again the data from experiments E1 to E5 and added experiments F1 and F2 (Supplementary Table  14, Supplementary Fig. 10). From these 7 experiments the 13 parameters for the enzymes Pyk, Ldh and G 3d were estimated by fitting the concentrations of DHAP, G3P, PYR and LAC. The estimated parameter values are given in Supplementary Table 15.

Parameterization of subsystem G and verification experiments for full glycolysis
As control and reference for the complete reaction system from GLC to either PYR/G3P or LAC we set up 2 ex periments using 1 U of every enzyme and either relying on Ldh or G3d for the NAD regeneration (experiments G1 and G2 in Supplementary Table 16). Here, we operated for the first time the full reaction system from GLC to the products LAC or G3P, respectively, and we observed a m ajor discrepancy in the FBP concentration (data not shown). As a result, we decided to re-estimate the parameters of Pfk and Ald again by fitting them to the concentrations of GLC, X6P and FBP, using (for consistency reasons) the experiments B1 to B6 and G1 and G 2. Newly estimated were the parameters maximal reaction velocities, affinities and equi libria while other parameters of Pfk, e.g. the Hill coefficient, were kept constant at their previously identified values. The resulting set of parameters for Pfk and Ald is summarized in Supplementary Table 17 and the simulations of the experiments G1 and G2 using this data set are displayed in Supplementary Fig. 11. With this final parameter adaption we had developed a fully parameterized model of the reaction system being able to describe the reactions from Glk to Ldh. The corresponding overall parameter set is summarized in Supplementary Table 18.

Optimized reaction system
The optimization of the reaction system was performed by adapting the activity of single enzymes for a gi ven total amount of activity with regards to one or two compound concentrations. For this the optimizer suggested an activity distribution and simulated the experimental outcome of a constant feed experiment with initial enzyme injection for a chosen experimental length. The finally reached concentration at steady-state was taken as criterion and provided to the optimizer that then suggested another distribution. As starting point for the optimization an equi-activity scenario with every enzyme having the same specific activity was chosen. As parameters for the enzymatic reaction the actually available best parameter set from the stepwise identification of the reaction system was used.

Optimization of upper part of reaction system
With the parameter set of Supplementary Table 7, which had been identified for the upper part of the reaction system, we proceeded to optimize system performance of this section. Specifically, we wanted to optimize the distribution of a gi ven amount of enzyme activities (given total number of units enzyme supplied) to obtain maximum product formation, here specifically for DHAP. DHAP is rapidly produced by Ald in an equilibrium-limited reaction from FBP, which in turn is formed in an i rreversible reaction from F6P. Therefore, in conceptual terms, the relevant optimization problem is the formation of FBP, from which DHAP can simply be formed by addition of sufficient Ald to allow rapid establishment of equilibrium. In consequence we optimized FBP formation by distributing a total of 3 U of enzyme activity over Glk, Pgi and P fk. We optimized for the experimental conditions of constant feed experiments (2 mM GLC and 2 mM ATP) and an experimental time of 60 min, in which the system reached steady-state. The fact that we used non-stoichiometric amounts of ATP (2 moles of ATP are needed to convert 1 mole of GLC to FBP) was of no particular concern, as ultimately ATP was expected to be regenerated in the entire reaction system by including the lower part of glycolysis. Optimization results as shown in Supplementary Table  19 and Figure 3A.

Optimization of lower part of cascade
The lower reaction cascade with NAD regeneration by G3d was optimized for a constant feed experiment (1 mM FBP, 2 mM ADP, 1 mM NAD and 2 mM P i ) and m aximum G3P production as the objective. The total amount of activity for the seven-enzyme system was 7 U. Parameters for Ald, Gdh and Pgk were extracted from Supplementary  Fig. 3B.

Optimization of entire cascade
The whole cascade, with G3d instead of Ldh, was optimized for the concentrations of G3P as end product and P YR for ATP regeneration for a 60 m in constant feed experiment (2 mM GLC, 2 mM ATP, 1 mM NAD and 2 mM P i ). As parameters we used the final parameter set from Supplementary Table 18. Two experiments were designed for 20 and 40 U of total enzymatic activity and results are shown in Supplementary Table 21 and Fig. 3C.

Minimization of NAD concentration in the feed
As a further application of the model, we investigated to what extent the NAD concentration in the feed could be r educed without affecting the performance of the reaction system. Although NAD concentration was already sub-stoichiometric in the scenarios investigated so far, we were interested in further decreasing the NAD requirement as it is by far the most expensive feed compound. We used the model with the final parameter set from Supplementary Table 18 to simulate the effect of decreasing NAD concentrations on the performance of an equi-activity distributed system (1 U of every enzyme) with constant feeds (2 mM GLC, 2 mM ATP and 2 mM P i ), which contained different concentrations of NAD, ranging from 0.1 mM to 2 mM. Simulation results are given in Supplementary Fig. 12. Then, we optimized the distribution of enzymatic activities for the above described reaction conditions and 0 .25 mM NAD for a 60 min continuous feed experiment (concentration as