Integrating Flux Balance Analysis into Kinetic Models to Decipher the Dynamic Metabolism of Shewanella oneidensis MR-1

Shewanella oneidensis MR-1 sequentially utilizes lactate and its waste products (pyruvate and acetate) during batch culture. To decipher MR-1 metabolism, we integrated genome-scale flux balance analysis (FBA) into a multiple-substrate Monod model to perform the dynamic flux balance analysis (dFBA). The dFBA employed a static optimization approach (SOA) by dividing the batch time into small intervals (i.e., ∼400 mini-FBAs), then the Monod model provided time-dependent inflow/outflow fluxes to constrain the mini-FBAs to profile the pseudo-steady-state fluxes in each time interval. The mini-FBAs used a dual-objective function (a weighted combination of “maximizing growth rate” and “minimizing overall flux”) to capture trade-offs between optimal growth and minimal enzyme usage. By fitting the experimental data, a bi-level optimization of dFBA revealed that the optimal weight in the dual-objective function was time-dependent: the objective function was constant in the early growth stage, while the functional weight of minimal enzyme usage increased significantly when lactate became scarce. The dFBA profiled biologically meaningful dynamic MR-1 metabolisms: 1. the oxidative TCA cycle fluxes increased initially and then decreased in the late growth stage; 2. fluxes in the pentose phosphate pathway and gluconeogenesis were stable in the exponential growth period; and 3. the glyoxylate shunt was up-regulated when acetate became the main carbon source for MR-1 growth.


Introduction
Cell metabolisms are highly dependent on environmental conditions, so the metabolic state often shifts during the cultivation period [1,2,3]. Characterizing the transience of metabolic fluxes is important for understanding how cells responded to environmental changes. Bioprocess models (e.g., a Monod-based kinetic model) [4] have been widely applied to predict microbial dynamics, but they cannot directly obtain the intracellular flux distributions. On the other hand, flux balance analysis (FBA) profiles the rates of enzymatic reactions based on stoichiometric mass balance, knowledge of reaction constraints, and measurements of inflow/outflow fluxes [5,6]. As an underdetermined model, FBA requires an objective function (e.g., ''maximizing growth rate'') for flux calculation. However, since cells may show suboptimal metabolism and reprogram their metabolic fluxes under different environmental conditions, the commonly used objective function is insufficient to describe cell physiologies [7,8,9]. Furthermore, FBA assumes steady-state metabolic conditions, and thus is unable to directly analyze the transience of cell metabolism [10,11,12].
This study developed an FBA framework that integrates Monod kinetics and FBA to decipher the dynamic metabolism of MR-1 ( Figure 1). MR-1 is a facultative anaerobic bacterium, which not only plays an important ecological role in carbon cycling and metal reduction, but also has been widely used for in situ bioremediation and microbial fuel cell applications [13,14,15]. MR-1 has a diverse carbon utilization capability and shifts its metabolism during batch cultivation [16]. MR-1 uses lactate for initial growth and produces acetate and pyruvate. In the late growth stage, MR-1 metabolizes less energy-favorable pyruvate and acetate. To describe such kinetic behavior, we used unsegregated Monod equations to simulate cell growth, lactate utilization, and metabolite secretion and reuse. The standard Monod model was incorporated into a genome-scale FBA model, iSO783 [17], to formulate the dynamic FBA (dFBA) framework [11], which enabled quantitative predictions of the MR-1 metabolism.

Monod model
MR-1 growth displayed an apparent lag phase (,7.1 h) in 30 mM lactate medium (0.1% inoculation). By incorporating a time delay function for the lag growth phase, a standard Monod model consisting of four ordinary differential equations was built to describe the extracellular metabolite curves and growth kinetics ( Figure 2). The parameters of the Monod model were estimated by fitting the experimental data. Table 1 indicates that the lactatebased biomass yield was higher than that for either pyruvate or acetate, confirming the preferential utilization of lactate as an energy-favorable carbon substrate for MR-1. Similarly, the lactatebased growth rate (m max ) was much higher than that for either pyruvate or acetate, indicating that lactate was the major carbon substrate for biomass growth at the early growth stage. Table 1 lists the rate coefficients (k pl , k ap and k al ) for waste products (pyruvate and acetate) synthesis and reuse, which indicates that MR-1 quickly consumed lactate, producing significant metabolic overflows to the waste products. Such a strategy illustrates an advantageous ecological niche for MR-1 in competing for favorable carbon sources. Finally, although our standard Monod model reasonably well described MR-1 growth data, its results showed some lack-of-fit with statistical analysis (Table S1). Such a discrepancy was possibly due to the model's simplification and to measurement noises. In this study, the kinetic model represents a compromise between complexity and practical simplicity.

