 Methodology
 Open Access
 Published:
Efficient estimation of the maximum metabolic productivity of batch systems
Biotechnology for Biofuelsvolume 10, Article number: 28 (2017)
Abstract
Background
Production of chemicals from engineered organisms in a batch culture involves an inherent tradeoff between productivity, yield, and titer. Existing strategies for strain design typically focus on designing mutations that achieve the highest yield possible while maintaining growth viability. While these methods are computationally tractable, an optimum productivity could be achieved by a dynamic strategy in which the intracellular division of resources is permitted to change with time. New methods for the design and implementation of dynamic microbial processes, both computational and experimental, have therefore been explored to maximize productivity. However, solving for the optimal metabolic behavior under the assumption that all fluxes in the cell are free to vary is a challenging numerical task. Previous studies have therefore typically focused on simpler strategies that are more feasible to implement in practice, such as the timedependent control of a single flux or control variable.
Results
This work presents an efficient method for the calculation of a maximum theoretical productivity of a batch culture system using a dynamic optimization framework. The proposed method follows traditional assumptions of dynamic flux balance analysis: first, that internal metabolite fluxes are governed by a pseudosteady state, and secondly that external metabolite fluxes are dynamically bounded. The optimization is achieved via collocation on finite elements, and accounts explicitly for an arbitrary number of flux changes. The method can be further extended to calculate the complete Pareto surface of productivity as a function of yield. We apply this method to succinate production in two engineered microbial hosts, Escherichia coli and Actinobacillus succinogenes, and demonstrate that maximum productivities can be more than doubled under dynamic control regimes.
Conclusions
The maximum theoretical yield is a measure that is well established in the metabolic engineering literature and whose use helps guide strain and pathway selection. We present a robust, efficient method to calculate the maximum theoretical productivity: a metric that will similarly help guide and evaluate the development of dynamic microbial bioconversions. Our results demonstrate that nearly optimal yields and productivities can be achieved with only two discrete flux stages, indicating that neartheoretical productivities might be achievable in practice.
Background
Microbial bioconversion plays a critical role in efforts to enable sustainable production of commodity chemicals from renewable feedstocks. The economic feasibility of the integrated biorefinery concept therefore hinges sharply on the productivity, yield, and titers that can be achieved by a given microbial host [1]. Flux balance modeling has emerged as an important tool in guiding experimental efforts in strain design by predicting the effects of gene knockouts and overexpression on metabolite yields [2]. In developing engineered organisms for optimal performance in batch cultures, a tradeoff is often encountered between the productivity and yield that can be obtained via metabolic interventions [3]. As a result, many existing strategies for strain design such as OptKnock or OptForce involve designing mutations to achieve the highest possible yield while maintaining growth viability [4, 5]. Other strategies, including DySScO, have specifically addressed the importance of productivity through dynamic simulations [6, 7]. Overall, approaches tend to follow the principle of designing static networks with minimum metabolic functionality to achieve desired product yields [8]. While these methods are computationally and experimentally tractable, optimum productivity is likely achieved by a dynamic strategy, in which the partition of resources between biomass and product formation varies with time [9].
Experimental studies have investigated the use of multistage fermentations to increase productivity and/or yield. These techniques range from simple manipulations, such as changing from aerobic to anaerobic conditions [10, 11], to more complex genetic toggle switches [12] or otherwise inducible gene expression [13]. Computational methods have also been developed to aid in optimizing twostage fermentation systems. Gadkar et al. optimized the flux profile of glycerol kinase to maximize the productivity of glycerol from glucose [14]. Similarly, Anesiadis and coworkers proposed an extension to Gadkar’s work, where parameters for a quorumsensing toggle switch are tuned to provide the optimal repression for a target gene to maximize ethanol productivity [15]. By choosing a single reaction target as their control strategy and ensuring that unmodulated fluxes are constrained to maximize a growth objective, these strategies follow similar, wellknown techniques such as MOMA [16] or ROOM [17] in calculating an expected biological response from an experimental perturbation. However, a method to calculate a theoretically optimal global limit for a batch culture’s metabolic productivity does not currently exist.
The calculation of maximum theoretical yield is well established in the metabolic engineering literature, and its use helps guide strain design and pathway selection for static knockout and/or overexpression efforts. Maximum yields are useful even though they are not physically realizable: they reveal the cofactor balancing and pathway split ratios necessary to achieve optimum carbon conservation. In an analogous fashion, an estimate of the maximum theoretical productivity of a batch culture system would prove equally useful for the growing field of designing dynamic metabolic interventions. However, due to the computational burden associated with optimizing over all feasible flux profiles, no generalized framework for such an estimate currently exists.
In this article, we present an efficient method for calculating the maximum theoretical productivity of a batch system based on dynamic optimization [18]. Dynamic optimization is a mathematical approach for solving problems involving a differential equation for a desired endpoint objective, in this case the maximum product production in the minimum time. Rather than guessing values for the metabolic state and repeatedly integrating the system forward in time, the approach involves dividing the time domain into a set of finite elements, and within each, representing the dynamic system by a collection of interpolating polynomials. The polynomials are constrained to be continuous between each finite element, and at each of a set of defined points, constraints are imposed to ensure the derivatives of the interpolating polynomials are consistent with maximum substrate uptake rates and maintenance requirements. In this manner, optimal metabolite profiles are found in a single optimization using a largescale nonlinear programming solver. Dynamic optimization has previously been used in metabolic modeling, including in dynamic flux balance analysis (DFBA) [19] and in calculating optimal control for fedbatch fermentations [20]. However, its widespread usage has been limited by the number of variables that can be simultaneously considered. We therefore remove the explicit mass balance constraints by calculating elementary flux modes [21], and reduce the dimensionality of the optimization problem using yield analysis [22]. Our method accounts explicitly for an arbitrary number of fermentation stages, and allows metabolic fluxes to change continuously over the course of the entire simulation. The method therefore prioritizes finding reasonable upper bounds for the productivities and yields of a given bioprocess over the faithful representation of a cell’s response to genetic perturbation. Since productivity and yield cannot be simultaneously maximized, we demonstrate how this method can be easily extended to calculate the complete productivity vs. yield Pareto frontier. We apply our method to succinic acid production in two microbial hosts: engineered Escherichia coli and native Actinobacillus succinogenes. Succinic acid is a precursor to commodity chemicals in several industries, and a promising intermediate in the development of sustainable chemical production routes [23]. We show that the method can be useful in strain choice by comparing optimal productivity–yield surfaces for two organisms, and further demonstrate that nearly optimal yields and productivities can be achieved with only two discrete flux stages.
Methods
Dynamic flux balance analysis
In flux balance analysis (FBA) models, intracellular metabolic reactions are represented by a stoichiometric matrix \(\mathbf {S}\), such that \(S_{ij}\) represents the quantity of metabolite j produced (or consumed) by reaction i. Fluxes through each reaction are represented by the vector \(\mathbf {v}\). It is assumed that the time scales associated with intracellular metabolite equilibria are much faster than those associated with cell growth or changes to external metabolite concentrations, and therefore that metabolites in the cell can be modeled using a pseudosteadystate approximation, \(\mathbf {Sv} \approx \mathbf {0}\).
FBA models can be extended to consider the dynamics of substrate consumption and biomass formation by including specific bounds on exchange fluxes and allowing the accumulation or depletion of external metabolites. The dynamic system for cell growth considered in this paper is
where \(\mathbf {x}(t)\) represents the concentration of external boundary species, chosen such that \(x_0(t)\) is the concentration of dry cell weight (DCW, in grams), and \(\mathbf {v}(t)\) represents the flux through each reaction in the cell (in \(\mathrm {mmol}\; \mathrm {g}_{\mathrm {DCW}}^{1}\; \mathrm {h}^{1}\)), chosen such that \(v_0(t)\) through \(v_N\) represent the exchange flux for species \(x_0\) through \(x_\text{N}\), respectively. The flux through the first reaction is therefore the specific growth rate (\(v_0(t) \equiv \mu (t)\)), and is in units of h^{−1}. The resulting dynamic flux balance analysis (DFBA) model takes the form of an ordinary differential equation with an embedded linear program, for which efficient methods for the solution of the initial value problem have been developed [24]. However, unlike in dynamic flux balance analysis, in this study we do not impose the optimality of a given cellular objective at each moment in time.
The objective considered in this paper is to maximize productivity of the desired metabolite, \(x_\text{p}\), by finding the optimum intracellular flux profiles \(\mathbf {v}(t)\) and final fermentation time \(t_\text{f}\), subject to steadystate constraints, reaction reversibility, and substrate uptake:
The optimization in Eqs. 1 and 2 therefore takes the form of an optimal control problem, which can be solved by discretizing in time and solving the problem using dynamic optimization. In previous applications of dynamic optimization in DFBA models, intercellular fluxes are either optimized directly [20] or replaced with representative input–output reactions found via pathway analysis [19]. As each dynamic variable requires a number of parameters to fully specify its shape over the course of the fermentation, minimizing the number of modeled fluxes helps ensure the optimizations converge. In this study, we take advantage of the productivity objective to select Paretooptimal pathways, greatly reducing the dimensionality of the optimization.
Calculation of elementary flux modes
To select the optimal metabolic pathways, we calculate elementary flux modes (EFMs) for both of the considered metabolic networks. An EFM is a vector \(\mathbf {r}\) in the right nullspace of \(\mathbf {S}\) (\(\mathbf {S}\mathbf {r} \equiv 0\)), such that no other elementary mode has nonzero entries that are a subset of the nonzero entries of \(r_i\) [21, 25]. EFMs contain the important property that any vector in the right nullspace of \(\mathbf {S}\), i.e., any feasible steadystate flux, can be expressed as a nonnegative combination of the elementary modes \(\mathbf {r}\) [26]. Because of this property and because we are only interested in the effect of an elementary mode on the maximum productivity, we can condense the complete set of EFMs to a convex hull in the projection in which we are interested [22]. We further restrict our analysis to the Pareto front of EFMs which allow optimum productivity, removing inefficient modes without affecting the optimal solution.
Dynamic optimization
We find the flux profiles that achieve the maximum productivity via orthogonal collocation on finite elements, a method for solving endpoint problems involving a dynamic system without an embedded ordinary differential equation (ODE) integrator [18]. In the method, the dynamic system is represented by a series of algebraic constraints that implement an implicit Runge–Kutta method. The resulting sparse nonlinear program is then solved via a largescale nonlinear programming (NLP) solver. Following the approach of orthogonal collocation [27], we represent the state variables \(\mathbf {x}(t)\) from Eq. 1 as a collection of Lagrange interpolating polynomials. A summary of the dimensions and variables used in this method are presented in Tables 1 and 2.
The time domain t of the fermentation is divided into \(N_K\) finite elements with a nondimensionalized internal time \(\tau \in [0, 1]\). Within each finite element, we represent the state vector by a polynomial of degree \(N_\text{D}\) at \(N_\text{D} + 1\) collocation points, denoted \(\tau _n\). We use Gauss–Radau collocation for its stiff decay, and therefore for \(N_\text{D} = 5\), we set \(\tau _k \in \left\{ 0, 0.06, 0.28, 0.58, 0.86, 1.00\right\}\):
To model changes in the flux distribution, the finite elements are allocated between \(N_\text{F}\) distinct fermentation stages. Within each stage, we assume the fractional distribution of elementary modes is held constant, while the total flux through the system is allowed to vary. The dynamic system can thus be calculated by
in which \(\mathbf {R}\) is a matrix of shape \((N_R, N_X)\) that contains the chosen elementary modes; \(\mathbf {y}_l\) contains the fractional distribution of each elementary mode in the given stage; \(a_{jk}\) is the timevarying combined flux through all elementary modes; and \(x_{0jk}\) represents the current biomass concentration at the given collocation point. In addition to changes in fractional EFM distribution, the step size \(h_l\) is also optimized by the nonlinear program and allowed to vary between fermentation modes.
Solution of the problem involves the constrained optimization over the variables found in Table 2 of the productivity objective
Orthogonal collocation imposes a number of constraints during the optimization process (Eqs. 6, 7, 8, 9, 10). The first of these constraints is that state variable profiles must obey system dynamics, accomplished by ensuring that the derivative of the interpolating polynomial defined in Eq. 3 is equal to the dynamics defined in Eq. 4 at each collocation point. Limits are also placed on the specific uptake and secretion rates at each collocation point in accordance with biological measurements, and the sum of the fractional expression of each elementary mode for the given fermentation stage is constrained to unity.
Furthermore, the finite elements are constrained to be continuous:
Bounds are also placed on each of the optimization variables:
Finally, two additional constraints are imposed in order to avoid numerical instability and trivial solutions. The relative change in step size between fermentation modes is constrained to a factor of 10,
and the fermentation must consume at least \(80\%\) of the initial glucose:
Implementation
Metabolic models are specified using cobrapy [28]. EFMs for each metabolic model are calculated via efmtool [26]. All calculations were performed in Python, using the CasADi library [29] to implement the orthogonal collocation method. IPOPT [30] was used to solve the resulting NLP problem.
Results and discussion
We demonstrate the usefulness of the method by calculating the maximum theoretical productivities of succinic acid from glucose for two organisms, E. coli and A. succinogenes. In a developing a DFBA model, we integrate knowledge of the corecarbon metabolic pathways present in the organism together with experimental data on biomass production and substrate uptake rates. An overview of the computational method is provided in Fig. 1.
Corecarbon models of A. succinogenes and E. coli
Actinobacillus succinogenes
A stoichiometric model of corecarbon metabolism in A. succinogenes was created based on genomic evidence and a previous metabolic reconstruction [31]. The model consists of 72 metabolites and 89 mass and chargebalanced reactions. A reaction describing A. succinogenes cell growth was constructed from a number of experimental measurements. First, a biomass yield was calculated from measured values of the specific growth rate and glucose uptake rate, \(0.414\;\mathrm {h}^{1}\) and \(9.5\; \mathrm {mmol}\), respectively [32]. As the models after reduction only map the input–output relationships of glucose to product, the internal details of the biomass function could be approximated without introducing significant error into the method. The stoichiometry of precursor metabolites for biomass synthesis was therefore adapted from an E. coli corecarbon model [33], and scaled to match the known cell composition [34]. Energetic requirements of biomass synthesis, including ATP, NADH, and NADPH demand, were adapted from literature values to match the desired biomass yield [35]. A nongrowthassisted maintenance value of \(4.7\;\mathrm {mmol}\; \mathrm {g}_{\mathrm {DCW}}^{1} \mathrm {h}^{1}\) was also enforced, which was determined from the value estimated by Lin and coworkers [36] by finding the amount of ATP that can be produced from 0.308 g of glucose. Dynamic constraints on substrate uptake were also imposed. Substrate and product inhibition on growth rate in A. succinogenes have previously been quantified with a Han–Levenspiel model [36]. Using the calculated biomass yield of \(0.044\; \mathrm {g}_{\mathrm {DCW}}\)mmol^{−1} glucose^{−1}, the experimentally determined parameters were converted to constraints on glucose uptake:
The parameters used in Eq. 11 are presented in Table 3. Additionally, the lower bound of the flux through nongrowthassociated ATP maintenance was constrained to the estimated value:
No oxygen uptake was allowed to simulate anaerobic growth.
Escherichia coli
A corecarbon model for E. coli central metabolism was taken from Orth et al. [33], which consists of 72 metabolites and 94 reactions. Dynamic models for glucose uptake in E. coli typically assume Michaelis–Menten kinetics, in which high concentrations of glucose ultimately saturate the import mechanisms at their maximum value [37–39]. However, this assumption leads to the erroneous result that maximum productivity is achieved at an infinite initial glucose concentration. Since no suitable literature model for substratelevel growth inhibition from the considered substrates could be found, glucose uptake kinetics were adapted from those used for A. succinogenes. Additional bounds on ATP maintenance and oxygen uptake were adapted from those determined by Feist and coworkers [40]:
Calculation and selection of elementary models for A. succinogenes and E. coli
Elementary flux modes were calculated to reduce the dimensionality of the optimization and to alleviate the necessity of enforcing stoichiometric equilibrium constraints during the dynamic optimization. For A. succinogenes, 4763 EFMs were calculated for growth on glucose and normalized by the glucose uptake rate. After reducing EFMs to the considered boundary species, duplicate EFMs were removed. The boundary species considered for A. succinogenes included biomass, ATP, glucose, succinate, formate, acetate, and pyruvate. Since any feasible flux state can be expressed as a nonnegative combination of elementary flux modes, the flux space in the reduced dimensionality can be spanned without loss of generality by only those EFMs at the vertices of a convex hull. Furthermore, as optimal succinate productivity will be achieved using flux modes on the Pareto front of cell growth, ATP production, and succinate secretion, the reduced set of EFMs were further reduced to a final set of 22 Paretooptimal modes. As EFMs are normalized by glucose consumption, glucose need not be explicitly included in the Pareto surface. For E. coli, the additional aerobic growth modes led to a total of 100,273 EFMs. The boundary species for E. coli were the same as those used for A. succinogenes, with the addition of oxygen and the omission of pyruvate, which was not found to be present in any optimal growth mode. In addition to biomass, ATP, and succinate production, the convex hull of optimal modes in E. coli must also consider oxygen consumption. After reduction, 137 modes were kept for E. coli.
The selected EFMs for each organism are shown in Fig. 2, which plots 2D slices of the 3D (A) and 4D (B) yield surfaces. In A. succinogenes, the chosen elementary modes can span the entirety of the yield space for the considered boundary species. In E. coli, modes that consume high amounts of oxygen while yielding low amounts of ATP are omitted by the algorithm, as they are unlikely to be utilized in the optimal productivity solution.
Orthogonal collocation
Optimum productivities are found via a dynamic optimization framework. To further reduce the dimensionality of the optimization, and to allow the effects of discrete fermentation stages to be explicitly simulated, we divide the fermentation time into a number of discrete stages. Within each stage, fluxes are represented by a total flux profile through all elementary modes along with a fractional breakdown of flux through each elementary mode. An example optimum solution for maximum productivity from a 3stage fermentation in A. succinogenes is shown in Fig. 3. For clarity, this example uses only 4 finite elements per fermentation stage, while all other calculations use a minimum of 20 total finite elements evenly divided between stages. The figure demonstrates how the step sizes \(\mathbf {h}\) are optimized to achieve the desired dynamic profile. Additionally, since EFMs are scaled by glucose uptake, the activity parameter smoothly tracks the maximum allowable uptake rate. ATP production is enforced at each collocation point by requiring that the flux through the ATP boundary reaction is greater than the nongrowthassociated maintenance requirement. This constraint has the effect of requiring the cells to choose pathways that produce extra ATP for cellular demands, instead of funneling all carbon towards biomass or product production. The selected EFMs in Fig. 3b demonstrate that optimal productivity is achieved by prioritizing cell growth early in the fermentation and succinate production in later stages. The elementary mode that represents the maximum theoretical yield for succinate on glucose, 1.71 mol succinate per mol glucose, is used only partially in the last fermentation stage, highlighting that constraints on cell maintenance prohibit the system from achieving the maximum theoretical yield.
Effect of increasing the number of fermentation stages
In addition to reducing the dimensionality of the problem, splitting the fermentation time into discrete stages with independent EFM expression allows a systematic investigation into the effect of stage count on maximum theoretical productivity. Relative flux ratios within each stage are fixed, and therefore represent consistent enzyme expression. Optimal productivities achieved for a varying number of fermentation stages are shown in Fig. 4 for A. succinogenes and E. coli. The righthand side of Fig. 4 plots the dynamic constraints imposed on the solution, including ATP maintenance production rate, maximum glucose uptake rate, and the oxygen uptake rate (for E. coli). These constraints are normalized to specific uptake or production rates, and therefore ATP and oxygen uptake constraints remain constant over the fermentation, while the glucose uptake varies based on external concentrations. In all cases, the solutions closely track the maximum allowable glucose uptake rate. In A. succinogenes, the transition from one to two stages allows the cell to divide its succinate production strategy into separate growth and product formation stages. The addition of subsequent stages permits slightly higher productivities by the gradual transition from growth to product formation. Additionally, as the flux ratios are fixed within each fermentation stage, the fraction of input carbon diverted to ATP production is also fixed. Thus, as lower extracellular glucose concentrations result in lower glucose uptake rates, the cell must dynamically allocate a higher percentage of energy towards meeting its ATP maintenance requirement. Higher numbers of stages therefore allow the cell to more efficiently meet the maintenance constraint, as demonstrated by lower excess ATP produced with additional stages.
In E. coli, maximum succinate production using a single stage is achieved microaerobically. With the addition of a second stage, the optimal fermentation time drops appreciably as succinate production is divided into an aerobic growth phase followed by a microaerobic production phase. The addition of further stages beyond two has little effect on either the succinate productivity or the concentration profiles, serving mainly to keep ATP production closer to the constrained minimum.
Pareto surfaces of yield versus productivity
While separate estimates of maximum theoretical yield and productivity can provide information on the economic feasibility of a bioprocess, in optimized settings one often cannot be increased without decreasing the other. We therefore show how the given method can be easily extended to calculate a full productivity–yield Pareto surface: the envelope in the multiobjective optimization on which productivity cannot be increased without sacrificing yield. The surface is found by first calculating the maximum productivity (\(\mathcal {P}_{\max }\), defined by Eq. 5) via the method described previously.
The nonlinear program is then solved again with yield as the objective,
in which \(\mathbf {x}_\text{g}\) represents the glucose concentration, while holding productivity constrained to a fraction of the maximum productivity,
Computational efficiency in the repeated optimizations is improved by using IPOPT’s warm solve method, which preserves the optimal solution as a starting guess for the next iteration.
Figure 5 plots the resulting yield–productivity surfaces for 50 linearly spaced \(\alpha\) values. By varying the number of allowed fermentation stages, these surfaces reveal the performance gains which can be achieved by allowing additional metabolic flexibility. For A. succinogenes, higher yields can be obtained due to the lower specific ATP maintenance requirement. Additionally, as succinic acid is the naturally predominant fermentation product in A. succinogenes, high yields are obtained even at close to the maximum productivity. In E. coli, notable performance gains are achieved by moving to a 2stage fermentation, as it enables efficient usage of aerobic growth modes. In both cases, moving beyond two distinct flux modes does not substantially increase the yields or productivities that can be achieved.
The calculated optimal surfaces are compared to data on succinate yield and productivity that have been achieved experimentally from both wildtype and engineered organisms, compiled in Tables 4 and 5. For consistency with the assumptions of the method, we limit the data to experimental values in batch cultures on grown on glucose without fedbatch operation or biomass recycling. The measured values of biomass yields and glucose uptake kinetics used in modeling vary from experiment to experiment, and thus the estimated productivity–yield surfaces include a degree of uncertainty. However, the data largely fall within the predicted envelope of singlestage processes, confirming the accuracy of the assumptions made in each model. These results illustrate how the proposed method can be useful for guiding experimental effort: as wildtype A. succinogenes fermentations naturally fall close to the optimal productivity–yield surface, further genetic manipulation of this organism is unlikely to yield significantly improved performance. Similarly with E. coli, improving succinate productivity beyond the values previously recorded will likely require a twostage process involving aerobic growth and anaerobic succinate production.
Conclusions
This study represents a computationally efficient method for determining the maximum theoretical productivity for a batch culture system. As the fields of metabolic engineering and synthetic biology continue to develop techniques for the dynamic manipulation of metabolism, our methodology will enable experimental efforts to be focused on where the greatest improvements can be expected. While this study has focused on finding globally optimal solutions, future implementations might search for optimal strategies using only experimentally tractable EFM selections (i.e., ones that are growthoptimal for 1 or 2 gene knockouts). The method can also easily be generalized for conversion of multiple substrates as long as appropriate experimental data exist to fit detailed models to substrate uptake kinetics. It could also be extended to explicitly optimize final titer instead or in addition to yield and productivity. Overall, this work emphasizes the need for better empirical models of substrate and product exchange rates and growth kinetics in designing dynamic metabolic interventions.
Abbreviations
 DFBA:

