Optimizing Cofactor Specificity of Oxidoreductase Enzymes for the Generation of Microbial Production Strains

Abstract

Central oxidoreductase enzymes (eg, dehydrogenases, reductases) in microbial metabolism often have preferential binding specificity for one of the two major currency metabolites NAD(H) and NADP(H). These enzyme specificities result in a division of the metabolic functionality of the currency metabolites: enzymes reducing NAD⁺ to NADH drive oxidative phosphorylation, and enzymes reducing NADP⁺ to NADPH drive anabolic reactions. In this work, we introduce the computational method OptSwap, which predicts bioprocessing strain designs by identifying optimal modifications of the cofactor binding specificities of oxidoreductase enzyme and complementary reaction knockouts. Using the Escherichia coli genome-scale metabolic model iJO1366, OptSwap predicted eight growth-coupled production designs with significantly greater product yields or substrate-specific productivities than designs predicted with gene knockouts alone. These designs were identified for the production of L-alanine, succinate, acetate, and D-lactate under modeled conditions. Simulations predicted that production of L-alanine and D-lactate can be strongly coupled to growth by knocking out three reactions and swapping the cofactor specificity of one oxidoreductase reaction, while growth coupling was not predicted with four or fewer reaction knockouts under identical conditions. A succinate production design and an acetate production design were predicted to have higher maximum growth rates and higher substrate-specific productivities than designs predicted solely with reaction knockouts. The OptSwap formulation can be readily extended to additional organisms, and the constraints enforcing oxidoreductase specificity swaps can be extended to target other specificity sets of interest.

Introduction

In microorganisms, anabolism and catabolism are precisely controlled by regulation of flux through metabolic enzymes. This regulation of anabolic and catabolic activity is an important element of the adaptive system that enables a microorganism to find the optimal phenotype for growth and reproduction in its environment. The currency metabolites NAD(H) and NADP(H), which carry and transfer reducing equivalents, play unique roles in metabolism. The primary role of the reduced respiratory cofactor NADH is to transfer electrons to oxygen via the electron transport chain, generating the proton gradient that is used for oxidative phosphorylation of adenosine diphosphate (ADP) to adenosine triphosphate (ATP).^1
–3 Concurrently, the reduced cofactor NADPH donates electrons to anabolic reactions and drives biosynthetic pathways in the cell.^2,3 Despite the chemical similarity between NAD(H) and NADP(H), many central dehydrogenase and reductase enzymes in the cell preferentially catalyze reduction or oxidation of a specific carrier.^4

–7 Thus, the functional separation of these electron carriers and the specificity of enzymes to one electron carrier allow the system to direct resources to energy production or anabolism on a whole-cell scale.

Cellular control over the direction of reducing equivalents to NADH or NADPH is facilitated both by the specificity and activity of dehydrogenase reactions and by the activity of transhydrogenase enzymes, which transfer reducing equivalents between the two cofactors. Two transhydrogenases are encoded in the genome and expressed in Escherichia coli. The soluble transhydrogenase encoded by the gene sthA catalyzes the reaction [NAD⁺]+[NADPH]→[NADP⁺]+[NADH]. The membrane-bound transhydrogenase encoded by pntAB couples reduction of NADP⁺ with inward proton translocation, catalyzing the reaction [NADP⁺]+[NADH]+[H⁺ _out]↔[NADPH]+[NAD⁺]+[H⁺ _in]. It has been reported that 35–45% of the NADPH required for biosynthesis is generated by the transhydrogenase encoded by pntAB during aerobic batch growth on glucose.² The large contribution of this transhydrogenase enzyme to NADPH production reflects the importance of transhydrogenases for balancing cofactor generation to meet the needs of the cellular environment.^2,8 Therefore, any strategy to modify cofactor production should consider the role that transhydrogenases play in maintaining homeostasis.

Oxidoreductase specificity for NAD(H) or NADP(H) is a central parameter in determining the direction of cellular resources. Thus, strategies have been developed to modulate the cofactor specificity of dehydrogenase enzymes to increase production of desirable cellular products. Protein engineering has been used to switch the carrier specificity of dehydrogenase enzymes by modifying the amino acid residues of the nucleotide-phosphate binding site. In E. coli, the enzymes dihydrolipoadmide dehydrogenase (a component of the pyruvate dehydrogenase complex and the 2-oxoglutarate dehydrogenase complex) and isopropylmalate dehydrogenase were engineered to prefer NADP(H) over NAD(H) as cofactor, and the enzyme isocitrate dehydrogenase was engineered to prefer NAD(H) over NADP(H).^4,6,9,10 Furthermore, the cofactor specificity of isocitrate dehydrogenase was reversed by protein engineering, and a competition study was performed in which the modified strain outcompeted wild-type when grown on glucose.¹¹ Despite these efforts, dehydrogenase enzymes modified by protein engineering have not been exploited for metabolic engineering purposes.

An alternative strategy to engineer cofactor specificity is to replace native dehydrogenase reactions in central carbon metabolism with heterologous dehydrogenases that have specificity for the opposite cofactor. These non-native enzymes may be more efficient than enzymes produced using protein engineering because the heterologous enzymes have benefited from evolutionary optimization. To this end, the NAD(H)-dependent glyceraldehyde-3-phosphate dehydrogenase in E. coli was replaced with the NADP(H)-dependent glyceraldehyde-3-phosphate from Clostridium acetobutylicum to increase lycopene yield.¹² In Saccharomyces cerevisiae, the native NAD(H)-dependent glyceraldehyde-3-phosphate dehydrogenase was replaced with the NADP(H)-dependent enzyme from Kluyveromyces lactis to increase fermentation of D-xylose to ethanol.¹³ Thus, researchers have exploited homologous swaps to boost production phenotypes.

The scientific interest shown in modifying electron carrier availability and the yield improvements seen in these modified strains suggest that a computational approach to model oxidoreductase specificity modifications could guide future experimental work and ultimately affect metabolic engineering. Determining which oxidoreductase specificity modifications are likely to have the greatest impact on the system and how oxidoreductase specificity changes can be paired with reaction knockouts most effectively can be answered with an in silico modeling procedure.