Link kinetic model to FBA
To resolve the flux dynamics, the static optimization approach (SOA) divided the cultivation phase into numerous pseudo-steady states so that a conventional genome-scale MR-1 framework (iSO783, containing 774 reactions and 634 metabolites) was able to calculate the flux distributions [17] in each five-minute time interval. Such dFBA model consisted of ,400 mini-FBAs. To avoid repeated and tedious measurements of biomass and metabolite concentrations for each mini-FBA, we used the Monod model to determine the inflow/outflow fluxes of lactate, acetate, and pyruvate in each time interval. The mini-FBAs could be resolved by an objective function of ''maximizing growth rate'', but this function severely overestimated the actual biomass growth

Author Summary
This study integrates two modeling approaches, a Monod kinetic model and genome-scale flux balance analysis, to analyze the dynamic metabolism of an environmentally important bacterium (S. oneidensis MR-1). The modeling results reveal that MR-1 metabolism is suboptimal for biomass growth, while MR-1 continuously reprograms the intracellular flux distributions in adaption to nutrient conditions. This innovative dFBA framework can be widely used to investigate transient cell metabolisms in response to environmental variations. Furthermore, the dFBA is able to simulate metabolite-labeling dynamics in 13 C-tracer experiments, and thus can serve as a springboard to advanced 13 C-assisted dynamic metabolic flux analysis by using labeled proteinogenic amino acids to improve flux results.  ( Figure 3). To account for the suboptimal metabolic features [7,8,9], we used a dual-objective function in dFBA: a combination of ''maximize growth rate'' and ''minimize overall flux''. By appropriately weighing both objectives, we explored the tradeoffs between optimal cell growth and minimal enzyme usage. Specifically, the Monod model determined the transient growth rate for each time interval, which tuned the weights in the dualobjective functions for the mini-FBAs so that the biomass growth curve simulated by dFBA was in agreement with experimental observations. Figure 3 showed that the optimal dual-objective function in mini-FBAs was time-dependent. In general, these dual-objective functions were invariable before the carbon substrate switched from lactate to acetate/pyruvate. When lactate became scarce, the weight of ''minimizing overall flux'' in the objective function increased significantly, indicating an intracellular reduction of enzyme synthesis and minimization of intracellular fluxes.

Dynamic flux distributions in MR-1
The dynamic flux distributions in MR-1 were calculated using the bi-level optimization ( Figure 4). The carbon flows to extracellular acetate and pyruvate were high when lactate was sufficient (,33% and ,25% of the lactate uptake flux before the carbon substrate switch, respectively). Fluxes into the gluconeogenesis pathway, reductive PP pathway, and ED pathway were mainly for biomass synthesis, and remained approximately constant during the exponential growth phase. In the middle log phase (22,25 hours), when the growth rate reached the maximum, fluxes in the oxidative TCA cycle reached a peak (e.g., ,6 mmol/g DCW/h for succinyl-CoA synthetase) to generate energy and building blocks. When lactate became scarce (25,30 hrs), MR-1 had to utilize its waste metabolites (acetate and pyruvate). During this metabolic shift, most intracellular fluxes started to decrease. In the late log growth phase (30,34 hrs), it was also observed that the glyoxylate shunt was up-regulated compared to TCA cycle fluxes after acetate became the main carbon source for MR-1 growth. The glyoxylate shunt reduced the oxidation of carbon substrate for CO 2 production by diverting the carbon flow into a shorter metabolic route. The glyoxylate shunt activity was further confirmed by in vitro enzyme assays at both the mid-log phase (malate synthase activity was 0.1860.11 mmol/g DCW/min) and the late-log phase (malate synthase activity was 0.3760.17 mmol/g DCW/min).

