Analysis of metabolic network disruption in engineered microbial hosts due to enzyme promiscuity.

Increasing understanding of metabolic and regulatory networks underlying microbial physiology has enabled creation of progressively more complex synthetic biological systems for biochemical, biomedical, agricultural, and environmental applications. However, despite best efforts, confounding phenotypes still emerge from unforeseen interplay between biological parts, and the design of robust and modular biological systems remains elusive. Such interactions are difficult to predict when designing synthetic systems and may manifest during experimental testing as inefficiencies that need to be overcome. Transforming organisms such as Escherichia coli into microbial factories is achieved via several engineering strategies, used individually or in combination, with the goal of maximizing the production of chosen target compounds. One technique relies on suppressing or overexpressing selected genes; another involves introducing heterologous enzymes into a microbial host. These modifications steer mass flux towards the set of desired metabolites but may create unexpected interactions. In this work, we develop a computational method, termed Metabolic Disruption Workflow (MDFlow), for discovering interactions and network disruptions arising from enzyme promiscuity - the ability of enzymes to act on a wide range of molecules that are structurally similar to their native substrates. We apply MDFlow to two experimentally verified cases where strains with essential genes knocked out are rescued by interactions resulting from overexpression of one or more other genes. We demonstrate how enzyme promiscuity may aid cells in adapting to disruptions of essential metabolic functions. We then apply MDFlow to predict and evaluate a number of putative promiscuous reactions that can interfere with two heterologous pathways designed for 3-hydroxypropionic acid (3-HP) production. Using MDFlow, we can identify putative enzyme promiscuity and the subsequent formation of unintended and undesirable byproducts that are not only disruptive to the host metabolism but also to the intended end-objective of high biosynthetic productivity and yield. As we demonstrate, MDFlow provides an innovative workflow to systematically identify incompatibilities between the native metabolism of the host and its engineered modifications due to enzyme promiscuity.

Introduction

Integrating heterologous synthesis pathways within microbial hosts has been instrumental in the biomanufacturing of industrial products such as biofuels, polymers, pharmaceuticals, therapeutics, flavors and chemical commodities (Lee et al., 2008; Madison and Huisman, 1999; Nakamura and Whited, 2003; Trantas et al., 2015; George et al., 2015). One strategy to improve yield is to use well-established metabolic engineering techniques such as gene deletion, promoter engineering, media optimization, etc. (Lee et al., 2009). Another strategy is to directly engineer enzymatic properties such as activity, selectivity, inhibition-resistance, and solubility (Yoshikuni et al., 2008). Using one or more of these strategies has proven effective in the development of strains with desired target yields, productivity, and titers.

Often, such metabolic engineering strategies yield unexpected enzyme-compound interactions. Some interactions can be beneficial for the survival of the host. For instance, Patrick et al. documented 41 rescue instances where the lethality of an essential protein deletion was suppressed by overexpression of one noncognate E. coli gene, attributing some of them to catalytic promiscuity and substrate ambiguity (Patrick et al., 2007). In other cases, beneficial interactions can come at a cost. The overexpression or knockout of enzymes can result in interactions that are disruptive for growth and maintenance by siphoning off key metabolic intermediates like pyruvate, acetyl-CoA, and NADH. For example, while seeking to suppress lethality of inactivating the pyridoxal-5-phosphate (PLP) cofactor synthesis pathway, Kim et al. experimentally identified a four-step serendipitous pathway in E. coli that restored the strain’s ability to grow on glucose at the expense of consuming essential intermediates in the native serine biosynthetic pathway (Kim et al., 2010), producing toxic byproducts (Kim and Copley, 2012).

The presence of high concentrations of heterologous enzymes and metabolites within microbial cells causes unexpected promiscuous interactions with host enzymes and metabolites. For example, the short-chain dehydrogenase YMR226C used to produce 3-hydroxypropionic acid (3-HP) in yeast is associated with 15 known substrates (Fujisawa et al., 2003; Jessen et al., 2015). In E. coli strains featuring the malonyl-CoA pathway for 3-HP synthesis, significant quantities of lactate and acetate are produced, even after lactate dehydrogenase (ldh) and pta were knocked out (Rathnasingh et al., 2012). Another instance of promiscuous activity can be observed for pathways intended for butanol production: the promiscuity of the bifunctional butyryl-CoA dehydrogenase (AdhE2) enzyme with substrate acetyl-CoA often results in concomitant synthesis of ethanol with butanol (Inui et al., 2008; Atsumi et al., 2008a; Nielsen et al., 2009). In yet another example, Liao and colleagues leveraged promiscuity of ketoacid decarboxylase (KVD) and alcohol dehydrogenase (ADH) enzymes to synthesize a spectrum of alcohols from branched-chain amino acid metabolic intermediates (Atsumi et al., 2008b). Yet, a caveat of this promiscuous activity is that no single alcohol can be made alone. That is, pyruvate itself is a ketoacid, which can be converted to ethanol by the same promiscuous activity of KVD and ADH. Thus, isobutanol synthesis is also coupled to ethanol synthesis due to enzyme promiscuity (Trinh et al., 2011).