dynamic flux balance analysis
 FBA:

flux balance analysis
 EFMs:

elementary flux modes
 ODE:

ordinary differential equation
 NLP:

nonlinear programming
 DCW:

dry cell weight
 ATP:

adenosine triphosphate
 NADH:

nicotinamide adenine dinucleotide
 NADPH:

nicotinamide adenine dinucleotide phosphate
References
 1.
Davis R, Tao L. Tan ECD, Biddy MJ, Beckham GT, Scarlata C, Jacobson J, Cafferty K, Ross J, Lukas J, Knorr D, Schoen P. Process design and economics for the conversion of lignocellulosic biomass to hydrocarbons: diluteacid and enzymatic deconstruction of biomass to sugars and biological conversion of sugars to hydrocarbons. Technical Report NREL/TP510060223. 2013. doi:10.2172/1107470.
 2.
Feist AM, Herrgård MJ, Thiele I, Reed JL, Palsson BØ. Reconstruction of biochemical networks in microorganisms. Nat Rev Microbiol. 2009;7(2):129–43. doi:10.1038/nrmicro1949.
 3.
Holtz WJ, Keasling JD. Engineering static and dynamic control of synthetic pathways. Cell. 2010;140(1):19–23. doi:10.1016/j.cell.2009.12.029.
 4.