As a tool for investigating metabolic networks in silico, constraint-based reconstruction and analysis (COBRA) methods have been shown to enable accurate prediction of bacterial behavior under many conditions.^14,15 COBRA methods utilize genome-scale models (GEMS), which are built by pairing the reconstructed metabolic network of an organism with governing constraints based on physico-chemical conservations (ie, known reaction stoichiometries), spatial limitations, and environmental parameters. Governing constraints are formulated as a set of linear inequalities that enclose the solution space available to metabolic fluxes.¹⁶ Solution spaces can be examined by optimizing for objectives using flux balance analysis (FBA) and other linear programming methods.¹⁷ The optimal solutions predicted by FBA can match in vivo behavior for cases in which the metabolic network of the cell is optimized for the same objective (eg, growth). A powerful approach to achieve in vivo optimality is adaptive laboratory evolution (ALE). ALE optimizes the genotype of the organism, and the result is often a match between the observed growth phenotypes and the model predictions.^18

–21 Growth-coupled strains—strains in which growth is directly coupled to the production of a given molecule—have drawn attention because these strains are predicted to produce high yields of the target molecule after ALE selecting for cell growth.²² Computational algorithms that identify growth-coupled designs have proliferated and include OptKnock, RobustKnock, and OptGene.^17,23
–25 Thus, in silico COBRA methods coupled with ALE constitute a strategy for rational engineering of high-yield production strains. A COBRA method for optimizing cofactor specificity has not been reported. Previous investigations have explored the importance of cofactor balancing for microbial production strains, but de novo strain design has not been explored using non-biased cofactor optimization.^8,13,26

In this work, we present OptSwap, an in silico, constraint-based modeling technique for generating strategies to optimize the production of native cellular compounds by modifying the electron carrier specificity of oxidoreductase reactions in the metabolic network. A mixed-integer linear programming (MILP) method is presented that optimizes growth-coupled product yield by pairing oxidoreductase specificity swaps with reaction knockouts. Utilizing OptSwap, we predict novel, growth-coupled designs for the production of valuable compounds by E. coli, and we compare these designs to solutions that are predicted utilizing reaction knockouts alone.

Materials and Methods

Modeling and Computational Tools

The iJO1366 metabolic construction of E. coli K-12 MG1655 was used for all simulations in this work.²⁷ As described previously, formate-hydrogen lyase (FHL) and the oxidative stress reactions catalase (CAT), cytosolic superoxide dismutase (SPODM), and periplasmic superoxide dismutase (SPODMpp) were constrained to zero.²⁷ The pyruvate:ferredoxin oxidoreductase (POR5) reaction was made irreversible, as supported by biochemical data.²⁸

FBA, parsimonious flux balance analysis (pFBA), flux variability analysis (FVA), and RobustKnock were implemented in MATLAB (The MathWorks Inc., Natick, MA) as described in the literature.^24,29
–31 All simulations were performed using MATLAB and the COBRA Toolbox software packages with TOMLAB/CPLEX (Tomlab Optimization Inc., San Diego, CA) and Gurobi (Gurobi Optimization, Inc., Houston, TX) LP/MILP solvers.³²

Substrate uptake rates for the solitary carbon substrates in each simulation were constrained to a maximum uptake rate of 20 mmol/gDW/h. For aerobic simulations, the oxygen uptake rate was set to a maximum of 20 mmol/gDW/h. These values were chosen based on experimental observations of aerobic and anaerobic growth of E. coli.^33,34

Model Reduction and Selection of Reaction Set for Knockouts

The reaction set available for knockout was restricted based on a previously reported method for model reduction and target reaction selection.²² After setting the bounds for the primary carbon source and oxygen exchange, we performed the following procedure, as previously reported.²² First, to generate a “reduced model,” reactions that could not be utilized under the conditions of the simulation were removed, and upper and lower bounds of all reactions were set to maximum and minimum obtainable values, as determined by FVA. Non-gene associated reactions, spontaneous reactions, transport reactions, reactions in peripheral metabolic pathways, and reactions acting on high-carbon containing molecules were then removed from the knockout reaction set. Finally, sets of coupled reactions were identified by examining the null space of the stoichiometric matrix, and only one reaction in each correlated set was considered for knockout.³⁵

With glucose as the substrate, the reactions available for knockout during optimization numbered 217 anaerobically and 243 aerobically. With D-xylose as the substrate, the reactions available for knockout numbered 216 anaerobically and 238 aerobically (Supplementary Table S1 ; Supplementary Data are available online at www.liebertpub.com/ind). Both the proton-translocating transhydrogenase reaction THD2pp (pntAB) and the energy-independent transhydrogenase reaction NADTRHD (sthA) were constrained to zero for the transhydrogenase-knockout design ΔpntAB ΔsthA.

Selection of Reaction Set for Cofactor Specificity Swaps

The set of oxidoreductase reactions available for modification in the OptSwap procedure was determined by finding the oxidoreductase reactions in the high-flux backbone.³⁶ In the metabolic model iJO1366, all reactions that utilize NAD(H) or NADP(H) as a cofactor were located. After pFBA optimization for flux, the reactions were sorted by flux magnitude through the biomass objective function under conditions of aerobic and anaerobic growth on glucose and D-xylose minimal media. The reactions with highest flux under each set of conditions were selected. D-Lactate dehydrogenase, malic enzymes, and L-1,2-propanediol oxidoreductase were added to the set based on interest in the literature.^37
–39 Malate dehydrogenase was removed from the pool of oxidoreductase enzymes that can be swapped with OptSwap because the NADPH–specific malate dehydrogenase allowed non-physiological loops to form in the flux solutions. Under anaerobic conditions, NADH:oxidoreductase I was removed from the pool of reactions during model reduction because the reaction cannot carry flux during simulations of anaerobic growth. Thus, 22 oxidoreductase enzymes were chosen under aerobic conditions and 21 oxidoreductase enzymes under anaerobic conditions (Table 1).^40

–43

Table 1.

Oxidoreductase Reactions Targeted for Analysis with OptSwap^a

					pFBA, RELATIVE FLUX SCALED TO CARBON UPTAKE