Simulation of dynamic 13 C-labeling in proteinogenic amino acids
In 13 C-labeled tracer experiments, dFBA can be used to predict the isotopomer dynamics in slow turnover metabolites, such as proteinogenic amino acids. During MR-1 growth with [3-13 C] lactate, the dynamic metabolism led to variations of labeling patterns in intracellular metabolites (biomass precursors) so that the isotopic labeling in proteinogenic amino acids was continuously changing during cell growth [18]. Here, we simulated the time-integrative labeling patterns in proteinogenic amino acids based on fluxes from dFBA. The predicted isotopomer labeling patterns of five proteinogenic amino acids Figure 3. Prediction of growth rates (h 21 ). Blue #: growth rate determined by the Monod model. Red #: dFBA prediction using the objective function (maximization of growth rate). Green %: dFBA prediction using dual-objective functions (maximization of growth rate and minimization of overall flux). Yellow e: the weight of the dual-objective functions that predicted the measured growth rates. Note: the summation of the square of fluxes (gv 2 ) was a very large number (total 774 fluxes), so the magnitude of weight w was small. doi:10.1371/journal.pcbi.1002376.g003 (Ala, Ser, Glu, Asp, and Gly, at t = 24 and 30 h) are illustrated in Figure 5A and Figure S1. Compared with the experimental measurements, the labeling patterns predicted by dFBA are consistent with the measured labeling patterns, but some lack-of-fit still exists. One of the limitations of the FBA model is that the intracellular pathway is treated as unidirectional, so the effect of exchange fluxes on isotopomer data is neglected. Considering that some in vivo reactions could be bi-directional, we imple- mented exchange coefficients for four pathways (e.g., the anaplerotic pathway: pyruvateRmalate) in the model to improve the simulation of 13 C-labeling (Table 2). After introduction of the exchange coefficients, the measured and the simulated isotopomer data for proteinogenic amino acids matched (R 2 = 0.9619, Figure 5B). Discussion dFBA models have been developed to describe the dynamic metabolism of E.coli [10], Saccharomyces cerevisiae [19], Lactococcus lactis [12], and even for a more complicated coculture system of E.coli and Saccharomyces cerevisiae [20]. In this study, we developed dFBA for analyzing metabolic states of S. oneidensis MR-1. The time-dependent inflow/outflow fluxes for dFBA can be constrained by either a Monod model or other empirical models (such as polynomial-fitting to the measurement data [21]). The Monod model is suitable to uncover kinetic properties of a scale-up bioprocess and empowers the dFBA to correlate the bioprocess parameters (such as nutrient concentrations and inhibition coefficients) with intracellular metabolism analysis. The integration of the Monod model and dFBA can decipher and predict cell metabolisms in response to batch fermentation conditions.
To describe biological realities, a physiologically reasonable objective function is important for FBA. For E.coli metabolism, 11 objective functions have been systemically investigated under different cultivation conditions [7]. It turns out that no single objective function can describe metabolic states accurately for all conditions. A recent study of MR-1 indicated that futile cycles may be operational, in which less energetically efficient enzymes are expressed at higher levels than their counterparts and decrease biomass yield [17]. Such suboptimal metabolic features in MR-1 make the conventional objective functions difficult to use in predicting actual cell physiology. To bridge the gap between the in silico simulations and experimental observations, we assigned dualobjective functions to resolve mini-FBAs. Using dual-objective functions, dFBA accurately predicted the elevated flux ratio of the glyoxylate shunt (represented by malate synthase activity) to the oxidative TCA cycle (represented by fumarase activity) when acetate started to be used as the main carbon substrate ( Figure S2). The up-regulation of the glyoxylate shunt and down-regulation of the oxidative TCA pathways were consistent with a previous 13 Cmetabolic flux analysis of MR-1 [18]. In comparison, this transient metabolic shift in the glyoxylate shunt could not be captured by a single objective function, such as maximal biomass growth. Moreover, our dFBA results showed the weight of the two objective functions remained relatively constant when lactate was sufficient. At the early stage of MR-1 growth, such a pseudo-steady state has been experimentally verified by previous isotopomerbased analysis [18]. Under nutrient scarcity conditions, MR-1 metabolism may reduce synthesis and usage of enzymes to achieve a compromise between minimization of general physiological activities and maintenance of essential cellular functions [22].
The dFBA model can also simulate time-dependent isotopomer enrichment in proteinogenic amino acids. In turn, the isotopomer results ( Figure 5) can be used to validate and improve the dFBA model predictions. For example, our dFBA model predicted low fluxes through malic enzyme during the exponential growth because these pathways may reduce biomass production, while the genetic analysis indicates a high functionality of malic enzyme [17]. In the dynamic isotopomer simulations, we found that the fitting of isotopomer labeling patterns was significantly improved by introducing the bi-directional fluxes through the pathway Malate « CO 2 +Pyruvate, while keeping the net flux minimal. Such reversible reactions suggest metabolic flexibility. The activity of malic enzyme was also confirmed by in vitro enzyme assays at both the mid-log phase (malic enzyme activity was 0.9060.18 mmol/g DCW/min) and late-log phase (malic enzyme activity was 1.7360.81 mmol/g DCW/min).
Proteinogenic amino acids are abundant in biomass and can easily be measured by GC-MS. Complementing this instrumental data, 13 C-metabolic flux analysis (MFA) offer analytic insight into the cell metabolisms in fermentation processes [3,23,24]. However, the turnover rate of protein is much slower than that of intracellular metabolites, so 13 C-MFA is useful only for analyzing the steady-state central metabolism. To perform 13 C-MFA of dynamic flux distributions, the fast turnover metabolites have to be extracted and analyzed at multiple time points [1,3], which requires significant sampling efforts and high-sensitivity analytical measurement of low-abundance/unstable metabolites. Moreover, the calculation of dynamic fluxes with isotopomer data formulates an inverse nonlinear optimization problem, which is computationally challenging. Due to insufficient methods for analyzing low abundance metabolites [25], as well as limitations in computational algorithms, dynamic 13 C-MFA cannot resolve the flux distributions in a large-scale metabolic network. To overcome these difficulties, this study illustrates a proof-of-concept approach that exploits the synergy between proteinogenic-amino-acid-based 13 C-MFA and genome-scale dynamic flux balance analysis.
In our dFBA, the Monod model is solved first independently of the FBA. As an alternative approach, we also tested to integrate the kinetic models with FBA (integrative Flux Balance Analysis, iFBA). iFBA simultaneously optimizes the kinetic model parameters and solves the dynamic cell metabolism in MR-1 (Supplementary Text S1). We found that iFBA also requires a dual objective function, a weighted combination of ''maximizing growth rate'' and ''minimizing overall flux'', to improve the model accuracy, similar to the dFBA approach. Such results indicate that it is difficult to use a single objective to describe the flux states under all growth conditions, while the time-dependent trade-off objective functions are effective for analyzing the dynamic suboptimal metabolism.
In conclusion, as in other FBA studies, the dFBA framework proposed in this study links macroscopic bioprocess variables (such as nutrient concentrations) to microscopic intracellular metabolism analysis. It predicts metabolic responses under dynamic culture conditions, and reveals the impact of the kinetic parameters (such as m max ) on intracellular flux distributions. Furthermore, dFBA can identify the objective functions that are possibly used by microorganisms in adaption to environmental variations. Finally, by simulating and comparing the isotopomer labeling patterns of different metabolites, the proposed dFBA framework can potentially improve dynamic flux resolutions by incorporating the isotopomer data from labeled proteinogenic amino acids.