Burgard AP, Pharkya P, Maranas CD. Optknock: a bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng. 2003;84(6):647–57. doi:10.1002/bit.10803.
 5.
Ranganathan S, Suthers PF, Maranas CD. Optforce: an optimization procedure for identifying all genetic manipulations leading to targeted overproductions. Plos Comput Biol. 2010;6(4):1000744. doi:10.1371/journal.pcbi.1000744.
 6.
Zhuang K, Yang L, Cluett WR, Mahadevan R. Dynamic strain scanning optimization: an efficient strain design strategy for balanced yield, titer, and productivity. DySScO strategy for strain design. BMC Biotechnol. 2013;13(1):8. doi:10.1186/14726750138.
 7.
Hanly TJ, Henson MA. Dynamic metabolic modeling of a microaerobic yeast coculture: predicting and optimizing ethanol production from glucose/xylose mixtures. Biotechnol Biofuels. 2013;6(1):44. doi:10.1186/17546834644.
 8.
Ruckerbauer DE, Jungreuthmayer C, Zanghellini JUR. Design of optimally constructed metabolic networks of minimal functionality. PLoS ONE. 2014;9(3):e92583. doi:10.1371/journal.pone.0092583.
 9.
Venayak N, Anesiadis N, Cluett WR, Mahadevan R. Engineering metabolism through dynamic control. Curr Opin Biotechnol. 2015;34:142–52. doi:10.1016/j.copbio.2014.12.022.
 10.