KEY OXIDOREDUCTASE ENZYMES IN E. coli	GENE SYMBOL	MODEL REACTION ABBREVIATION	NATIVE ELECTRON CARRIER IN E. coli	PAST STUDIES EXPLORING ALTERNATIVE COFACTOR USES	GLUCOSE ANAEROBIC	GLUCOSE AEROBIC	D-XYLOSE ANAEROBIC	D-XYLOSE AEROBIC
Glyceraldehyde-3-phosphate dehydrogenase	gapA	GAPD	NADH	Non-native enzyme^12,13	1.89	1.68	1.59	1.38
Acetaldehyde dehydrogenase	adhE OR mhpF	ACALD	NADH		0.93	0.0	0.78	0.0
Ethanol dehydrogenase	adhP OR adhE	ALCD2x	NADH		0.93	0.0	0.78	0.0
Glutamate dehydrogenase	gdhA	GLUDy	NADPH	Non-native enzyme⁴⁰	0.21	0.62	0.15	0.55
Glucose-6-phosphate dehydrogenase	zwf	G6PDH2r	NADPH		0.15	0.45	0.11	0.40
6-Phosphogluconate dehydrogenase	gnd	GND	NADPH		0.15	0.45	0.11	0.40
FAD reductase	fre	FADRx	NADH		0.06	0.16	0.04	0.14
Phosphoglycerate dehydrogenase	serA	PGCD	NADH		0.04	0.12	0.03	0.11
Isocitrate dehydrogenase	icd	ICDHyr	NADPH	In silico structural analysis^10,41 Protein engineering^5,11,38,42	0.03	0.08	0.02	0.07
Aspartate-semialdehyde dehydrogenase	asd	ASAD	NADPH		0.03	0.08	0.02	0.07
Methylene tetrahydrofolate dehydrogenase	folD	MTHFD	NADPH		0.02	0.06	0.02	0.06
Acetohydroxy acid isomeroreductase(2-Acetolactate)	ilvC	KARA1	NADPH		0.02	0.06	0.02	0.06
Homoserine hydrogenase	metL OR thrA	HSDy	NADPH		0.02	0.05	0.01	0.05
3-Isopropylmalate dehydrogenase	leuB	IPMD	NADH	Protein engineering^4,41	0.01	0.03	0.01	0.03
Shikimate dehydrogenase	aroE	SHK3Dr	NADPH		0.01	0.03	0.01	0.02
Dihydrodipicolinate reductase	dapB	DHDPRy	NADPH		0.01	0.03	0.01	0.02
NADH:ubiquinone oxidoreductase I	nuoF, nuoA–C, nuoE, nuoG–N	NADH16pp	NADH	Directed evolution⁴³	0.0	1.98	0.0	1.98
Pyruvate dehydrogenase	lpd , aceE, aceF	PDH	NADH	Protein engineering^6,9	0.0	0.37	0.0	0.66
L-1,2-Propanediol oxidoreductase	fucO	LCARR	NADH		0.0	0.0	0.0	0.0
D-Lactate dehydrogenase	ldhA	LDH_D	NADH		0.0	0.0	0.0	0.0
Malic enzyme (NADH)	maeA	ME1	NADH		0.0	0.0	0.0	0.0
Malic enzyme (NADPH)	maeB	ME2	NADPH		0.0	0.0	0.0	0.0

Reactions were chosen by identifying oxidoreductase reactions with highest predicted flux in the genome-scale model during aerobic and anaerobic growth on glucose and D-xylose minimal media. Relative flux values are reported from pFBA simulation optimizing flux through the biomass objective function under these conditions. The enzymes L-1,2-propanediol oxidoreductase, D-lactate dehydrogenase, and malic enzymes were also included in the analysis because interest has been shown in these enzymes in the literature. Malate dehydrogenase was removed from the reaction pool because thermodynamically infeasible loops were created in the presence of the swapped enzyme and, under aerobic conditions, NADH:ubiquinone oxidoreductase I was removed from the model during the model reduction procedure.

MILP Formulation

OptSwap is a bi-level MILP problem based on RobustKnock with new constraints that enforce swaps of the cofactor specificity of oxidoreductase reactions (Fig. 1).²⁴ The following procedure was used to incorporate these constraints into the RobustKnock problem. First, for each oxidoreductase enzyme in the pool of the reactions that can be “swapped,” a reaction with opposite specificity—either NAD(H) or NADP(H)—was added to the model. New Boolean decision variables were utilized. The variable s_d represents the on/off state of the native oxidoreductase reactions, and t_d represents the on/off state of the swapped reactions, where a value of 0 means the reaction is off and a value of 1 means the reaction is on. D is the set of oxidoreductase reaction pairs (native and swapped): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} s_d \in \{ 0 , 1 \} \quad \forall d \in D \qquad\qquad\qquad{\textit{Equation} \ 1} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}t_d \in \{ 0 , 1 \} \quad \forall d \in D \qquad\qquad\qquad{\textit{Equation} \ 2}\end{align*} \end{document}

Fig. 1.

The OptSwap formulation for optimizing cofactor specificity of major metabolic enzymes [NAD(H) vs. NADP(H)] coupled to reaction knockouts. Constraints are added to the MILP problem to enforce swapping the cofactor specificity of reactions catalyzed by oxidoreductase enzymes (A). Production envelopes show the ability to growth couple a product of interest, which is not possible with native cofactor specificity (B). G6PDH: glucose-6-phosphate dehydrogenase; G6P: glucose-6-phosphate; 6PGL: 6-phophogluconolactone.

These variables are present in addition to the RobustKnock Boolean variable y_e _, which represents the on/off state of all reactions that can be knocked out in the model. E is the set of reactions that can be knocked out. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}y_e \in \{ 0 , 1 \} \quad \forall e \in E \qquad\qquad\qquad{\textit{Equation} \ 3}\end{align*} \end{document}

Second, a constraint was added to the outer problem that requires either the native or the swapped reaction to be knocked out for each oxidoreductase reaction pair in D. This is the constraint that forces an oxidoreductase swap: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}s_e + t_d = 1 \ \,\forall d \in D \qquad\qquad\qquad{\textit{Equation} \ 4}\end{align*} \end{document}

Third, a constraint was added to the outer problem to limit the number of swaps to be less than or equal to the parameter L: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\mathop\sum_{d \in D} ( 1 - s_d ) \leq L \qquad\qquad\qquad{\textit{Equation} \ 5}\end{align*} \end{document}