Materials and Methods
Culture conditions, analytical methods, and isotopomer analysis S. oneidensis MR-1 (ATCC 70050) was first grown in LB medium in shake flasks overnight. A 0.1% inoculum volume was then cultured into modified MR-1 defined medium [26] in shake flasks (100 mL culture for each of triplicate experiments, shaken at 150 rpm) at 30uC. The initial carbon source was ,30 mM lactate. The growth curve was monitored by dried biomass weight. The concentrations of lactate and acetate in the medium were measured using enzyme kits (r-Biopharm, Darmstadt, Germany). The concentration of pyruvate in the medium was measured with the enzyme assay developed by Marbach and Weli [27].
To analyze the activity of malate synthase and malic enzyme, samples were taken at early middle log phase (biomass of MR-1 ,0.08 g/L) and late log phase (biomass of MR-1 ,0.23 g/L). The harvested cells were centrifuged and re-suspended in 100 mM Tris buffer. The samples were then ultra-sonicated for 5 min to release the enzymes. Malate synthase activity was gauged based on the reaction of CoASH with DTNB (Acetyl-CoA+glyoxylateRMala-te+CoASH; CoASH+DTNBRCoA-TNB+TNB), as described by Dixon and Kornberg [28]. In general, 20 mL acetyl-CoA (5 mM), 10 mL DTNB (10 mM), 50 mL cell extract, and 500 mL of a solution containing 0.1 M potassium phosphate and 10 mM MgCl 2 were mixed with water. The mixture was then added with 20 mL 100 mM glyoxylate. The difference in absorbance at OD 412 before and after glyoxylate addition reflected the activity of malate synthase, in which one unit DOD 412 corresponded to 70.6 nmol of TNB produced (in a 1 mL reaction solution). Furthermore, the activity of malic enzyme was detected based on increased absorbance at 340 nm due to the reduction of NAD + to NADH [29]. In brief, 400 mL 250 mM Tris HCl, 20 mL 50 mM MnCl 2 , 25 mL 40 mM NH 4 Cl, 100 mL 1 M KCl, 50 mL 20 mM NAD + , 100 mL 100 mM malate, and 50 mL cell extract were mixed with water (1 mL reaction solution). The change in absorbance at OD 340 reflected the activity of malic enzyme.
In the labeling experiment, MR-1 was first grown overnight in the LB medium in shake flasks. A 0.1% inoculum volume was then cultured into 100 mL of modified MR-1 defined medium at 30uC, with the initial carbon source as 30 mM [3-13 C] lactate (purity.98%) purchased from Cambridge Isotope Laboratories, Inc. (Andover, MA). The biomass was harvested at ,24 h (before lactate was depleted) and ,30 h (after the substrate had switched from lactate to waste products). To analyze the labeling pattern of proteinogenic amino acids, we hydrolyzed the biomass with 6 M HCl at 100uC. The isotopic analysis of proteinogenic amino acids was performed using a GC-MS based TBDMS method, as previously described [30,31,32]. Ions of [M-57] + (unfragmented amino acid) were used for the 13 C-simulations [33].