Vemuri GN, Eiteman MA, Altman E. Succinate production in dualphase Escherichia coli fermentations depends on the time of transition from aerobic to anaerobic conditions. J Ind Microbiol Biotechnol. 2002;28(6):325–32. doi:10.1038/sj.jim.7000250.
 11.
Andersson C, Hodge D, Berglund KA, Rova U. Effect of different carbon sources on the production of succinic acid using metabolically engineered Escherichia coli. Biotechnol Progr. 2007;23(2):381–8. doi:10.1021/bp060301y.
 12.
Soma Y, Tsuruno K, Wada M, Yokota A, Hanai T. Metabolic flux redirection from a central metabolic pathway toward a synthetic pathway using a metabolic toggle switch. Metab Eng. 2014;23:175–84. doi:10.1016/j.ymben.2014.02.008.
 13.
ValdezCruz NA, Caspeta L, Pérez NO, Ramírez OT, TrujilloRoldán MA. Production of recombinant proteins in E. coli by the heat inducible expression system based on the phage lambda pl and/or pr promoters. Microb Cell Fact. 2010;9(1):1. doi:10.1186/14752859918.
 14.
Gadkar KG, Doyle FJ III, Edwards JS, Mahadevan R. Estimating optimal profiles of genetic alterations using constraintbased models. Biotechnol Bioeng. 2004;89(2):243–51. doi:10.1002/bit.20349.
 15.