The knockout limitation constraint from RobustKnock ensures that the number of reaction knockouts is less than or equal to the parameter K: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\mathop\sum_{e \in E} ( 1 - y_e ) \leq K \qquad\qquad\qquad{\textit{Equation} \ 6}\end{align*} \end{document}

Fourth, a constraint was added to the outer problem to limit the number of interventions (oxidoreductase swaps and reaction knockouts) to be less than or equal to the parameter X: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\mathop\sum_{e \in E} ( 1 - y_e ) + \mathop\sum_{d \in D} ( 1 - s_d ) \leq X \qquad\qquad\qquad {\textit{Equation} \ 7}\end{align*} \end{document}

These three constraints on the number of knockouts and swaps can be included or excluded from the problem according to the desired simulation scenario. If the oxidoreductase set D is empty, then OptSwap reduces to the RobustKnock problem.

Fifth, a set of constraints was added to limit flux to zero for any native or swapped oxidoreductase reaction whose corresponding decision variable is equal to zero. The function x maps the set of oxidoreductase reaction pairs, D, to the corresponding fluxes and bounds for native oxidoreductase reactions in the model, and the function y maps the set of oxidoreductase reaction pairs to the corresponding fluxes and bounds for the non-native, swapped oxidoreductase reactions. Thus, when s_d =0, the flux v_x(d) is constrained to zero: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}s_d LB_{x ( d ) } \leq v_{x ( d ) } \leq s_d UB_{x( d ) } \quad \forall d \in D \qquad\qquad\qquad{\textit{Equation} \ 8}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}t_d LB_{y ( d ) } \leq v_{y ( d ) } \leq t_d UB_{y( d ) } \quad \forall d \in D \qquad\qquad\qquad{\textit{Equation} \ 9}\end{align*} \end{document}

These are present in addition to the knockout constraint from RobustKnock. The function z maps the set of reactions that can be knocked out, E, to the corresponding fluxes: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}y_e LB_{z ( e ) } \leq v_{z ( e ) } \leq y_e UB_{z( e ) } \quad \forall e \in E \qquad\qquad\qquad{\textit{Equation} \ 10}\end{align*} \end{document}

As an illustration of the variables s and t, consider the case in which s_d =0. Then, by Equation 4, t_d =1. Equation 9 reduces to LB_y(d)<=v_y(d)<=UB_y(d) , so flux through the swapped oxidoreductase reaction is constrained only by the lower and upper bounds—the reaction is “on.” For the same case, Equation 8 reduces to 0<=v_x(d)<=0, so flux through the native oxidoreductase reaction is constrained to zero—the reaction is “off.”

Thus, the final formulation of the OptSwap problem can be stated as follows: J is the set of all reaction fluxes, and I is the set of all metabolites in the model:

max min v _chemical

s.t. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\left[ \begin{matrix} \max v_{{ \rm biomass}}\hfill \\ { \rm s.t.}\hfill \\ \mathop\sum \limits_{j \in J} S_{ij}v_j= 0\hfill \forall i \in I\hfill \\ LB_j \leq v_j \leq UB_j\hfill \forall j \in J\hfill \\ s_d LB_{x ( d ) } \leq v_{x ( d ) } \leq s_d UB_{x ( d ) }\hfill \forall d \in D\hfill \\ t_d LB_{y( d ) } \leq v_{y ( d ) } \leq t_d UB_{y ( d ) }\hfill \forall d \in D\hfill \\ y_e LB_{z ( e ) } \leq v_{z ( e ) } \leq y_e UB_{z ( e ) }\hfill \forall e \in E\hfill\end{matrix} \right] \qquad\qquad\qquad Equation \ \it11\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} s_d \in \{ 0 , 1 \} \quad \forall d \in D \\ t_d \in \{ 0 , 1 \} \quad \forall d \in D \\ y_e \in \{ 0 , 1 \} \quad \forall e \in E \\ s_d + t_d = 1 \quad \forall d \in D \\ \mathop\sum \limits_{d \in D} \left( 1 - s_d \right) \le L \\ \mathop\sum \limits_{e \in E} \left( 1 - y_e \right) \le K \\ \mathop\sum \limits_{e \in E} \left( 1 - {y_e} \right) + \mathop\sum \limits_{d \in D} \left( 1 - s_d \right) \leq X\end{align*} \end{document}

To solve this bi-level max-min problem, we implemented the techniques that were used to simplify and solve RobustKnock.²⁴ The function dual_embed from the RobustKnock implementation was used to generate the dual of the inner problem and linearize the bilinear terms (Equations 8 –10). By the strong duality theory of linear programming, the outer problem and the dual of the inner problem could be integrated.²⁴ Finally, the minimization problem in the resulting min-max problem was converted to a maximization problem using the dual_embed function. After this second conversion and integration, the OptSwap problem is a max-max optimization that can be solved with standard MILP solvers (Supplementary Table S2).

A slight variation of the OptSwap MILP was also investigated. Equation 4 can be converted to an inequality so that oxidoreductase reactions can be swapped (s_d=0, t_d=1) or knocked out (s_d=t_d=0): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}s_d + t_d \leq 1 \quad \forall d \in D \qquad\qquad\qquad Equation \ \it12\end{align*} \end{document}

However, this more general formulation proved to be more computationally intensive and presented numerical challenges that we could not resolve with the TOMLAB/CPLEX solver. Results with this formulation were limited to run for 12 hours, and many simulations did not solve to optimality (Supplementary Table S4).

Results and Discussion

The in silico metabolic optimization procedure OptSwap was used to predict production strains that couple cellular growth to production of targeted products using a combination of oxidoreductase specificity swaps and reaction knockouts. To optimize for maximal growth-coupled yield, a MILP problem was developed based on the RobustKnock problem.²⁴ RobustKnock is a bi-level optimization problem similar to OptKnock, except that OptKnock maximizes production of a target molecule subject to maximizing biomass production, whereas the RobustKnock problem maximizes the minimum production of a target product subject to maximizing growth rate. This procedure ensures that non-unique solutions are not produced—a concern that was first addressed by tilting the objective function prior to the appearance of RobustKnock.²² OptSwap contains additional constraints so that both reaction knockouts and oxidoreductase swaps can be utilized to identify a growth-coupled design.