Monod model development
A multiple-substrate Monod model was developed to describe the cell growth, lactate consumption and secretion, and reuse of waste products (acetate and pyruvate).
where X is biomass (g DCW/L); LACT, ACT, and PYR are lactate, acetate, and pyruvate concentrations (mmol/L), respectively; m L , m A , and m P are the specific growth rates (h 21 ) on lactate, acetate, and pyruvate, respectively; k e is the endogenous metabolism rate constant (h 21 ); Y X/L , Y X/A , and Y X/P are the biomass yield coefficients (g DCW/mol substrate) of lactate, acetate, and pyruvate respectively; r P,L and r A,L are the production rates (mmol/L/h) of acetate and pyruvate from lactate, respectively. r A,P is the production rates (mmol/L/h) of acetate from pyruvate. S(t2t L ) is the dimensionless unit-step time delay function (S = 0 when t,t L ; S = 1 when t = t L ) which described the lag phase after inoculation.
The specific cell growth rate was described by Monod equations: where m max,L , m max,A , and m max,P are the maximal growth rates (h 21 ) for fully aerobic growth on lactate, acetate, and pyruvate, respectively; and K s,l , K s,a , and K s,p are Monod constants (mmol/ L) for lactate, acetate, and pyruvate, respectively. The acetate and pyruvate production rates are assumed to be proportional to the biomass [16], as expressed below: r P,L~kpl : LACT : X ð9Þ where k al and k pl are rate constants of acetate and pyruvate production, respectively (LN (hNg DCW) 21 ). k ap is the rate constant of acetate production from pyruvate (LN (hNg DCW) 21 ). The kinetic model (Equations 1,10) contained 14 kinetic parameters. To estimate the model parameters, an ordinary least squares (OLS) method [34] was applied to solve the inverse problem. OLS aimed to minimize the residual sum of the squares (R) between model simulations and experimental measurements, expressed as where g represents four dependent variables simulated by the kinetic model; b represents the vector of the parameters to be estimated; and Y is the vector of the measured value of the dependent variables. Since the scales of the dependent variables were different (e.g., the scale of the biomass measurement was ,1 g/L, while the lactate measurement .10 mmol/L), the direct application of OLS would overemphasize the fitting of dependent variables with large scales. Therefore, we normalized dependent variables by the corresponding maximum concentrations measured in the experiments. The ''ode23'' command in MATLAB (R2009a, Mathworks) was used to solve ODEs, and the ''fmincon'' command was used to find suitable settings of the parameters. Figure S3 is the histogram of normalized model residuals. The standard deviations of the parameters were derived from a bootstrap analysis, in which the experimental measurements were randomly re-sampled 1000 times and the corresponding parameters were simulated with the same parameter estimation approach. The 1000 re-sampling was found to be adequate since the variation of parameters converged to within a desired tolerance of 0.1%. The MATLAB code of parameter estimation in the Monod model was attached in Dataset S1.