Anesiadis N, Cluett WR, Mahadevan R. Dynamic metabolic engineering for increasing bioprocess productivity. Metab Eng. 2008;10(5):255–66. doi:10.1016/j.ymben.2008.06.004.
 16.
Segre D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA. 2002;99(23):15112–7. doi:10.1073/pnas.232349399.
 17.
Shlomi T, Berkman O, Ruppin E. Regulatory on/off minimization of metabolic flux changes after genetic perturbations. Proc Natl Acad Sci USA. 2005;102(21):7695–700. doi:10.1073/pnas.0406346102.
 18.
Biegler LT. An overview of simultaneous strategies for dynamic optimization. Chem Eng Process. 2007;46(11):1043–53. doi:10.1016/j.cep.2006.06.021.
 19.
Mahadevan R, Edwards JS, Doyle FJ III. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys J. 2002;83(3):1331–40. doi:10.1016/S00063495(02)739039.
 20.
Hjersted JL, Henson MA. Optimization of fedbatch Saccharomyces cerevisiae fermentation using dynamic flux balance models. Biotechnol Progr. 2006;22(5):1239–48. doi:10.1021/bp060059v.
 21.
Schuster S, Dandekar T, Fell DA. Detection of elementary flux modes in biochemical networks: a promising tool for pathway analysis and metabolic engineering. Trends Biotechnol. 1999;17(2):53–60. doi:10.1016/s01677799(98)012906.
 22.
