Strains used and culture conditions
The strains used were Escherichia coli BW25113 (lacIq rrn BT14 ΔlacZWJ16 hsdR514 ΔaraBADAH22 ΔrhaBADLD78) and its ppc (JWK3928), pck (JWK3366) and pykF (JWK1666) mutants. The M9 synthetic medium was used for both batch and continuous cultivations, where it consists of 48 mM Na2HPO4, 22 mM KH2PO4, 10 mM NaCl, 40 mM (NH4)zSO4 and 4 g/l of glucose. The following components were filter-sterilized separatory and then added (per liter of final volume): 1 ml 1 M MgSO4, 1 ml vitamin B1 (1 mgl stock), 1 ml 0.1 mM CaCl2, and 10 ml trace element solution containing (per liter): 0.55g CaCl2 1g FeCl3, 0.1 mg/l MnCl2.4H2O, 0.17 g ZnCl2, 0.043 g CuCl2.2H2O, 0.06 g CoCl2.2H2O, and 0.06 g Na2MoO4.2H2O). Batch and continuous cultures were conducted at 37°C in a working volume of 1 l in a 2l fermentor (M-100, Tokyo Rikakikai, Japan) equipped with pH and temperature sensors. The air flow rate was maintained at 1 l min-1 and 300-400 RPM, which ensured the aerobic condition where the dissolved oxygen concentration was kept at 3-4 ppm. The pH of the culture was controlled at 7.0 ± 0.1 by automatic addition of either 2.0 M HCl or 2.0 M NaOH with a pH controller. The CO2 and O2 concentrations in the off-gas were measured by an off-gas analyzer (DZX-2562, Able Co. Japan).
Measurement for extracellular metabolites concentrations
The 5 ml of samples were taken from the culture broth and rapidly quenched into 15 ml of 60% (v/v) aqueous methanol containing 70 mM HEPES, where quenching solution was kept at -80°C before used. The cells were separated from the culture by centrifugation at 10,000 ×g for 15 min at 0°C. To extract intracellular metabolites from the cell pellet, 500 ml of 50% methanol was added, and the cells were re-suspended by vortexing the mixture. Then 2 ml of 35% of perchloric acid was added, which was pre-cooled on ice, and vortexed again for 10 s. After one freeze-thaw cycle, proteins and cell fragments were removed by centrifugation at 12,000 ×g for 30 min at 0°C. The clear suparnatant was neutralized by adding collected supernatant, and then 895 μm of 5 M K2KO3 was added. Afterward, precipitated perchlorates were removed by another centrifugation step (12,000 ×g at 0°C for 10 min), the clear supernatant was collected, and they were stored as 200 μm aliquots at -20°C for further analysis [20, 21].
The intracellular metabolite concentrations were determined by enzymatic methods, where those were performed mainly as described by [22, 23] and others [3, 20, 21, 24] after some modifications. Enzymatic determinations of metabolites were performed in 50 mM tri ethanolamine buffer (TEA) at pH 7.6. For example, G6P concentration was determined in the presence of 5 mM MgSO4 and 0.48 mM NADP+ after addition of 0.7 U ml-1 of G6PDH. The measurements were made using Luminescence Spectrophotometer (LS55, Perkin Elmer, UK). The chemical determination was made at 37°C .
13C-labeling experiments and sample preparation
The biomass sample was kept on ice for 2-3 minutes, and the sample was centrifuged at 6,000 rpm at 2°C for 15 minutes . The cell pellets were washed three times with 20 mM Tris-HCl at pH 7.6 and suspended in 10 ml of 6 M HCl. The mixture was then hydrolyzed at 105°C for 15 hours in a sealed glass tube. During acid hydrolysis, tryptophan and cysteine were oxidized, and glutamine and asparagine were deaminated. The hydrolysate was evaporated to dryness. The dried material was dissolved in milli-Q water and filtered through a 0.22 μm pore-size filter and evaporated to dryness. About 1.5 ml acetonitrile was added in the dried hydrolyte and incubated at room temperature overnight. After the color of liquid turned a color of pale yellow, it was filtrated through 0.22 μm pore-size filter. The filtrate was then derivatized by N-(tert-butyldimethylsilyl)-N-methyl-trifluoroacetamide (MTBST-FA) (Aldrich, USA) at 110C for 30 minutes and was transferred to GC-MS sample tube for analysis .
13C-labeling experiments were initiated after the culture reached steady state, which was inferred from the stable oxygen and carbon dioxide concentrations in the off-gas and stable OD in the effluent medium for at least twice as long as the residence time. The feed medium with 4 g/l of unlabeled glucose was then replaced by an identical medium containing 0.4 g of [U-13C] glucose, 0.4 g of [1-13C] glucose, and 3.2 g of naturally labeled glucose per liter. Biomass samples for GC-MS analysis were taken after one residence time. Sample preparation, analytical procedures for GC-MS analysis, and flux computation are given elsewhere [16–18, 25].
Metabolic flux analysis
Preparation of biomass hydrolysates and recording of the GC-MS spectra (PerklinElmer, Germany) were made as described previously [25–27]. The program Turbomass Gold (Perklin Elmer,Germany) was used for peak assignment and MS data processing. The main idea is to perform isotopomer balance on carbon atoms in order to track the fate of the labelled carbon atoms from the substrate . Isotopomer balance enables us to determine the isotopomer distributions of the intracellular metabolites in the central metabolic network. Since the isotopomer distribution of amino acids can be inferred from the isotopomer distributions of their corresponding precursors, the GC-MS signals for the amino acids can then be simulated. A set of flux distributions is then determined by minimizing the differences between the experimental and simulated data. For the estimation of GC-MS signals using isotopomer balance, three types of corrections were made to take into account the effects of natural abundance, non-steady state condition, and skewing effect. First, the isotopomer distributions of the input substrate were corrected for natural abundance in the unlabeled glucose and impurity in the labeled glucose. Second, the isotopic steady state condition is normally not attained at the time of harvesting biomass. Thus, the simulated data have to be corrected based on the actual harvesting time by assuming first order washout dynamics . Finally, the simulated GC-MS data were corrected for the natural isotope abundance of O, N, H, Si, S, and C atoms in the derivatizing agent using the correction matrices .
Referring to Figure 1, the dynamic equations may be described based on mass balances as
where [·] denotes the concentration, and the abbreviations for the metabolites are explained in Nomenclature. μ is the specific growth rate and v
are the intracellular fluxes. The superscript "ex" means extra cellular. Note that the last term with μ in each equation denotes the dilution effect due to the increase in cell volume as the cell grows . As shown in Figure 1, GAP and DHAP were lumped together. A sequence of enzymatic reactions by GAPDH, Pgk, Pgm and Eno can be considered to be in equilibrium, and thus it was assumed to be one lumped reaction from GAP to PEP, for simplicity. The overall simulation procedure is given in Additional file 8.
Modeling for the cell growth rate
The estimation for the cell growth rate is by far important. In particular, it is critical for the prediction of the specific gene knockout mutant. The most typical equation for the cell growth rate is the Monod equation:
where S is the substrate concentration and S ≡ [Glcex ] in the present case. The drawback of this equation is that the cell will keep growing at μ
as far as S ≫ K
is usually quite small in the case of the utilization of glucose as the substrate. This does not reflect the fermentation characteristics in practice. Therefore, we introduced another term such as
where X is the cell concentration and X
is the final value of X in the batch culture, and n is the model parameter. X
may be set as S(0)Y
where S (0) is the initial substrate concentration, and Y
is the yield coefficient. Although this equation can be used to simulate the dynamic behavior in batch culture, Eq. (2b) may not be able express the cell growth rate for the mutants. We, therefore, introduce another idea by taking into account the effect of ATP generation on the cell growth rate based on experimental observation.
The ATP production is made either by substrate level phosphorylation and oxidative phosphorylation, where the reducing power of NADH and FADH2 can contribute in generating ATP via oxidative phosphorylation. The pathways involved in electron transfer and oxidative phosphorylation have variable stoichiometry due to the use of different dehydrogenases and cytochromes. Namely, the NADH dehydrogenase NDH-I encoded by nuo gene transports 2 H+/e-, while NDH-II encoded by ndh gene transports 0 H+/e- . Quinol then oxidases cytochromes cyt bO3 and cyt bd to transport 2 H+/e- and 1 H+/e- , respectively . The number of H+ transported into the cell by the membrane bound H+-ATPase to phosphorylate ADP to form ATP has been estimated to be 2-4, with 3 being most likely under aerobic condition . The P/O ratio is a non integer value in practice. Here we consider P/O ratio to be the model parameter to be tuned by fitting the experimental data. Thus the specific ATP production rate may be estimated by the following equation:
and is the specific ATP production rates via oxidative phosphorylation, which may be estimated by
where: (P/O) and (P/O)' are the P/O ratios for NADH and FADH2, respectively.
Figure 2a implies the linear relationship between the specific growth rate and the specific ATP production rate, and then Eq (2b) may be further modified as
(•) is the specific ATP production rate computed from the fluxes, and v
is the adjustable model parameter.
Kinetic equations for the enzyme kinetics
The kinetic equations for the various enzyme-catalyzed reactions considered in the present study are given in Table 1. The specific growth rate was assumed to be expressed by Eq.(1), where Eq.(1a) in Table 1 was used for the cell growth phase as state above, while Eq.(1b) was used for the late growth phase where acetate was used as a carbon source.
The glucose transport via phosphotransferase system (PTS) has been extensively investigated and is catalyzed by a sequence of enzymes such as EI, HPr, EIIAGlc and EIICBGlc. These were encoded by such genes as ptsHI, crr and ptsG. It may be assumed that EIICBGlc can be either phosphorylated first or bind with glucose first, although the phosphorylation of EIIB may be facilitated by binding EIIC to glucose. Based on these assumptions, the kinetic rate equation may be derived for the glucose transport through PTS as given by Eq.(2) in Table 1 , where Glcex and PEP are the substrates for the reaction, and PYR and G6P are the products, which inhibit the reaction rate as those are accumulated.
The phophoglucose isomerase (Pgi) is the first enzyme in glycolysis, and the two Pgi isoenzymes have been lumped into one rate equation . The phosphofructokinase (Pfk) is encoded by pfkA and pfkB in E. coli, where Pfk-1 encoded by pfkA is dominant, accounting for 90% of the total activity . Here, we, therefore, considered only Pfk-1, where it is inhibited by PEP , and its reaction rate may be expressed as given by Eq.(4) of Table 1, where F6P is the substrate for this reaction, and PEP inhibits vPfk as it accumulates . For FDP aldolase reaction, Eq.(5) of Table 1 was considered , where FDP is the substrate and GAP is the product. The minus term in the numerator is for the backward reaction. For GAPDH reaction, Eq.(6) of Table 1  was considered, where GAP is the substrate, and PEP is the product for this reaction. Note that the reactions from GAP to PEP were lumped into GAPDH reaction for simplicity. This reaction rate was assumed to be inhibited by NADH as [NADH]/[NAD] ratio increases. The pyruvate kinase (Pyk) reaction is catalyzed by two isoenzymes, PykI and PykII, encoded by pykF and pykA, respectively. PykI is activated by FDP and inhibited by ATP, whereas PykII is activated by AMP. Here we lumped these together, and the rate equations may be expressed as given by Eq.(7) of Table 1, where PEP and ADP are the substrates for this reaction and inhibited by ATP . Moreover, Eq.(7) indicates that FDP and AMP are the activators of the reaction.
It is quite important to correctly simulate the fluxes at the branch point of PEP, where the enzyme activities of Pyk and Ppc determine the fluxes. Based on the experimental data , it is known that the reaction catalyzed by Ppc exhibits a hyperbolic function with respect to PEP concentration, and the reaction rate is usually low without any activator, where AcCoA is a very potent activator, and FDP alone exhibits no activation, but it produces a strong synergistic activation with AcCoA. Based on these observations, Eq.(8) of Table 1 may be considered . Consider now the reaction in the reverse direction, namely Pck. The kinetics of the Pck in E. coli may be assumed to follow the rapid equilibrium mechanism . Thus, we may express the rate equation as given by Eq.(9) of Table 1, where OAA and ATP are the substrates for this reaction, and the reaction rate is inhibited by its products such as PEP and ADP .
For the rate equation of PDHc, the hyperbolic function (Michaelis-Menten type) was considered with respect to PYR, and it was assumed to be inhibited by NADH as [NADH]/[NAD] ratio increases as given by Eq.(10) of Table 1 . The acetate formation is of important concern in practice for the cell growth and the specific metabolite production in E. coli. The cells undergo a metabolism switch associated with the production and utilization of acetate. When enough glucose is present around the cell, the cell produces and excretes acetate. Once glucose is consumed or becomes low level, the cells utilize acetate instead of excreting it. This acetate-associated metabolic switch occurs just as the cells begin to decelerate growth, or fall into stationary phase in the batch culture. The first pathway for acetate production is Pta-Ack pathway, where the reaction catalyzed by Pta and Ack proceeds through an unstable, high energy, acetyl phosphate (AceP) intermediate. For the Pta reaction, Eq.(11) of Table 1 was considered [40, 41]. The rate equation for Ack was considered as expressed by Eq.(12) of Table 1 . During the stationary phase when glucose was depleted, the cells scavenge for the acetate by using the 2nd pathway of Acs. This reaction catalyzed by Acs proceeds through an enzyme-bound acetyladenylate (Acetyl-AMP) intermediate. Acs has high affinity to acetate while Ack-Pta has low affinity. The rate equation for Acs was assumed to be expressed by the simple Michaelis-Menten equation as given by Eq.(13) of Table 1 .
For citrate synthase reaction, Eq.(14) in Table 1  was considered, where AcCoA and OAA are the substrates, and the reaction is inhibited by NADH. Another important branch point is at isocitrate (ICIT) branch point in the TCA cycle, where ICDH (high affinity to ICIT, Km = 8 mM) and Icl (low affinity to ICIT, Km = 406 mM) compete for ICIT. This branch point is under both gene level regulation by aceBAK operon and the enzyme level regulation by the phosphorylation/dephosphorylation of ICDH. Here, we considered the rate equation for ICDH as given by Eq.(15) of Table 1 , where ICIT is the substrate, and the reaction is inhibited by NAD(P)H as [NAD(P)H]/[NAD(P)] ratio increases and as [αKG] increases. For KGDH reaction, Eq.(18) of Table 1 was considered, where αKG and CoA are the substrates, and the reaction is assumed to be inhibited by NADH as [NADH]/[NAD] ratio increases . For SDH reaction, Eq.(19) of Table 1 was considered, where SUC is the substrate and FUM is the product . For fumarase reaction, Eq.(20) of Table 1 was considered, where FUM is the substrate, and the reaction is inhibited by its product MAL . For MDH reaction, Eq.(21) of Table 1 was considered , where MAL is the substrate, and the reaction is inhibited by its products OAA and NADH as these concentrations increase. For malic enzyme reaction, Eq.(22) of Table 1 was considered .
As for the glyoxylate pathway, Eq. (16) of Table 1 was considered for Icl, ICIT being the substrate of reaction, inhibited by high accumulation of its products SUC and GOX . For malate synthase reaction, Eq.(17) of Table 1 was considered, where GOX and AcCoA are the substrates, and the reaction is inhibited by its product MAL .
For G6PDH and 6PGDH reactions in the PP pathway, Eqs.(23) and (24) of Table 1 were considered, where G6P is the substrate of the former reaction inhibited by NADPH, and 6PG is the substrate of the latter reaction inhibited by NADPH . Other equations for PP pathway such as Rpe, Rpi, TktA, TktB, and Tal were considered in the form of mass-action kinetic, given in Eqs. (25)-(29) in Table 1 .
Once the fluxes were determined, the specific CO2 production rate was assumed to be estimated by
and the total CO2 production rate was computed by the product of net production of CO2 and biomass concentration, i.e. [X].
Moreover, the specific NADPH production rate v
may also be estimated by
It should be noted that the flux determination in the PP pathway is not easy for the enzyme-based kinetic equations. In particular, it is not easy to accurately predict the flux of G6PDH at the branch point of G6P. Figure 2b shows the relationship between the specific growth rate μ and v
computed based on the flux estimates obtained from 13C-labeling experiments [14, 15], indicating the linear relationship between the two as:
Once we identified k
value from Figure 2b, we may compute k
after μ was computed. Then, the flux of G6PDH may be estimated.
As stated in introduction, it is not easy, even if it is not impossible, to express gene-level regulation by mathematical equations. Here, we alternatively introduced several rules for the simulation based on the experimental observations and the roles of global regulators such as given in Additional file 3: Table S1, where the rules 1-3 in Additional file 3: Table S1 came from the experimental observations, while Rule 4 and Rule 5 are known by gene level regulations [47, 48]. Rule 6 came from our previous analysis .
2KG: 2-Keto-D-gluconate; 6PG: 6-Phosphogluconolactone; ACE: Acetate; AcP: Acetyl phosphate; AcCoA: Acetyl-CoA; ASP: Aspartate; ADP: Adenosine diphosphate; ATP: Adenosine-5'-triphosphate; AMP: Adenosine monophosphate; DHAP: Dihydroxyacetone phosphate; E4P: Erythrose 4-phosphate; F6P: Fructose 6-phosphate; FDP: Fructose 1,6-bisphosphate; FUM: Fumarate; G6P: Glucose-6-phosphate; GAP: Glyceraldehyde 3-phosphate; ICIT: Isocitrate; NAD/NADH: Nicotinamide adenine dinucleotide; NADP/NADPH: Nicotinamide adenine dinucleotide phosphate; MAL: Malate; OAA: Oxaloacetate; P: Phosphate; PEP: Phosphoenolpyruvate; PYR: Pyruvate; R5P: Ribulose 5-phosphate; Ru5P: Ribose 5-phosphate; S7P: Sedoheptulose 7-phosphate; SUC: Succinate; X5P: Xylulose 5-phosphate
2KGDH: 2-Keto-D-gluconate Dehydrogenase; Ack: Acetate kinase; Acs: Acetyl coenzyme A synthetase; AcP: Acetyl phosphate; Aldo: Aldolase; CS: Citrate synthase; EI: Cytoplasmatic protein (enzyme1); EII, EIIB, EIIC,: Carbohydrate specific (enzyme II); EIIAGlc,/EIICBGlc: Enzyme for glucose; Eno: Enolase; Fum: Fumarase; G6PDH: Glucose-6-phosphate dehydrogenase; GAPDH: Glyceraldehyde 3-phosphate dehydrogenase; ICDH: Isocitrate dehydrogenase; HPr: Histidine containing protien; Icl: Isocitrate lyase; MDH: Malate dehydrogenase; Mez: Malic enzyme; Ms: Malate synthase; Pck: Phosphoenolpyruvate carboxykinase; PDH: Pyruvate dehydrogenase; Pfk/Pfk-1: Phosphofructokinase-1; Pgi: Phosphoglucose isomerase/Glucosephosphate isomerase; Pgk: Phosphoglycerate kinase; Pgm: Phosphoglycerate mutase; Ppc: PEP carboxylase; Pps: Phosphoenolpyruvate synthase; Pta: Phosphotransacetylase; PTS: Phosphotransferase system; Pyk: Pyruvate kinase; Rpe: Ribulose phosphate 3-epimerase; Rpi: Ribulose 5-phosphate 3-isomerase; SDH: Succinate dehydrogenase; Tal: Transaldolase; TktA: TransketolaseI; TktB: TransketolaseII
aceBAK: operon which encodes the metabolic and regulatory enzymes of the glyoxylate bypass; cra: Catabolite repressor/activator; crr: Catabolite repressor; fadR: Fatty acid metabolism regulator; icdA: Isocitrate dehydrogenase A; iclR: Transcriptional repressor IclR; pckA: Phosphoenolpyruvate carboxykinase gene; pfkA, pfkB: Phosphofructokinase gene; ppsA: Phosphoenolpyruvate synthase gene; ptsG, ptsHI: Pts gene; pykA, pykF: Pyruvate kinase gene