Bi-level dFBA study
The growth phase was divided into 408 pseudo-steady-state intervals with instantaneous transitions between the two adjacent intervals [11]. At each pseudo-steady state (,five minutes) [35], a mini-FBA was formulated by a dual-objective function comprised of ''maximizing the growth rate'' and ''minimizing overall flux''; and inflow/outflow fluxes (for lactate, acetate, and pyruvate) derived from the Monod model. The inflow/outflow fluxes were calculated from: dACT dt~v act outflow : X ð13Þ where v lac_inflow , v act_outflow , and v pyr_outflow are transient lactate inflow flux, acetate outflow flux, and pyruvate outflow flux, respectively. At each pseudo steady state, the mini-FBA followed a bi-level optimization framework similar to ObjFind [36]. The internal optimization was an FBA with a combined objective function in which the weighting factor of ''minimizing overall flux'' ranged from zero to one. The difference between the transient growth rate simulated from the FBA and that derived from the Monod model was minimized in the external optimization, by tuning the weighting factor in the combined objective function. The bi-level optimization determined a trade-off between maximizing growth rate and maximizing enzyme efficiency under the specified growth conditions. The bi-level optimization was formulated as: The internal optimization was a typical quadratic programming (QP) problem and was solved using the CPLEX solver in the TOMLAB optimization toolbox (TOMLAB Optimization Inc, Seattle, WA) within MATLAB (R2009a). The external optimization problem (i.e., search for weight) was solved by a grid search. Since the QP problem in this study was convex, the locally searched results were also the global solutions [37]. The MATLAB code of bi-level dFBA was attached in Dataset S1.

Simulation of dynamic 13 C-labeling in proteinogenic amino acids
Our previous 13 C-MFA study of MR-1 showed that the labeling patterns in its proteinogenic amino acids changed during the latestage of batch growth [18]. Using the dynamic flux distributions from dFBA, we could now simulate dynamics of isotopomer labeling patterns in proteinogenic amino acids using the algorithm below.
IDV int t~t 0 j~Ð Since FBA neglects flux reversibility, we implemented exchange coefficients exch to account for the reversibility of four key metabolic pathways (ME2, GHMT, GLYCL, and SUCOAS in iSO783, Table 2). In Equation 16, v exch is the vector of exchanged fluxes in the reversible reactions; v is the vector of the transient fluxes at the t9 interval, simulated from the mini-FBA; p(t) are fluxes to proteinogenic amino acids at each time interval; Dt is the scale of a pseudo steady state (5 min); DX(t) is the biomass produced at each time interval; IMM is the isotopomer mapping matrices describing the carbon atoms transitions from reactants to products in a reaction; IDV(t) are the isotopomer distribution vectors of transient intracellular metabolites at each time interval, which is calculated based on the different labeling patterns of precursors from n pathways; IDV int | t = t9 are the isotopomer distribution vectors in proteinogenic amino acids at the end of the t9 time interval; MDV i,sim and MDV i,exp are the mass distribution vectors for each of the five proteinogenic amino acids, as simulated and as measured by GC-MS, respectively; M is the matrix for converting IDV to MDV; s i is the standard deviation of the GC-MS measurement (error ,2%), which is assumed to be constant in this study. The carbon transitions in the reactions involved were given in Supplementary Text S2, and the SBML file for the metabolic model of Shewanella oneidensis MR-1 was provided in Supplementary Text S3. The ''fmincon'' command in MATLAB was used to find the exchange coefficients, and the ''nlparci'' command in MATLAB (Dataset S1) was used to find the asymptotic confidence intervals of the exchange coefficients. These exchange coefficients significantly improved dFBA simulation of the labeling patterns in proteinogenic amino acids.

Supporting Information
Dataset S1 MATLAB codes used for mathematical modeling. (RAR) Figure S1 Experimental observed and simulated isotopomer labeling patterns [M-57] + in key proteinogenic amino acids. The standard error for GC-MS measurement was ,0.02. Area plot: dynamic isotopomer simulation (case 1: simulation without considering reaction reversibility; case 2: simulation considering reaction reversibility). Bar plot: comparison of experimental data to simulated isotopomer labeling patterns (case 1: without considering reaction reversibility; case 2: considering reaction reversibility). (DOC) Figure S2 Flux ratio of malate synthase (MALS) and fumarase (FUM) in dynamic metabolism of Shewanella oneidensis MR-1. Blue &: time profiles of flux ratio using ''maximizing growth rate'' as the objective function in dFBA; red m: time profiles of flux ratio using dual-objective function in dFBA: a combination of ''maximize growth rate'' and ''minimize overall flux''. The entire growth of MR-1 was divided into three phases. In phase I, lactate was mainly used as the carbon substrate. In phase II, lactate, acetate and pyruvate were used as the carbon substrates. In phase III, acetate was used as the carbon substrate.