Song HS, Ramkrishna D. Reduction of a set of elementary modes using yield analysis. Biotechnol Bioeng. 2009;102(2):554–68. doi:10.1002/bit.22062.
 23.
Song H, Lee SY. Production of succinic acid by bacterial fermentation. Enzyme Microb Tech. 2006;39(3):352–61. doi:10.1016/j.enzmictec.2005.11.043.
 24.
Höffner K, Harwood SM, Barton PI. A reliable simulator for dynamic flux balance analysis. Biotechnol Bioeng. 2012;110(3):792–802. doi:10.1002/bit.24748.
 25.
Jevremovic D, Trinh CT, Srienc F, Boley D. On algebraic properties of extreme pathways in metabolic networks. J Comput Biol. 2010;17(2):107–19. doi:10.1089/cmb.2009.0020.
 26.
Terzer M, Stelling J. Largescale computation of elementary flux modes with bit pattern trees. Method Biochem Anal. 2008;24(19):2229–35. doi:10.1093/bioinformatics/btn401.
 27.
Biegler LT. Nonlinear Programming. doi:10.1137/1.9780898719383.ch10.
 28.
Ebrahim A, Lerman JA, Palsson BØ, Hyduke DR. Cobrapy: constraintsbased reconstruction and analysis for python. BMC Syst Biol. 2013;7(1):74. doi:10.1186/17520509774.
 29.