While ubiquitous (D’Ari and Casadesús, 1998; Nobeli et al., 2009; Khersonsky et al., 2006; Tawfik and D. S, 2010) and often observed, effects of enzyme promiscuity on the host metabolic network are often ignored during design and only identified during experimental studies. Predicting such interactions early in the design cycle could yield improved design outcomes and reduce experimental efforts. The prediction of enzymatic products due to substrate promiscuity has mainly relied on hand-curated rules that capture well-known enzymatic transformations. For example, a list of 50 reaction rules, each associated with one or more reaction, was previously defined to explore novel synthesis pathways (Cho et al., 2010; Campodonico et al., 2014). Another set of rules was applied repetitively to generate novel synthesis (Li et al., 2004) or degradation pathways (Finley et al., 2009). Further use of such rules allowed the compilation of over 130,000 hypothetical enzymatic reactions that connect two or more KEGG metabolites (Hadadi et al., 2016), and the compilation of predicted metabolic products into databases such as MINEs (Jeffryes et al., 2015). MyCompoundID (Li et al., 2013) utilizes a similar paradigm and generates products by the repeated application of addition or subtraction of common functional groups. BioTransformer (Djoumbou-Feunang et al., 2019) predicts derivatives by utilizing five separate prediction modules in concert with machine learning and a rule-based knowledge base. PROXIMAL (Yousofshahi et al., 2015) utilizes enzyme-specific reactant–product transformation patterns from the KEGG database (Moriya et al., 2010) as a lookup table to predict products for query molecules. The PROXIMAL algorithm was utilized to create organism-specific Extended Metabolic Models (EMMs) that extend reference metabolic models catalogued in databases to include putative products due to promiscuous native enzymatic activities on native metabolites (Amin et al., 2019; Hassanpour et al., 2020). Despite advances in predicting promiscuous products, however, these efforts have not been put forward in a systematic way to analyze metabolic network disruption in engineered microbial hosts due to enzyme promiscuity.

We develop in this paper a computational method, Metabolic Disruption Workflow (MDFlow), to analyze the disruptive impact of enzyme promiscuity on engineered microbial hosts in a systematic manner. We define metabolic disruption as off-target changes in host metabolism that arise because of enzyme-substrate interactions upon gene or pathway overexpression, where such interactions neither exist in the wild-type chassis organism nor are they expected due to presence of recombinant enzyme(s). Accordingly, this definition encompasses all enzyme promiscuity that arises because of adding heterologous enzymes and their chemical products to a host microbe. Therefore, MDFlow is designed to consider two different disruption scenarios. In “Scenario 1”, promiscuous activity is predicted in the context of overexpressed enzymes, whether heterologous or native, acting promiscuously on native host metabolites. Meanwhile, in “Scenario 2,” predictions are made by assuming that native enzymes exhibit promiscuous interactions with synthesis pathway metabolites introduced with engineering changes. Of course, in a biological system, both scenarios would occur simultaneously to some degree, and even higher-order interactions would be possible (e.g., subsequent use of promiscuous reaction products as substrates for additional transformations). Using PROXIMAL (Yousofshahi et al., 2015) and flux analysis (Orth et al., 2010; Segrè et al., 2002), MDFlow models and evaluates such promiscuous scenarios. We demonstrate the use of MDFlow to evaluate how engineered microbial hosts are impacted by enzyme promiscuity under three engineering strategies. First, MDFlow is used to explain how essential gene deletion can be suppressed by overexpression of another native gene. Next, MDFlow is used to identify multi-step interactions that may compensate for essential gene knockouts. Third, MDFlow is used to evaluate the potential disruption when a heterologous pathway is added to a microbial host. The first two cases represent Scenario 1 and are evaluated against experimentally verified data. The last case represents a simultaneous application of both Scenarios 1 and 2 and demonstrates that the choice of synthesis pathway can impact metabolic disruption scenarios.

This work is novel as it is the first to systematically investigate effects of heterologous and native enzyme promiscuity on host metabolism and its consequences on biocatalysts. Our method serves as the first computational tool that can assist metabolic engineers in (1) identifying sources of unexpected byproducts, (2) assessing the consequences of metabolic engineering on the host, and (3) quantifying pathway-host incompatibility using metabolic network disruption. Outcomes from this work will aid in future studies to design robust systems with more predictable behaviors and improved desired product yield.

Methods

MDFlow (Fig. 1) is an integrated method that combines several techniques to analyze cellular metabolic disruption that may occur due to enzyme promiscuity. To predict such putative interactions, MDFlow relies on PROXIMAL (Yousofshahi et al., 2015) to predict byproducts resulting from promiscuous activities under Scenarios 1 and 2. Stoichiometrically balanced reactions are then derived based on substrate, product, and the reaction associated with the promiscuous enzymatic transformation pattern. To assess the disruption impact in a systematic fashion, the host metabolic network model is incrementally modified and evaluated after each change – first to set a baseline, then augmented with the engineering strategy of interest. Quantitative flux analysis (Flux Balance Analysis (FBA) (Orth et al., 2010) and/or Minimization of Metabolic Adjustment (MOMA) (Segrè et al., 2002)) was then used to evaluate the impact of such changes on biomass growth rates or product yield. We provide a detailed overview of each of these steps.

Fig. 1

An overview of the four-step process used by MDFlow to identify and evaluate byproducts formed due to enzymes promiscuity for Scenarios 1 and 2. The original host metabolic model is progressively augmented with engineered modifications and predicted interactions. The updated models are evaluated using FBA and/or MOMA at different stages.

Step 1 – establish the baseline growth rate and/or target yield using FBA

To model the E. coli metabolic network, we built upon the iML1515 model published by Monk et al. (Monk et al., 2017). When evaluating gene knockout modifications, we used a derivative of that model, iML1428, that offers improved accuracy of lethality predictions in such experiments. iML1428 achieves this improvement by removing isozymes that are minimally expressed in glucose M9 media, thus preventing them from incorrectly compensating for the removal of essential genes. We used the conditions suggested by the authors of iML1515 for constraint-based modeling – the lower bounds of all exchange reactions were set to zero, except for glucose, oxygen, and all inorganic ions: the lower bound for each of those reactions was set to –10, –20, and –1000 mmol gDW–1 h–1, respectively. With the constraints configured for aerobic growth on glucose, we evaluated the baseline growth rate, and target metabolite yield (if applicable), of the host using FBA. Then, we set the lower bound on biomass growth to be equal to 10 % of the baseline growth rate to prioritize minimum growth of the strain required to ensure its long-term survival. This lower bound may be increased or decreased as appropriate for the application.