OptSwap was demonstrated using the iJO1366 genome-scale metabolic model of E. coli. ²⁷ The pool of oxidoreductase enzymes that could be swapped in the analysis was determined by identifying central, high-flux-carrying reactions. Where reaction knockouts were considered for use with OptSwap, a reaction selection workflow was utilized, as previously reported.²² All swap designs were based on the transhydrogenase-deficient ΔpntAB ΔsthA genotype to ensure that unconstrained transhydrogenase activity in the model did not counteract the effect of swapping oxidoreductase specificity.

Growth-coupled designs predicted by OptSwap with both oxidoreductase specificity swaps and reaction knockouts were compared with solutions from RobustKnock in which only reaction knockouts were considered. The designs were categorized by the number of interventions; an intervention is a reaction knockout or a dehydrogenase swap. Thirteen high-value products were investigated based on previous screening criteria for industrially significant compounds produced natively in E. coli.²² For each product, simulations were performed under conditions of glucose and D-xylose minimal media for both aerobic and anaerobic growth (Supplementary Table S3).

For the conditions considered in these simulations, OptSwap produced 10 designs in which coupling oxidoreductase specificity swaps with reaction knockouts resulted in superior predicted production rate compared to just reaction knockouts (Table 2). The eight designs in which OptSwap predicted significantly stronger growth coupling or higher substrate-specific productivity (yield×growth rate) were investigated in more detail (Fig. 2). OptSwap predicted a strongly growth-coupled design for L-alanine production that secretes no coproducts and generates L-alanine with high yield under anaerobic conditions on glucose minimal media. The design requires only three knockouts and one swap of the native enzyme (Fig. 3; Supplementary Fig. S6); it utilizes an NAD(H)-specific glutamate dehydrogenase in place of the native NADP(H)-specific reaction (ghdA). NAD(H)-specific glutamate dehydrogenase transfers electrons from NADH to glutamate. Subsequently, L-alanine transaminase (alaA or alaC) produces L-alanine and 2-oxoglutarate from glutamate and pyruvate. Knockouts of D-lactate dehydrogenase (ldhA) and ethanol dehydrogenase (adhP or adhE) prevent D-lactate and ethanol production. With the new NADH-specific glutamate dehydrogenase pathway in place, L-alanine secretion is energetically favorable to succinate secretion. Knockout of acetyl-CoA acetyltransferase (atoB) is also predicted to be necessary for growth coupling (Fig. 4).

Fig. 2.

Calculated production envelopes for OptSwap designs predicted to have significantly higher optimal production rate or substrate-specific productivity than designs with just reaction knockouts. The wild-type (solid grey line) and the 0pntAB 0sthA genotype (dashed black line) are shown. All designs were found by first knocking out transhydrogenases (0pntAB 0sthA). The optimized phenotypes for OptSwap (solid red line) and RobustKnock (dashed orange line) are compared. Each phenotype required four interventions (one intervention being either a reaction knockout or a oxidoreductase swap). The star (*) indicates that the same design was identified in both cases. The orange X indicates that no growth-coupled design was identified by RobustKnock under these conditions. OptSwap predicts growth-coupled designs for producing L-alanine under all four conditions and D-lactate aerobically on D-xylose substrate; growth coupling is not predicted for these products with just four or fewer reaction knockouts under identical conditions. The OptSwap designs for succinate and acetate are predicted to have higher yield than designs with just reaction knockouts. The designs for succinate production and acetate production on glucose substrate are also predicted to have higher substrate-specific productivities than the designs predicted by RobustKnock.

Fig. 3.

Network diagrams showing the shift in reduced cofactor usage between wild-type and the L-alanine production design. Simulated wild-type flux distribution during anaerobic fermentation in E. coli. The fluxes shown are unique solutions calculated using pFBA, optimizing flux through the biomass objective function (A). L-alanine production design; red X's indicate reaction deletions, and the orange box indicates the oxidoreductase swap (B). Reactions are shown with arrows pointing in the direction of flux. Knockouts of ethanol dehydrogenase (adhP or adhE), acetyl-CoA acetyltransferase (atoB), and lactate dehydrogenase (ldhA) reactions prevent ethanol and D-lactate from being produced. With the NAD(H)-specific glutamate dehydrogenase (gdhA_nadh) in place, L-alanine is predicted to be the energetically favorable final fermentation product.

Fig. 4.

The predicted effect of combinatorial interventions on the design for anaerobic production of L-alanine on glucose minimal media. Glutamate dehydrogenase must be swapped for high yields of L-alanine to be produced at high growth rates; dashed blue envelope shows production with the native glutamate dehydrogenase. With fewer than three reaction knockouts and the swapped glutamate dehydrogenase, high yield of L-alanine is possible, but L-alanine production is not growth-coupled. The final set of three reactions knockouts and one dehydrogenase swap causes a strong coupling between L-alanine production and cell growth. LDH_D: D-lactate dehydrogenase; ACACT1r: acetyl-CoA acetyltransferase; ALCD2x: ethanol dehydrogenase; GLUDy: glutamate dehydrogenase.

Table 2.

Calculated Maximum Optimal Production Rate for Designs Selected by OptSwap and Designs Selected by RobustKnock^a

MAXIMUM OPTIMAL PRODUCTION RATE, OptSwap/RobustKnock (mmol/gDW/h)
	ANAEROBIC								AEROBIC
	GLUCOSE				D-XYLOSE				GLUCOSE				D-XYLOSE
	INTERVENTIONS				INTERVENTIONS				INTERVENTIONS				INTERVENTIONS
	1	2	3	4	1	2	3	4	1	2	3	4	1	2	3	4