Andersson J. A generalpurpose software framework for dynamic optimization. Arenberg Doctoral School; 2013.
 30.
Wächter A, Biegler LT. On the implementation of an interiorpoint filter linesearch algorithm for largescale nonlinear programming. Math Progr. 2005;106(1):25–57. doi:10.1007/s101070040559y.
 31.
McKinlay JB, Laivenieks M, Schindler BD, McKinlay AA, Siddaramappa S, Challacombe JF, Lowry SR, Clum A, Lapidus AL, Burkhart KB, Harkins V, Vieille C. A genomic perspective on the potential of Actinobacillus succinogenes for industrial succinate production. BMC Genom. 2010;11(1):680. doi:10.1186/1471216411680.
 32.
McKinlay JB, Zeikus JG, Vieille C. Insights into Actinobacillus succinogenes fermentative metabolism in a chemically defined growth medium. Appl Environ Microb. 2005;71(11):6651–6. doi:10.1128/AEM.71.11.66516656.2005.
 33.
Orth JD, Fleming RMT, Palsson BØ. Reconstruction and use of microbial metabolic networks: the core Escherichia coli metabolic model as an educational guide. EcoSal Plus. 2010; 4. doi:10.1128/ecosalplus.10.2.1.
 34.
der Werf MJV, Guettler MV, Jain MK, Zeikus JG. Environmental and physiological factors affecting the succinate product ratio during carbohydrate fermentation by Actinobacillus sp. 130z. Arch Microbiol. 1997;167(6):332–42. doi:10.1007/s002030050452.
 35.
McKinlay JB, ShacharHill Y, Zeikus JG, Vieille C. Determining Actinobacillus succinogenes metabolic pathways and fluxes by nmr and gcms analyses of 13clabeled metabolic product isotopomers. Metab Eng. 2007;9(2):177–92. doi:10.1016/j.ymben.2006.10.006.
 36.
Lin SKC, Du C, Koutinas A, Wang R, Webb C. Substrate and product inhibition kinetics in succinic acid production by Actinobacillus succinogenes. Biochem Eng J. 2008;41(2):128–35. doi:10.1016/j.bej.2008.03.013.
 37.