Step 2 – evaluate direct impact of the engineering change using flux analysis

We implemented intended changes for metabolic engineering in the context of addition or removal of reactions and metabolites in the network. To construct a synthesis pathway for transforming a metabolite within a host into a target compound, we added each synthesis step to the stochiometric matrix (S-matrix) of the metabolic model as a new reaction, along with a demand reaction for the target metabolite. To model a gene knockout, we enacted the effect of deleting the gene from the model by setting flux to zero for all inactive reactions, identified based on their individual gene-protein-reaction (GPR) rules (Schellenberger et al., 2011). The new “engineered” model was then evaluated using FBA to demonstrate the impact of introduced changes. At this stage in the workflow, it is important to verify that the model’s prediction reflects expectations or experimental data.

Step 3 – predict promiscuous interactions using PROXIMAL

The modified model from the previous step was used to predict promiscuous byproducts for Scenarios 1 and 2. Each added reaction was assumed reversible unless indicated otherwise. For each scenario, PROXIMAL first created a lookup table of all known biotransformation operators based on catalogued reactions within the model. These operators encode molecular transformation patterns associated with an enzyme and the reactions it catalyzes. Given a query molecule, PROXIMAL first identified operators in the lookup tables that can act on the molecule. Then, the query molecule was mapped to possible byproducts. An operator acts on a query molecule if its reaction center and its nearest neighbor atom(s) exactly match those of the native substrate, as encoded in the lookup tables. Each potential byproduct was reported in the form of a mol file, which was then used to identify if the potential byproduct is a known metabolite. The mol file was matched to either a metabolite in the model, a KEGG ID, or a PubChem ID using InChIKeys (Heller et al., 2015) generated by an open-source chemical toolbox RDKit (RDKit).

For Scenario 1, biotransformation operators were derived from the overexpressed enzymes along the synthesis pathway and applied to all native metabolites in the model. For Scenario 2, biotransformations were derived from native host enzymes within the model and applied to pathway metabolites in the engineered pathway. Applying PROXIMAL operators resulted in a list of byproducts for each scenario. For each predicted byproduct, we developed a new balanced enzymatic reaction based on the catalyzing enzyme’s reaction pattern. A reaction template with suitable cofactors was obtained from reaction(s) associated with each enzymatic biotransformation. A reaction is balanced when the number of atoms for each element on the reactant side matches those on the product side. Reactions were verified to be balanced using ChemPy (Dahlgren, 2018). If a reaction was not balanced, it was discarded and not considered for further analysis as it violated the assumptions of FBA.

Step 4 – evaluate network disruption using flux analysis

To predict the effect of promiscuity, the engineered model network was augmented with balanced promiscuous reactions. The metabolites in these reactions were either mapped to existing metabolites in the model or added to the model along with the corresponding exchange reactions allowing unlimited export (but not import) of a metabolite from the host. Each predicted reaction was first added to the engineered model separately. The flux range of each promiscuous reaction was calculated by minimizing (“min”) and maximizing (“max”) its flux. Since the new reactions did not have any particular direction associated with them, and the most disruptive one was not yet known at this stage, flux ranges calculations were performed twice: once assuming a forward (“fwd”) flux direction and another assuming the reverse (“rev”) flux. As a result, four flux values were estimated for each reaction: vfwd_min, vfwd_max, vrev_min, and vrev_max. The maximum reverse flux (vrev_max) and minimum forward flux (vfwd_min) may be trivial solutions; however, they could be non-zero for reactions that were required for the model to maintain the minimum growth rate. The added reactions were either new to the model or formed biotransformation routes that were catalyzed differently than those already in the model. This updated network represented the “disrupted” model.

To model promiscuous activity, we assumed a non-zero flux through the added reaction (vfwd_added, vrev_added) based on a pre-defined percentage p, referred to as a coupling percentage, of the calculated maximum flux. As in case of all reactions when performing FBA, the absolute maximum for any reaction was set at ±1,000 ​mmol/h/gDW. The constraints imposed on the reaction can be described by the following inequalities – the first one pertains to the forward flux through the reaction, while the second describes the reverse flux:–1000 ​≤ ​

Once constraints were set, the metabolic network disruption caused by the developed reactions could be evaluated. The disrupted yield value was calculated and compared with the value of the yield of the undisrupted engineered model. The direction of each reaction was then fixed to the one that caused maximum disruption. For the purposes of subsequent analysis, reactions were considered irreversible and their individual flux was constrained by only one of the two inequalities. Disruption in yield was then evaluated under various coupling percentages, with a random subset of developed reactions added to the undisrupted model. The results were then placed on a scatter plot relating the extent of metabolic disruption to the intensity of promiscuous activity for visualization purposes.

Implementation