Ethanol	22.3/22.3	34.7/34.7	36.4/36.4	37.6/37.6	23.4/23.4	31.9/31.9	32.3/32.3	32.4/32.4	2.6/2.6	18.0/18.0	24.8/24.8	24.9/25.0	0.0/0.0	13.8/13.8	22.4/22.4	23.3/22.6
Formate	36.5/36.5	36.9/36.9	37.2/40.3	34.6/40.7	30.6/30.6	31.5/31.5	32.5/33.5	29.5/33.8	29.8/29.8	36.3/38.5	39.1/39.1	39.1/39.1	25.2/25.2	29.1/29.1	32.5/32.5	34.8/34.8
Succinate	0.5/0.5	12.1/12.1	17.8/17.8	25.4/25.4	14.9/14.9	20.0/20.0	21.9/21.9	21.9/21.9	1.5/1.5	1.6/1.4	5.1/5.1	10.1/9.8	1.3/1.3	6.1/6.1	12.7/12.7	18.0/18.0
Acetate	17.7/17.7	16.1/24.8	25.1/26.8	31.6/29.8	19.7/19.7	19.7/20.3	19.7/21.4	24.5/22.1	25.1/25.1	32.2/36.6	33.3/37.3	42.4/40.0	20.8/20.8	32.6/32.6	38.4/38.4	39.0/39.0
D-Lactate	0.0/0.0	34.7/34.7	34.9/35.1	37.6/37.6	0.0/0.0	31.9/31.9	31.9/32.0	32.2/32.2	0.0/0.0	13.5/13.5	17.1/17.1	25.0/25.0	0.0/0.0	9.0/9.0	12.6/12.6	19.2/0.0
2-Oxoglutarate	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/7.0
L-Alanine	0.0/0.0	0.0/0.0	0.0/0.0	33.5/0.0	0.0/0.0	0.0/0.0	0.0/0.0	31.7/0.0	0.0/0.0	0.0/0.0	0.0/0.0	13.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	8.5/0.0
Glycerol	0.0/0.0	0.0/0.0	0.0/0.0	0.1/0.1	0.0/0.0	0.0/0.0	5.5/5.5	10.0/10.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	3.3/3.3
L-Serine	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0
Pyruvate	0.0/0.0	0.0/0.0	0.0/12.0	0.0/21.0	0.0/0.0	0.0/0.0	0.0/11.4	0.0/16.4	0.0/0.0	0.0/0.0	0.0/19.3	17.0/23.2	0.0/0.0	0.0/0.0	0.0/13.8	11.6/17.7
Fumarate	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0
L-Malate	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0
L-Glutamate	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0	0.0/0.0

Results are arranged by number of interventions, in which one intervention is either one knockout or one oxidoreductase swap. For OptSwap, any combination of swaps and knockouts was allowed. The cases in which OptSwap predicted growth-coupled designs with higher maximum optimal production rate are in bold, and the cases in which the RobustKnock designs are predicted to have higher maximum optimal production rate are in italics. Ten OptSwap designs were identified with higher maximum optimal production rate than designs with just reaction knockouts.

For anaerobic production of L-alanine from D-xylose and for aerobic production of L-alanine from glucose and D-xylose, slightly different designs were predicted. Phenotypes for aerobic production of L-alanine are predicted to have lower yield and higher maximum growth rate than the anaerobic designs (Fig. 2). In all four L-alanine production designs, swapping the cofactor specificity of glutamate dehydrogenase results in predicted phenotypes with strong growth coupling, while growth coupling is not predicted in any designs with four or fewer reaction knockouts. However, previously reported simulations suggest that growth-coupling of L-alanine production may be possible with five reaction knockouts in E. coli.²²

OptSwap predicted that anaerobic production of acetate can be increased with four interventions (Fig. 2). The acetate designs rely on replacing the native NAD(H)-dependent ethanol dehydrogenase (adhP or adhE) with a NADP(H)-dependent ethanol dehydrogenase. On both glucose and D-xylose substrates, maximum optimal acetate production is predicted to be higher with the OptSwap design than with just four or fewer reaction knockouts. Furthermore, simulation predicted that acetate production on glucose substrate with the OptSwap design can be obtained at higher growth rates, and thus with higher substrate-specific productivity.

Simulations predicted that succinate can be produced aerobically with glucose as the substrate using three knockouts and one dehydrogenase swap. The OptSwap design for succinate production utilizes a NAD(H)-dependent glutamate dehydrogenase in place of the native enzyme, coupled to three reaction knockouts. The result is a small increase in the predicted succinate production rate at maximum growth rate, but a large increase in the predicted substrate-specific productivity, so this design is when substrate-specific productivity is desirable.

The OptSwap design for aerobic production of D-lactate on D-xylose substrate is predicted to be strongly growth-coupled, whereas RobustKnock does not predict any growth-coupling for D-lactate with four or fewer interventions under identical conditions. The OptSwap design swaps the native NAD(H)-dependent ethanol dehydrogenase (adhP or adhE) with a NADP(H)-dependent ethanol dehydrogenase. Three knockouts were predicted: acetate kinase (tdcD or ackA or purT), ATP synthase (atpA–H), and pyruvate formate lyase (pflA and pflB). ATP synthase is not predicted to be essential by the metabolic model for growth on either of the substrates considered in these simulations (glucose amd D-xylose). Experimental evidence indicates that ATP synthase knockout strains grow slowly on glucose in vivo.⁴⁴

A number of cases exist in which RobustKnock predicted a design with greater optimal yield than the OptSwap solution (Table 2). In these cases, the RobustKnock solution contains a knocked out oxidoreductase reaction. In the formulation of OptSwap presented here, oxidoreductase reactions cannot be knocked out—Equation 4 enforces swaps and disallows knockouts. An alternative OptSwap formulation was considered in which oxidoreductase reaction knockouts are allowed. The equality constraint (Equation 4) is replaced with an inequality (Equation 12) that allows swaps or knockouts. However, the solution space of the more-general MILP problem caused challenges for the MILP solver used for this work. Numerical imprecision caused non-real solutions, and solution times increased to more than 12 hours. However, when the new problems solved to optimality (for a number of substrate/product combinations tested), they identified the same solutions found with RobustKnock (Supplementary Table S4).

Conclusions

This study presents a computational method to predict optimal cofactor specificity of oxidoreductase reactions in the genome-scale metabolic model. The designs identified by OptSwap are non-intuitive solutions to produce desirable cellular products—solutions that are not possible with just reaction knockouts. Growth-coupled designs for L-alanine and D-lactate were identified by OptSwap using four interventions, in which reaction knockouts could not be used to identify growth-coupled production phenotypes under identical conditions. The anaerobic production designs for L-alanine utilize a NAD(H)-specific glutamate dehydrogenase coupled to reaction knockouts that limit the pathways available for anaerobic fermentation in order to produce very high yields of L-alanine at maximum growth rate. For succinate and acetate, yield improvements were smaller, but large increases in substrate-specific productivity were predicted with OptSwap when compared to designs with just reaction knockouts.