Hanly TJ, Urello M, Henson MA. Dynamic flux balance modeling of S. cerevisiae and E. coli cocultures for efficient consumption of glucose/xylose mixtures. Appl Microbiol Biotechnol. 2011;93(6):2529–41. doi:10.1007/s0025301136281.
 38.
Song HS, Ramkrishna D. Cybernetic models based on lumped elementary modes accurately predict strainspecific metabolic function. Biotechnol Bioeng. 2010;108(1):127–40. doi:10.1002/bit.22922.
 39.
Harwood SM, Höffner K, Barton PI. Efficient solution of ordinary differential equations with a parametric lexicographic linear program embedded. Numer Math. 2015. doi:10.1007/s0021101507603.
 40.
Feist AM, Henry CS, Reed JL, Krummenacker M, Joyce AR, Karp PD, Broadbelt LJ, Hatzimanikatis V, Palsson B. A genomescale metabolic reconstruction for Escherichia coli K12 MG1655 that accounts for 1260 ORFs and thermodynamic information. Mol Syst Biol. 2007;3(121):121. doi:10.1038/msb4100155.
 41.
Li J, Jiang M, Chen KQ, Ye Q, Shang LA, Wei P. Effect of redox potential regulation on succinic acid production by Actinobacillus succinogenes. Bioprocess Biosyst Eng. 2010;33(8):911–20. doi:10.1007/s004490100414x.
 42.
Salvachúa D, Mohagheghi A, Smith H, Bradfield MFA, Nicol W, Black BA, Biddy MJ, Dowe N, Beckham GT. Succinic acid production on xyloseenriched biorefinery streams by Actinobacillus succinogenes in batch fermentation. Biotechnol Biofuels. 2016;9(1):121. doi:10.1186/s1306801604251.
 43.
Guettler MV, Jain MK, Rumler D. Method for making succinic acid, bacterial variants for use in the process, and methods for obtaining variants. US5573931 A. 1996.
 44.
Liu YP, Zheng P, Sun ZH, Ni Y, Dong JJ, Zhu LL. Economical succinic acid production from cane molasses by Actinobacillus succinogenes. Bioresour Technol. 2008;99(6):1736–42. doi:10.1016/j.biortech.2007.03.044.
 45.
Xi YL, Chen KQ, Dai WY, Ma JF, Zhang M, Jiang M, Wei P, Ouyang PK. Succinic acid production by Actinobacillus succinogenes nj113 using corn steep liquor powder as nitrogen source. Bioresour Technol. 2013;136:775–9. doi:10.1016/j.biortech.2013.03.107.
 46.
Jantama K, Haupt MJ, Svoronos SA, Zhang X, Moore JC, Shanmugam KT, Ingram LO. Combining metabolic engineering and metabolic evolution to develop nonrecombinant strains of Escherichia coli C that produce succinate and malate. Biotechnol Bioeng. 2008;99(5):1140–53. doi:10.1002/bit.21694.
 47.
Millard CS, Chao YP, Liao JC, Donnelly MI. Enhanced production of succinic acid by overexpression of phosphoenolpyruvate carboxylase in Escherichia coli. Appl Environ Microb. 1996;62(5):1808–10.
 48.
Stols L, Donnelly MI. Production of succinic acid through overexpression of NAD(+)dependent malic enzyme in an Escherichia coli mutant. Appl Environ Microb. 1997;63(7):2695–701.
 49.
Lin H, Bennett GN, San KY. Metabolic engineering of aerobic succinate production systems in Escherichia coli to improve process productivity and achieve the maximum theoretical succinate yield. Metab Eng. 2005;7(2):116–27. doi:10.1016/j.ymben.2004.10.003.
 50.
Sánchez AM, Bennett GN, San KY. Novel pathway engineering design of the anaerobic central metabolic pathway in Escherichia coli to increase succinate yield and productivity. Metab Eng. 2005;7(3):229–39. doi:10.1016/j.ymben.2005.03.001.
Authors’ contributions
PSJ, YJB, and MFC designed the method and wrote the manuscript. PSJ performed all numerical experiments. All authors read and approved the final manuscript.
Acknowledgements
We thank Gregg Beckham for his careful review of the manuscript.
Competing interests
The authors declare that they have no competing interests.
Availability of supporting data
The datasets generated during and/or analyzed during the current study are available in the github repository, https://github.com/pstjohn/productivitymaximization.
Funding
This work was funded by the US Department of Energy’s Bioenergy Technologies Office (DOEBETO), Contract No. DEAC3608GO28308 with the National Renewable Energy Laboratory.
Author information
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Flux balance analysis
 Dynamic optimizations
 Elementary flux modes
 Actinobacillus succinogenes
 Escherichia coli