The workflow was written in Python, targeting versions 3.7 and up, and is available on GitHub (https://github.com/HassounLab/MDFlow). Instructions for running the workflow are provided on the project’s GitHub page and in the README file accompanying the source code. Promiscuous reactions predicted by PROXIMAL (Yousofshahi et al., 2015) (https://hassounlab.cs.tufts.edu/proximal/) are provided in the repository. The FBA implementation and model handling logic were provided by the COBRApy package (Ebrahim et al., 2013), and all visualizations were built using matplotlib (Hunter, 2007), seaborn (Waskom, 2020), and RDKit (RDKit). All experiments were conducted on an Intel® Core™ i7-2600 machine with 8 ​GB of RAM, though the workflow does not require the use of any specific hardware platform.

Results

MDFlow explains how suppressors can rescue growth of essential gene deletions

We first demonstrated that MDFlow can be used to predict how enzymatic promiscuity suppresses the effect of gene deletions. We validated MDFlow against experimental data suggested by Patrick et al. (Patrick et al., 2007). The authors set out to identify and categorize multifunctional genes that enabled cells’ adaptability to genetic lesions. Within the Keio collection (Baba et al., 2006), the authors identified 107 single-gene knockout strains that were unable to grow on the M9-glucose medium. Patrick et al. found 21 strains that could be rescued via overexpression of one non-cognate gene, for a total of 41 unique suppression examples. Comparing structural superimpositions of deleted proteins and their suppressors, Patrick et al. observed significant structural homology between a handful of enzyme pairs and attributed several examples to substrate ambiguity and catalytic promiscuity (Patrick et al., 2007). The fact that an essential gene deletion was suppressed by overexpression of another native gene and that enzyme/substrate non-specificity played a role in the rescue suggests that the responsible mechanism in many instances could have been a single enzymatic reaction mediated by promiscuous activity of overexpressed enzymes – that is, activity following Scenario 1. In this test case, we applied MDFlow to each of the 41 unique knockout-suppressor pairs and evaluated its effectiveness in explaining the suppression of gene deletion via overexpression of another native gene. The observations made by Patrick et al. were thus used to both guide the workflow and validate it against experimental data.

To establish the baseline for validating rescue, we first determined if the deletion of each of the 20 essential genes responsible for the 41 pairs arrests the growth of the strain. To this end, we performed gene knockouts in iML1428 and measured the biomass growth rate reported by FBA. Out of the 20 lethal knockouts reported by the authors, we were able to confirm lack of growth for 9 of them, accounting for 15 out of the 41 unique knockout-suppressor pairs. Since we cannot computationally validate rescue without first establishing the lethality of the genetic lesion, we only focused on those 15 cases where the deletion led to no biomass growth. For each of the suppressor genes, we used PROXIMAL to generate promiscuous reactions due to Scenario 1 type interactions (e.g., overexpressed gene acting on native metabolites) and then applied FBA to determine their individual effect on the biomass, resulting in the complete recovery of 6 strains (ΔilvE, ΔglnA, ΔilvD, ΔhisH, ΔpabA, and ΔilvA). In the majority of the recovered cases, a promiscuous reaction replicated the deleted essential reaction either exactly (e.g., ΔilvA/tdcB, ΔpabA/pabB, and ΔhisH/hisF) or with different cofactors (e.g., ΔilvE/avtA and ΔglnA/asnB). Both pabA/pabB and hisH/hisF form well-characterized heterodimeric enzymes consisting of a larger subunit (PabB and HisF) and a smaller one (PabA and HisH). In both cases, the larger subunit has shown activity in the absence of the smaller subunit (Ye et al., 1990; Klem and Davisson, 1993). Therefore, overexpressing PabB and HisF could be reasonably expected to compensate for the deletions in the ΔpabA and ΔhisH strains. Meanwhile, in the ΔilvA/tdcB case, the two proteins are known to be isozymes of one another, thus having the same function (Patrick et al., 2007). Therefore, there is biological justification for the rescuer enzyme’s potential ability to exactly replicate the reactions lost by the deletion. From the perspective of our method, both the deleted gene and its suppressor had the same Enzyme Commission (EC) numbers, which made them interchangeable for the purposes of transformation pattern extraction. In the other 3 case sets, however, the outcome appears to be due to enzyme promiscuity – although the first two or three EC groups matched in two instances, the deleted protein and its replacement were not identical enzymes in all three. For ΔglnA/asnB, the promiscuous activity was responsible for creating a reaction edge between L-glutamate and L-glutamine, which appears to be the mechanism of recovery. For ΔilvE/avtA, the promiscuous reaction allowed the production of L-isoleucine, which is essential for survival as double mutants (ΔilvE ΔavtA) are known to require supplementation of both, isoleucine and valine, for growth (Whalen and Berg, 1982). Notably, the deletion of ilvE also caused the loss of an L-valine-producing reaction in the model, however, its essentiality was not confirmed through FBA. For ΔilvD/avtA, no single predicted reaction was responsible for growth recovery, however, including various combinations of two predicted reactions demonstrated rescuing effect. When using MOMA, the lethality and rescuing effects of knockouts and multicopy suppressors were identical as predicted when using FBA. A summary of multicopy suppression results can be found in Table 1 along with example compensating reactions predicted for each case. A complete set of predicted reactions is provided in Supplementary File 1.

Table 1

Confirmed-lethal single-gene deletion strains from Patrick et al. (Patrick et al., 2007) with the corresponding multicopy suppressors (MS) and representative subsets of predicted compensating reactions. Biomass growth rate in the disrupted network computed using FBA is presented as a percentage of the wild type strain growth rate.

Deletion	Deficient Reactions	MS	MS Biomass	Compensating Reactions
ΔilvA (4.3.1.19)	L-threonine → 2-oxobutanoate ​+ ​NH4+	tdcB (4.3.1.17, 4.3.1.19)	100.00%	L-threonine → 2-oxobutanoate ​+ ​NH4+
			100.47%	L-allo-threonine → 2-oxobutanoate ​+ ​NH4+
ΔilvE (2.6.1.42)	(S)-3-methyl-2-oxopentanoate ​+ ​L-glutamate → 2-oxoglutarate ​+ ​L-isoleucine3-methyl-2-oxobutanoate ​+ ​L-glutamate → 2-oxoglutarate ​+ ​L-valine	avtA (2.6.1.66)	100.00%	(S)-3-methyl-2-oxopentanoate ​+ ​L-alanine → L-isoleucine ​+ ​pyruvate
			100.79%	3 3-methyl-2-oxobutanoate ​+ ​2 ​L-alanine → 2 ​L-isoleucine ​+ ​3 pyruvate
ΔglnA (6.3.1.2)	ATP ​+ ​L-glutamate ​+ ​NH4+ → ADP ​+ ​L-glutamine ​+ ​H+ ​+ ​PO43-	asnB (6.3.5.4)	100.94%	AMP ​+ ​L-asparagine ​+ ​L-glutamate ​+ ​P2O74- → L-aspartate ​+ ​ATP ​+ ​L-glutamine ​+ ​H2O
			101.26%	L-glutamate ​+ ​NH4+ → L-glutamine ​+ ​H2O
			99.38%	AMP ​+ ​L-asparagine ​+ ​D-glucose 1-phosphate ​+ ​L-glutamate → ADP-glucose ​+ ​L-aspartate ​+ ​L-glutamine ​+ ​H2O
			(and others)
ΔilvD (4.2.1.9)	(R)-2,3-dihydroxy-3-methylbutanoate → 3-methyl-2-oxobutanoate ​+ ​H2O(R)-2,3-dihydroxy-3-methylpentanoate → (S)-3-methyl-2-oxopentanoate ​+ ​H2O	avtA (2.6.1.66)	100.23%	2 Pyruvate ​+ ​3 ​L-valine → 2 (S)-3-methyl-2-oxopentanoate ​+ ​3 ​L-alanine2 4-methyl-2-oxopentanoate ​+ ​3 ​L-alanine → 2 pyruvate ​+ ​3 ​L-valine
			100.23%	3 3-methyl-2-oxobutanoate ​+ ​2 ​L-alanine → 2 ​L-isoleucine ​+ ​3 pyruvate2 ​L-leucine ​+ ​3 Pyruvate → 3 3-methyl-2-oxobutanoate ​+ ​2 ​L-alanine
			(other combinations may be possible)
ΔhisH (4.3.2.10, 3.5.1.2)	L-glutamine ​+ ​phosphoribulosylformimino-AICAR-phosphate → AICAR ​+ ​erythro-imidazole-glycerol-phosphate ​+ ​L-glutamate ​+ ​H+	hisF (4.3.2.10)	100.02%	L-glutamine ​+ ​phosphoribulosylformimino-AICAR-phosphate → AICAR ​+ ​erythro-imidazole-glycerol-phosphate ​+ ​L-glutamate
			100.02%	L-glutamine ​+ ​phosphoribosylformiminoaicar-phosphate → AICAR ​+ ​erythro-imidazole-glycerol-phosphate ​+ ​L-glutamate
			100.05%	L-glutamine ​+ ​2 phosphoribulosylformimino-AICAR-phosphate → (S)-2-hydroxyglutarate ​+ ​2 AICAR ​+ ​2 erythro-imidazole-glycerol-phosphate
			(and others)
ΔpabA (2.6.1.85)	chorismate ​+ ​L-glutamine → 4-amino-4-deoxychorismate ​+ ​L-glutamate	pabB (2.6.1.85)	100.00%	chorismate ​+ ​L-glutamine → 4-amino-4-deoxychorismate ​+ ​L-glutamate
			100.00%	2 chorismate ​+ ​L-glutamine → 2 4-amino-4-deoxychorismate ​+ ​(S)-2-hydroxyglutarate
			100.00%	chorismate ​+ ​L-glutamine → 4-amino-4-deoxychorismate ​+ ​O-acetyl-L-serine
			100.00%	L-asparagine ​+ ​chorismate → 4-amino-4-deoxychorismate ​+ ​L-aspartate
			100.00%	L-glutamine ​+ ​isochorismate → 4-amino-4-deoxychorismate ​+ ​L-glutamate
			100.00%	L-glutamine ​+ ​prephenate → 4-amino-4-deoxychorismate ​+ ​L-glutamate

MDFlow can identify multi-step bypasses to essential gene functions

Because MDFlow anticipates interactions between native metabolites and enzymes due to multiple genes that are simultaneously overexpressed, MDFlow can be used to identify multi-step interactions where a native pathway contributes metabolites to an unexpected, serendipitous process. We validated MDFlow for predicting promiscuous interactions in the serendipitous pathway identified by Kim et al. (Kim et al., 2010). The authors identified a four-step chain of interactions that compensated for an essential gene knockout. The pdxB gene deletion disrupted the PLP synthesis pathway, resulting in its inability to grow on M9-glucose at 37 ​°C. Simultaneous overexpression of seven other genes, however, rescued the strain. Using genetic complementation experiments, the authors were able to separate those genes into groups and describe one of the rescuing pathways in detail. The serendipitous pathway, shown in Fig. 2, was found to bypass the knocked-out enzyme in the PLP pathway by diverting flux from serine biosynthesis. The interactions comprising this pathway were catalyzed by three enzymes – two of which exhibited either promiscuity (ThrB) or broad specificity (LtaE). The function of the third enzyme (YeaB) was unknown, and one of the interactions appeared to be non-enzymatic. Therefore, at least part of the pathway could have emerged due to promiscuous activity classified as Scenario 1 type interactions by our method.

Fig. 2

Serendipitous pathway discovered by Kim et al. that bypasses the deletion of an intermediate gene pdxB in the native PLP pathway by siphoning off material from the serine biosynthesis pathway. Circles highlight reaction steps that were predicted – or not predicted – by PROXIMAL. Pathway layout for the drawing was adapted from the authors’ original paper.

The new pathway bypasses the lesion by diverting a metabolite (3-phosphohydroxypyruvate, 3-PHP) from the serine production pathway and converting it to a metabolite (L-4-phosphohydroxythreonine, 4-PHT) downstream the PLP synthesis pathway via a series of 4 ad-hoc reactions (Kim et al., 2010). Assuming that the reactions comprising the pathway are due to promiscuity, discovering each step of the pathway entails predicting interactions between the candidate enzymes – that is, all native enzymes except the knocked out PdxB – and all available metabolites in the model – which includes the native metabolites as well as metabolites generated by promiscuous activity in previous steps. Repeating this process up to four times, in principle, would allow the extraction of the entire pathway. Unfortunately, this approach would be computationally intensive as it amounts to an exhaustive enumeration of all possible promiscuous reaction sequences. However, because the steps comprising the serendipitous PLP pathway are already known, verifying if the exhaustive search would eventually find the pathway can be done efficiently. At each step, instead of predicting interactions for all available metabolites, we only considered the small subset that overlaps with the ΔpdxB bypass pathway, significantly reducing the search space. Using MDFlow configured in this way, we were able to rediscover 3 out of the 4 reactions along this novel pathway. The elusive step, the LtaE-catalyzed transformation of glycolaldehyde to L-4-hydroxythreonine (4-HT), was not predicted due to the lack of an appropriate reaction pattern in PROXIMAL’s operator lookup table. This observation highlights an important shortcoming of rule-based methods such as BioTransformer (Djoumbou-Feunang et al., 2019) and PROXIMAL (Yousofshahi et al., 2015): their metabolism predictions are limited to the finite set of biotransformation rules encoded in their knowledgebase. Nonetheless, the results for this pathway and multicopy suppressors from (Patrick et al., 2007) show that our technique can be used to discover experimentally validated unexpected interactions that lead to the survival of the host.

To further illustrate the utility of MDFlow in predicting multi-step pathways, we used MDFlow to discover two-step promiscuous pathways in iML1515 that may compensate for the deletion of pdxB. Such a two-step pathway transforms a metabolite within iML1515 to one of the following metabolites that are downstream from 4-phosphoerythronate: 2-oxo-3-hydroxy-4-phosphobutanoate, 4-hydroxy-L-threonine, pyridoxal 5-phosphate. Intermediate metabolites at each step were predicted by PROXIMAL using biotransformations specific to E. coli as listed in the KEGG database. Then, we retained pathways that only had intermediates with known PubChem IDs. We further filtered these results based on thermodynamic feasibility of the pathways as estimated by eQuilibrator (Flamholz et al., 2011; Noor et al., 2012; Noor et al., 2013; Noor et al., 2014). There were 21 such feasible two-step pathways. There were 13 and 8 pathways that terminated on pyridoxal 5-phosphate and 2-oxo-3-hydroxy-4-phosphobutanoate, respectively, whereas none that terminated on 4-hydroxy-L-threonine. As the iML1428 model does not capture the lethality of knocking out pdxB, we were unable to evaluate the impact of these two-step pathways on growth rates using FBA. A complete set of predicted two-step pathways is provided in Supplementary File 2.

MDFlow predicts metabolic disruption during 3-hydroxypropionic acid (3-HP) biosynthesis adversely affects yield

We used our workflow to analyze metabolic disruption through Scenarios 1 and 2 for two different synthetic 3-HP pathways, an important precursor metabolite to useful derivatives such as acrylic acid, 1,3-propanediol, and malonic acid (Chen and Nielsen, 2013; Della Pina et al., 2011; Kumar et al., 2013). Several groups have reported production of 3-HP in various organisms (Jiang et al., 2009; Huang et al., 2013), including E. coli (Raj et al., 2008; Kim et al., 2014). Using E. coli as a host, we first added a two-step 3-HP synthesis pathway (Rathnasingh et al., 2012) (Fig. 3A) catalyzed by malonyl-CoA reductase (MCR) (Cheng et al., 2016) and 3-hydroxyacid dehydrogenase (YdfG). Then, separately, we added a three-step 3-HP pathway (Wang et al., 2014) (Fig. 3B) comprising of panD, gabT, and ydfG. The 3-HP yield of the baseline unmodified model was estimated to be 0.90 ​mol/mol of glucose, and the addition of the new pathways increased that to 1.62 ​mol/mol of glucose for the first pathway and 1.66 ​mol/mol for the second. These figures correspond to the maximum theoretical yield of 3-HP assuming no disruption. When using MOMA instead of FBA to evaluate the 3-HP yield, we found that the added reactions do not increase the yield of 3-HP over the baseline. We believe this is a direct consequence of the MOMA approach: when the engineering change is the addition of reactions, the problem of minimizing flux redistribution has a trivial solution that allows no flux through the new reactions. As a result, we do not believe MOMA is applicable to this particular case. Thus, we conducted all subsequent analysis using FBA only. After the addition of a synthesis pathway to E. coli, our workflow utilized PROXIMAL to predict derivatives for both Scenarios 1 and 2. For both Scenarios, identities of predicted derivatives were looked up in either iML1515, KEGG, or PubChem. As not all predicted derivatives were identifiable through this method, we considered only derivatives that are documented in at least one of the three databases. Example reactions for each Scenarios 1 and 2 are shown in Fig. 3, panels C and D, respectively. Other predicted reactions are given in Supplementary File 3.

Fig. 3

Overview of the two heterologous 3-HP pathways integrated into the E. coli model and the method used to construct putative promiscuous reactions for each scenario. (A) Pathway 1, which converts malonyl-CoA into 3-HP via two reactions catalyzed by malonyl-CoA reductase (MCR) and YdfG. (B) Pathway 2 that produces 3-HP from L-aspartate using PanD and GabT in addition to YdfG. Developed reactions examples of Scenario 1 (C) and Scenario 2 (D). Both panels (C) and (D) are divided in three sections: i. the native reaction catalyzed by the potentially promiscuous enzyme, ii. the RDM pattern showing the rction center in red where the biotransformation occurs, and iii. the developed balanced reaction indicating the reactants, products, and the promiscuous enzyme. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Balanced predicted reactions due to each scenario were first evaluated using FBA to determine the range of fluxes they can potentially sustain in the context of the E. coli metabolic model. The model was then augmented with all predicted reactions, along with the 3-HP synthesis pathway. We then assumed that each of the predicted reactions exhibits a certain minimum activity, measured as a fraction of its maximum possible flux. We refer to this fraction as the coupling percentage of a given reaction.

As not all enzymes and metabolites in E. coli are present in high concentrations, not all predicted promiscuous reactions occur at appreciable levels to disrupt the metabolic network. In addition, the level of promiscuous activity is not constant across the network and is likely to vary from reaction to reaction. To address these issues, we performed a two-stage randomized probabilistic analysis. For each scenario, we assumed that only a portion of the predicted promiscuous reactions act on their target molecules. First, a random mean coupling percentage pmean was sampled from a normal distribution with μ ​= ​1 % and σ ​= ​μ/3. The standard deviation was chosen in such way to place 99.7 % of the samples within a range strictly above 0 %, and in the remaining 0.3 % of the cases, the coupling percentage was resampled repeatedly to obtain a value above 0 %, still. Then, we randomly selected 10 % of the developed reactions to exhibit promiscuous activity. For each of the selected reactions, a coupling percentage was chosen randomly from another normal distribution with μ ​= ​pmean and σ ​= ​μ/3. The selected coupling percentage was then used to set the minimum required flux for the reaction according to one of the inequalities presented earlier in the methods section, while the lower and upper bounds of all other (non-selected) developed reactions were set to zero. Sampling the mean coupling percentage thus allowed us to vary the extent of the overall promiscuous activity, while sampling individual reaction-specific coupling percentages allowed each reaction to behave independently from any other. 10,000 FBA runs were performed, each time selecting a different set of promiscuous reactions and re-sampling the mean and individual reaction coupling percentages. In each run, the FBA objective was to maximize the yield of 3-HP while maintaining biomass growth of at least 10 % of wildtype. The analysis was also repeated assuming 25 % and 50 % of the developed reactions to have promiscuous activity. Additionally, another run was performed, where the total number of active reactions was fixed at 10 and each Scenario was either allocated the entire set (10) or half of it (5). This was done to account for the fact that Scenario 2 reactions were much more plentiful than Scenario 1 reactions. Avoiding overrepresentation of the former reactions enabled qualitative comparison of the two scenarios.

The results in Fig. 4 (A, B, C, and D) for both pathways and scenarios exhibit similar trends. Higher disruption tends to be correlated with more active reactions and with higher coupling percentages. As the mean coupling percentage increases, so does the disruption, and the relationship is almost linear in all but Pathway 1, Scenario 1 (Fig. 3A). The mean coupling percentage was sampled from a normal distribution and the resulting distribution of disruption values followed a very similar bell curve. The same effect occurs with an increase in the number of active reactions – note the mean of the distribution of disrupted yields shifts towards increased disruption with greater reaction activities. In fact, Scenario 2 has between 1.5- and 2-times the number of reactions compared to Scenario 1 (e.g., 34 vs. 17 for pathway 1 and 91 vs. 60 for pathway 2), which leads to a proportional difference in magnitude of disruption observed in the experiments as well, and when the total number of reactions is set to a fixed value as opposed to a percentage, the difference between the two Scenarios almost vanishes (Fig. 4F).

Fig. 4

Comparison of simulation runs for the two pathways and both Scenarios, demonstrating the effect of mean coupling percentage and fraction of active promiscuous reactions (10, 25, 50 %) on the yield disruption of 3-HP. (A, D) Pathways 1 and 2 with only Scenario 1-type promiscuous interactions incorporated. (B, E) Pathways 1 and 2 with Scenario 2 interactions. (C, F) Comparison of Scenario 1 (S1) and 2 (S2) interactions for each pathway given a fixed total number of promiscuous reactions. The scatter plots document the results of each individual experiment, while the distribution plots on the right represent the probability density of obtaining a given disruption under the specified conditions.

Scenario 1 for Pathway 1 was an exception to this trend, however. In that instance, only ~17 % of predicted reactions had any effect on 3-HP yield disruption when measured individually – suggesting that only a small fraction of sampled reactions would have a disruptive impact at any given instant under those conditions. As a result, there is a disproportionately large number of instances where promiscuous reactions led to a low magnitude of disruption, particularly for low coupling percentages (Fig. 4A). The comparison of the two Scenarios for that Pathway further demonstrates the difference in behavior, with Scenario 1 being significantly less disruptive than Scenario 2 on average (Fig. 4C). This outcome highlights the downside of considering Scenario 1 alone: when few enzymes are being evaluated, the results are highly sensitive to the set of transformation patterns that can be derived from them – subject to knowledge limitations and enzymatic behaviors. Under more realistic biological circumstances, a more diverse range of enzymes participates in metabolism, allowing for more representative results. These observations imply that as promiscuous activity is intensified, it directly competes with 3-HP synthesis pathway flux causing significant disruption. However, such high activity of promiscuous reactions is unlikely under physiological settings, so these exercises likely overestimate the actual extent of disruption.

Overall, Pathway 1 was found to experience comparatively less disruption compared to pathway 2 under the same circumstances, possibly due to a lower number of enzymes, fewer intermediates, and/or lower promiscuous activity of MCR relative to PanD and GabT. In light of our analysis and considering the similar yield of Pathways 1 and 2, the design of Pathway 1 may be the preferred design option.

Discussion

In metabolic engineering, a number of strategies are employed to produce a target metabolite of interest, including the introduction of heterologous enzymes and selective overexpression and deletion of certain genes. After such interventions or combinations thereof, it is not uncommon for unexpected and/or undesirable metabolite product profiles to arise from various interactions between the introduced and native machinery in the host. Using MDFlow, it is possible to predict such interactions in the form of single- and multi-step pathways enabled by promiscuous enzymatic activity. The method utilizes PROXIMAL to construct reactions arising from promiscuity and relies on FBA to assess their impact on yield or biomass growth rate in pre- and post-modification hosts. The results for the single-gene deletions that were rescued via single- and multi-step enzymatic pathways, which were experimentally validated in prior studies, provide evidence of the ability of MDFlow in predicting metabolic network disruption. The results for the 3-HP case illustrate how MDFlow can help identify promiscuous interactions early in the design cycle. Further, our sampling-based FBA analysis shows that promiscuity can cause unexpected byproducts and results in yield disruption. Importantly, MDFlow can be used to explain byproducts often observed but not well explained in the literature.

Our Scenario 1 and 2 disruption classification has direct correspondence to the network inference classification proposed by Kim et al. (Kim and Copley, 2012). Interference is classified into three groups: those due to (i) heterologous metabolites in new pathways interfering with native metabolism, (ii) native metabolites interfering with a heterologous pathway, and (iii) heterologous pathway intermediates being diverted by promiscuous activity of native enzymes. MDFlow identifies the same interactions as long as they’re caused by enzyme substrate promiscuity – with groups (i) and (iii) corresponding to Scenario 2 predictions and group (ii) represented by Scenario 1-type interactions.

The computational methodology in MDFlow can be further enhanced. It currently uses PROXIMAL to predict promiscuous byproducts; it is possible to use alternatives for promiscuous product prediction. We selected PROXIMAL because we have established confidence in its capabilities in predicting organism-specific enzymatic transformations. Our earlier study of promiscuity using PROXIMAL on non-engineered E. coli allowed the discovery of 17 putative enzymatic reactions that explained metabolomics measurement (Yousofshahi et al., 2015). Regardless of the tool, however, there always remains the issue of false positives. In our work with PROXIMAL, we discarded byproducts that were not documented prior in PubChem, KEGG, or iML1515/1428. Using machine-learning tools that evaluate the likelihood of compound-enzyme interactions, such as EPP-HMCNF (Visani et al., 2020), might provide further confidence in such predictions. Importantly, better prediction of enzymatic products and their likelihood can improve the workflow’s ability to uncover unexpected interactions and evaluate their impact on the engineered host.

Metabolism simulation aspect of MDFlow also offers opportunities for improvement. In certain applications, MOMA can be an attractive alternative to FBA. While FBA utilizes linear programming to maximize an objective function, typically yield or biomass, MOMA approaches cell modeling from the perspective of minimizing redistribution of metabolic fluxes compared to the wildtype conditions. MOMA can provide improved correlation of predictions with experimental flux data over the steady-state modeling provided by FBA (Segrè et al., 2002). Additionally, the model itself can be tailored to the specific circumstances of an experiment. For example, isozymes that are minimally expressed in glucose M9 media were removed from the iML1515 model to create iML1428. As a result of this reduction in degrees of freedom, the derivative context-specific model tends to offer more accurate lethality predictions than iML1515 in gene knockout experiments (Monk et al., 2017). Such adjustments can be informed by manual investigation of the host’s metabolic network, by leveraging detailed kinetic models (Khodayari and Maranas, 2016), or via techniques such as 13C Metabolic Flux Analysis (Long and Antoniewicz, 2019). Of course, depending on the application, it may be helpful to consider using more complex techniques that can capture translation and regulation aspects of underground metabolism, though we believe – in the three cases we have considered – the overexpressed genes are the dominating factor. Knowing more about the biological sample, such as the concentrations of metabolites and enzymes can shed further light on the amount of disruption.

Another area of improvement is the approach used to determine the direction and maximum flux limits of the predicted reactions. We rely on the preset 1% average mean coupling percentage to estimate the limits of all reactions, which may not be representative of the higher or lower actual mean coupling percentage under the conditions of a given experiment. Future studies may re-evaluate directions and flux limits in the context of each experiment individually. Despite these limitations, the presented results are promising and call for further design exploration of the impact of enzyme promiscuity on engineered microorganisms.

Conclusion

We presented MDFlow, a method for quantitatively evaluating the side effects of engineered modifications on the host metabolic network resulting from enzyme promiscuity. Without mitigation, such side effects may lead to either unexpected behaviors during experiments or failure to take advantage of potentially beneficial interactions. By combining PROXIMAL and flux analysis in a streamlined workflow, MDFlow is capable of both discovering new interactions and evaluating their effects without the need for costly, time-consuming studies of in vivo experiments. MDFlow can be used at all stages of the metabolic engineering process. Prior to any design work, the method can be applied to expose native background promiscuous activity, revealing potentially interfering enzymes and metabolites that may be present but not documented in the model, resulting in building Extended Metabolic Models (EMM), as we demonstrated in prior works (Amin et al., 2019; Hassanpour et al., 2020). The workflow can then actively guide the pathway construction process by providing feedback on engineering decisions as they are being made, thus helping to identify modifications that create more robust strains. And finally, MDFlow may be leveraged to compare a fully designed pathway to other candidate pathways for the same target compound: when implemented, a pathway with less predicted metabolic disruption may have a higher yield. MDFlow is a first systematic automated analysis step towards the evaluation of underground metabolism and its interaction with engineered cellular machinery.

Author statement

Conceptualization: NUN, SH.

Methodology: NUN, SH, VP, SA, VEG, TN.

Software: VP, SA.

Investigation: VP, SA.

Data Curation: VP, SA, VEG, TN.

Funding acquisition: SH, NUN.

Writing - Original Draft: VP, SA.

Writing - Review & Editing: VP, SH, NUN, TN.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.