The unique anaerobic L-alanine production design reported here is simpler than a previously reported design requiring seven knockouts and one gene insertion.³⁷ In that study, the native D-lactate dehydrogenase of E. coli was replaced with alanine dehydrogenase from Geobacillus stearothermophilus. To ensure growth coupling of L-alanine, the authors also disabled fermentation pathways for ethanol, acetate, formate, D-lactate, L-lactate, and succinate. By utilizing the predictive power of constraint-based modeling, we predicted a design with high yield and fewer genetic manipulations.

The OptSwap in silico design strategies can be broadly applied to any product in the genome-scale model. A previous analysis demonstrated that growth coupling in genome-scale models is most easily achieved for a class of products related to the native fermentation pathways in E. coli.²² However, techniques like OptSwap may expand the scope of products that can be growth coupled by exploring a broader set of solutions. The knockout and swap designs predicted by OptSwap are candidates for experimental validation and can be built with current technologies. More complex designs than those presented here have been successfully implemented in E. coli, including the strain for production of L-alanine.³⁷

The main limitation of OptSwap is that it is more computationally demanding than previous methods, including RobustKnock and OptKnock. The greater complexity of the solution space means that finding an optimal solution is more difficult. MILP problems are nondeterministic polynomial time (NP)-complete, and solution time depends on the structure of the problem and the solving method. The solution times increased with the parameters K, L, and X (Supplementary Table S5), and adopting a more general problem increases solution times further. While the simpler formulation of OptSwap reported here is computationally tractable, the more general problem allowing oxidoreductase knockouts could not be solved in many cases within 12 hours. CPLEX is a sophisticated MILP solver, and tweaking solver parameters did yield improvements in performance. OptSwap can be solved by choosing the less general formulation utilized for these results and restricting interventions (L or X) to be four or less. Improvements in MILP solvers and computational power will reduce these challenges in the future, and even more complex MILP problems will be solvable.

OptSwap adds a new level of complexity to the COBRA methods for identification of useful production phenotypes in microorganisms.¹⁷ The constraints placed on the presence/absence of oxidoreductase reactions with specificity for different cofactors allow direct investigation and manipulation of cofactor pools that are central to directing cellular resources in microorganisms. Changing cofactor specificity alters the metabolic solution space. While reaction knockouts always decrease the size of the solution space, changing cofactor specificity can have more complex effects and may even increase the solution space in the direction of an objective function.

OptSwap can be readily implemented in metabolic models of other organisms. For example, the experimental findings with a dehydrogenase swap in S. cerevisiae could be further investigated with the genome-scale metabolic model of that organism.^{13
,
45} The constraints that enforce dehydrogenase specificity swaps can also be extended to target other metabolic specificity sets (eg, nucleoside triphosphates) and can be incorporated into other in silico design strategies.

Footnotes

Acknowledgments

We would like to thank Daniel Zielinski, Ali Ebrahim, and Bernhard Palsson for their assistance and insight on the project. We would like to acknowledge the Novo Nordisk Foundation for providing funding for the project.

Author Disclosure Statement

No competing financial interests exist.

References

Russell

, Cook

. Energetics of bacterial growth: Balance of anabolic and catabolic reactions. Microbiol Rev, 1995; 59,1.

Sauer

, Canonaco

, Heri

et al. The soluble and membrane-bound transhydrogenases UdhA and PntAB have divergent functions in NADPH metabolism of Escherichia coli. J Biol Chem, 2004; 279,8:6613–6619.

Gottschalk

. Bacterial Metabolism, 2nd. New York: Springer-Verlag, 1986.

Lunzer

, Miller

, Felsheim

, Dean

. The biochemical architecture of an ancient adaptive landscape. Science, 2005; 310,5747:499–501.

Rodríguez-Arnedo

, Camacho

, Llorca

, Bonete

M-J

. Complete reversal of coenzyme specificity of isocitrate dehydrogenase from Haloferax volcanii. Protein J, 2005; 24,5:259–266.

Bocanegra

, Scrutton

, Perham

. Creation of an NADP-dependent pyruvate dehydrogenase multienzyme complex by protein engineering. Biochemistry, 1993; 32,11:2737–2740.

Kim

, Lee

, Nam

, Kim

. Alteration of reducing powers in an isogenic phosphoglucose isomerase (pgi)-disrupted Escherichia coli expressing NAD(P)-dependent malic enzymes and NADP-dependent glyceraldehyde 3-phosphate dehydrogenase. Lett Appl Microbiol, 2011; 52,5:433–440.

Fuhrer

, Sauer

. Different biochemical mechanisms ensure network-wide balancing of reducing equivalents in microbial metabolism. J Bacteriol, 2009; 191,7:2112–2121.

Guest

, Abdel-Hamid

, Auger

et al. Physiological effects of replacing the PDH complex of E. coli by genetically engineered variants or by pyruvate oxidase. Gordon

, Mulchand

. Thiamine: Catalytic Mechanisms in Normal and Disease States. New York: CRC Press, 2003; 387–407.

10.

Hurley

, Chen

, Dean

. Determinants of cofactor specificity in isocitrate dehydrogenase: Structure of an engineering NADP⁺→NAD⁺ specificity-reversal mutant. Biochemistry, 1996; 2960,95:5670–5678.

11.

Zhu

, Golding

, Dean

. The selective cause of an ancient adaptation. Science, 2005; 307,5713:1279–1282.

12.

Martínez

, Zhu

, Lin

et al. Replacing Escherichia coli NAD-dependent glyceraldehyde 3-phosphate dehydrogenase (GAPDH) with a NADP-dependent enzyme from Clostridium acetobutylicum facilitates NADPH dependent pathways. Metab Eng, 2008; 10,6:352–359.

13.

Verho

, Londesborough

, Penttilä

, Richard

. Engineering redox cofactor regeneration for improved pentose fermentation in Saccharomyces cerevisiae. Appl Environ Microbiol, 2003; 69,10:5892.

14.

Feist

, Palsson

BØ

. The growing scope of applications of genome-scale metabolic reconstructions using Escherichia coli. Nat Biotechnol, 2008; 26,6:659–667.

15.

McCloskey

, Palsson

BØ

, Feist

. Basic and applied uses of genome-scale metabolic network reconstructions of Escherichia coli. Mol Syst Biol, 2013; 9,661:1–15.

16.

Price

, Reed

, Palsson

BØ

. Genome-scale models of microbial cells: Evaluating the consequences of constraints. Nat Rev Microbiol, 2004; 2,11:886–897.

17.

Lewis

, Nagarajan

, Palsson

BØ

. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nat Rev Microbiol, 2012; 10,4:291–305.

18.

Ibarra

, Edwards

, Palsson

BØ

. Escherichia coli K-12 undergoes adaptive evolution to achieve in silico predicted optimal growth. Nature, 2002; 420:186–189.

19.

Fong

, Nanchen

, Palsson

BØ

, Sauer

. Latent pathway activation and increased pathway capacity enable Escherichia coli adaptation to loss of key metabolic enzyme. J Biol Chem, 2006; 281,12:8024–8033.

20.

Fong

, Palsson

BØ

. Metabolic gene-deletion strains of Escherichia coli evolve to computationally predicted growth phenotypes. Nat Genet, 2004; 36,10:1056–1058.

21.

Fong

, Burgard

, Herring

et al. In silico design and adaptive evolution of Escherichia coli for production of lactic acid. Biotechnol Bioeng, 2005; 91,5:643–648.

22.

Feist

, Zielinski

, Orth

et al. Model-driven evaluation of the production potential for growth-coupled products of Escherichia coli. Metab Eng, 2010; 12,3:173–186.

23.

Burgard

, Pharkya

, Maranas

. Optknock: A bilevel programming framework for identifying gene knockout strategies for microbial strain optimization. Biotechnol Bioeng, 2003; 84,6:647–657.

24.

Tepper

, Shlomi

. Predicting metabolic engineering knockout strategies for chemical production: accounting for competing pathways. Bioinformatics, 2010; 26,4:536–543.

25.

Patil

, Rocha

, Förster

, Nielsen

. Evolutionary programming as a platform for in silico metabolic engineering. BMC Bioinformatics, 2005; 6:308.

26.

Ghosh

, Zhao

, Price

. Genome-scale consequences of cofactor balancing in engineered pentose utilization pathways in Saccharomyces cerevisiae. PLoS One, 2011; 6,11:e27316.

27.

Orth

, Conrad

, Na

et al. A comprehensive genome-scale reconstruction of Escherichia coli metabolism—2011. Mol Syst Biol, 2011; 7,535:535.

28.

Blaschkowski

, Neuer

, Ludwig-Festl

, Knappe

. Routes of flavodoxin and ferredoxin reduction in Escherichia coli. Eur J Biochem, 1982; 123,3:563–569.

29.

Kauffman

, Prakash

, Edwards

. Advances in flux balance analysis. Curr Opin Biotechnol, 2003; 14,5:491–496.

30.

Lewis

, Hixson

, Conrad

et al. Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. Mol Syst Biol, 2010; 6:390.

31.

Mahadevan

, Schilling

. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng, 2003; 5,4:264–276.

32.

Becker

S a

, Feist

, Mo

et al. Quantitative prediction of cellular metabolism with constraint-based models: The COBRA Toolbox. Nat Protoc, 2007; 2,3:727–738.

33.

Varma

, Boesch

, Palsson

. Stoichiometric interpretation of Escherichia coli glucose catabolism under various oxygenation rates. Appl Environ Microbiol, 1993; 59,8:2465–2473.

34.

Varma

, Palsson

BØ

. Metabolic flux balancing: Basic concepts, scientific and practical use. Nat Biotechnol, 1994; 12:994–998.

35.

Palsson

BØ

. Systems Biology: Properties of Reconstructed Networks. Cambridge, UK: Cambridge University Press, 2006.

36.

Almaas

, Kovács

, Vicsek

et al. Global organization of metabolic fluxes in the bacterium Escherichia coli. Nature, 2004; 427,6977:839–843.

37.

Zhang

, Jantama

, Moore

et al. Production of L -alanine by metabolically engineered Escherichia coli. Appl Microbiol Biotechnol, 2007; 77,2:355–366.

38.

Wang

, Wang

, Zheng

et al. Biochemical properties and physiological roles of NADP-dependent malic enzyme in Escherichia coli. J Microbiol, 2011; 49,5:797–802.

39.

Stols

, Donnelly

. Production of succinic acid through overexpression of NAD(+)-dependent malic enzyme in an Escherichia coli mutant. Appl Environ Microbiol, 1997; 63,7:2695–2701.

40.

Marx

, Eikmanns

, Sahm

et al. Response of the central metabolism in Corynebacterium glutamicum to the use of an NADH-dependent glutamate dehydrogenase. Metab Eng, 1999; 1,1:35–48.

41.

Kalinina

, Gelfand

. Amino acid residues that determine functional specificity of NADP-and NAD-dependent isocitrate and isopropylmalate dehydrogenases. Proteins, 2006; 64,4:1001–1009.

42.

Yaoi

, Miyazaki

, Oshima

et al. Conversion of the coenzyme specificity of isocitrate dehydrogenase by module replacement. J Biochem, 1996; 119,5:1014–1018.

43.

Auriol

, Bestel-Corre

, Claude

J-B

et al. Stress-induced evolution of Escherichia coli points to original concepts in respiratory cofactor selectivity. Proc Natl Acad Sci, 2011; 108,4:1278–1283.

44.

Baba

, Ara

, Hasegawa

et al. Construction of Escherichia coli K-12 in-frame, single-gene knockout mutants: The Keio collection. Mol Syst Biol, 2006; 2. 10.1038/msb4100050

45.

Heavner

, Smallbone

, Barker

et al. Yeast 5—an expanded reconstruction of the Saccharomyces cerevisiae metabolic network. BMC Syst Biol, 2012; 6,1:55.

Optimizing Cofactor Specificity of Oxidoreductase Enzymes for the Generation of Microbial Production Strains—OptSwap

Abstract

Introduction

Materials and Methods

Modeling and Computational Tools

Model Reduction and Selection of Reaction Set for Knockouts

Selection of Reaction Set for Cofactor Specificity Swaps

MILP Formulation

Results and Discussion

Conclusions

Footnotes

Acknowledgments

Author Disclosure Statement

References