Overview

Data from LCA and studies from literature on Vanillin production in E.coli from PET suggests metabolic engineering approaches to increasing yields of high-value products (HVPs) from upcycling PET is essential to make these processes economically viable. We explored potential pathways for improving the downstream production yield of products which could be converted into HVPs such as BHB. Our plan was to depolymerize PHB, which in itself is a high-value product, into D-BHB; therefore we explored pathways to increase PHB production in P.putida and E.coli. According to Zhang et al. (2014), we identified gene and protein targets of which we hypothesized would increase the production of Acetyl-CoA and NADPH, which downstream would result in an increase in production of PHB via the PHB biosynthesis pathway. Our main aim was to increase NADPH and Acetyl-CoA yield. We utilized FBA, a constraint base metabolic modeling method, to determine the distribution of flux across reactions and pathways in the organisms, which helps predict the organisms' growth rate and production of certain metabolites. The objective reactions were G6PD and PDH. Glucose-6-phosphate dehydrogenase (G6PD) is responsible for generating most of NADPH in the pentose phosphate (PP) pathway. (Lim et a., 2002) Similarly, research showed that pyruvate dehydrogenase (PDH) expressed in a yeast cytosol enabled production of Acetyl-CoA, playing a key role in the production process. (Zhang et al., 2020) We utilized the preexisting models iJN1463 for Pseudomonas putida KT2440 and iEC1356_Bl21DE3 for E. coli BL21(DE3). Using the model of Pseudomonas putida KT2440, we determined that zwf increased NADPH yield. The team then modelled PHGHD production by serA using ODEs through MATLAB, which is one of the enzymes involved in the FBA. A similar process was done with zwf and through enzyme kinetics, In order to obtain upper bound of reaction for FBA we modelled PHGHD production by serA using ODEs through MATLAB. Similar process was done with zwf and through enzyme kinetics, the team was able to obtain the upper bound of reaction that could be used for the FBA.

Flux Balance Analysis

Flux balance analysis (FBA) is a constraint-based metabolic modeling method that allows for the simulation of the metabolism of an organism. It effectively determines the distribution of fluxes across reactions and pathways in different organisms. Furthermore, FBA can help predict an organism’s growth rate and the production of a certain metabolite after engineering the organism. (Rowe, 2018)

FBA can be simulated using various methods, such as COBRA Toolbox and COBRApy. In the modeling, COBRApy was used. This Python package can be used for the creation of a model with metabolic reactions, metabolites, and genes. Moreover, the package can be used on an existing model. (Ebrahim, 2013) An existing model from the BiGG Model database, iJN1463, was used to optimize the process of converting TPA to PHB. This model comprises 2,153 metabolites, 2,927 reactions, and 1,462 genes in Pseudomonas putida KT2440.

FBA allows the optimization process to be done with an objective reaction. The objective function refers to the goal that the model aims to optimize. Our optimization process had two objectives: G6PD (Glucose 6- phosphate dehydrogenase) and PDH (Pyruvate dehydrogenase). These two reactions produce Acetyl-CoA and NADPH, which are key metabolites needed in the conversion process. In this case, the objectives refer to the reactions in which the model aims to send the most flux. Flux is a representation of the rate of flow of metabolites through its metabolic network, and it is expressed in mmol/gDW/h.

The weighing of the two objectives can be adjusted. Weight is expressed with a decimal number between 0 and 1, inclusive. For our purposes, 3 different weightings were used: 0.5 (G6PD) and 0.5 (PDH), 0.3 (G6PD) and 0.7 (PDH), and 0.7 (G6PD) and 0.3 (PDH).

Objective

The main aim of this process is to increase the NADPH and Acetyl-CoA yield by editing the model. FBA allows knocking out a single or double gene, increasing certain genes, and more. The following list was obtained from Wet Lab, and FBA was conducted to determine which change effectively increased the yield of NADPH and Acetyl-CoA.

1. Increase serA, serB, and serC
2. Increase serA
3. Increase zwf
4. Increase zwf and serA
5. Increase sdaA
6. Increase fbaA
7. Knock out pgi
8. Coexpression of sdaA, serA, and pgk

Limitations

Throughout the process, we encountered a major problem. Increasing a certain gene required increasing the bounds of the reaction associated with that gene. However, the upper bounds that could be increased to were unknown. Therefore, the only step we could take at this stage of modeling was to set a theoretical lower and upper bound and determine whether the modification was effective or not. However, the FBA was still run with pre-determined bounds to gauge the theoretical impact of adjusting the expression of genes.

Results

In summary, an increase in serA and zwf was successful in significantly increasing the NADPH yield, zwf resulting in a higher yield of the two. However, the acetyl-CoA yield stayed at about the same level. The other changes made to the model had minimal effect on NADPH and Acetyl-CoA yield. The team moved onto determining maximum concentration for each targets: zwf and serA. Below is the workflow of modelling PHGDH production by serA.

Modelling PHGDH production by serA

Objective

The protein 3-phosphoglycerate dehydrogenase, also known as PHGDH, is coded for by the gene serA and is an enzyme catalyzing the oxidation of 3-phospho-D-glycerate to 3-phosphonooxypyruvate, which starts the L-serine biosynthesis pathway, the pathway synthesizing the non-essential amino acid L-serine (UniProt, n.d.). This means the production of PHGDH, and hence the transcription and translation rate of the gene serA, is a major factor in the production of this amino acid and, thus, that changing the rate at which PHGDH is produced can greatly change the rate of L-serine production.

This study intends to find the maximum protein production rate of 3-phosphoglycerate dehydrogenase (PHGDH), which is coded for by the serA coding sequence under the T7 promoter. We aimed to estimate the production rate by calculating ordinary differential equations (ODEs) through the MATLAB SIMBiology software. These values obtained form our ODEs would be used to determine the upper reaction bounds of the reaction (oxidation of 3-phospho-D-glycerate to 3-phosphonooxypyruvate) and be incirporated into our FBA, and thus model how the system would respond to our metablic enginnering, with respects to Acetyl-CoA and NADPH production. The basis of this model is a genetic construct that uses the T7 promoter (BBa_I719005) to produce PHGDH from the serA gene (BBa_K2086001) for enhanced production of polyhydroxybutyrate (PHB).

ODEs & Matlab

Both transcription into mRNA and translation into protein of T7 RNA polymerase (RNAP) and serA must be modeled to predict the expression of the target protein, D-3-phosphoglycerate dehydrogenase. Four ODEs were created as a slightly simplified version of the full model. These equations assume steady-state lacI repression, meaning that IPTG fully derepresses the lac operon, ignore resource competition, and treat the reactions that represent binding between lacI and IPTG and between T7 RNAP and the T7 promoter as negligible.

Lac promoter-regulated production of T7 RNA polymerase. The production of T7 RNA polymerase must be measured because the binding of IPTG to lacI initiates lac operon-regulated transcription of T7 RNA polymerase, ultimately allowing serA to be transcribed under the T7 promoter. The rate of change of T7 RNA polymerase mRNA can be described as:

$\frac{d\lbrack mRNA_{T7 RNAP} \rbrack}{dt}$ = $k_{plac}$ * $C_{N}$ * $\frac{\lbrack IPTG \rbrack ²}{K_{IPTG} ² + \lbrack IPTG \rbrack ²}$ - $d_{T7mRNA}$ * $\lbrack mRNA_{T7 RNAP} \rbrack$

where $k_{plac}$ is the lac promoter transcription rate, $K_{IPTG}$ is the half-saturation constant of IPTG, $d_{T7mRNA}$ is the degradation rate of T7 RNA polymerase mRNA, and $C_{N}$ is the plasmid copy number. Assume [IPTG] is abundent and not limiting.

Translation of the T7 RNA polymerase can then be modeled by creating a first-order ODE for the change in T7 RNA polymerase over time in the form:

$\frac{d\lbrack T7 RNAP \rbrack}{dt}$ = $k_{RBS}$ * $\lbrack mRNA_{T7 RNAP} \rbrack$ - $d_{T7prot}$ * $\lbrack T7 RNAP \rbrack$

T7 promoter-regulated production of D-3-phosphoglycerate dehydrogenase. Expression of T7 RNA polymerase allows it to bind to the T7 promoter, activating transcription of the serA gene under the control of the T7 promoter and thereby translation of the target protein, PHGDH. Transcription of the serA gene can be modeled by the ODE in which RNA polymerase acts as the activator, yielding:

$\frac{d\lbrack mRNA_{serA} \rbrack}{dt}$ = $k_{pT7}$ * $C_{N}$ * $\frac{\lbrack T7 RNAP \rbrack ¹}{K_{T7 RNAP} + \lbrack T7 RNAP \rbrack ¹}$ - $d_{SmRNA}$ * $\lbrack mRNA_{serA} \rbrack$

In this case, the Hill coefficient is 1 because T7 RNA polymerase exhibits non-cooperative DNA binding, meaning it does not require multiple binding events or additional transcription factors to initiate transcription (Kuzmine et al., 2003). Thus, translation of serA mRNA into PHGDH can be modeled by the ODE:

$\frac{d\lbrack serA \rbrack}{dt}$ = $k_{RBS}$ * $\lbrack mRNA_{serA} \rbrack$ - $d_{Sprot}$ * $\lbrack serA \rbrack$

Figure 1: SIMBiology diagram of serA expression using the T7 expression system including binding reactions in green, transcription in yellow, translation in blue, and degradation in black.

Table 25: The reactions involved in the production of PHGDH based on the SIMBiology model.

Figure 2: Full ODEs with additional reactions as defined in Figure 2 generated with MATLAB SIMBiology. Variables in blue represent species, and variables in green represent parameters.

Table 26: Parameters for the ODEs from literature in µM or molecules/min-1 unless otherwise specified.

Initial values of lacI, IPTG, and the T7 promoter also affect the behavior of the model. The initial conditions for all other species are set as 0 because both lac- and T7-regulated transcription are initially repressed by lacI and lack of T7 RNA polymerase expression. For further investigation, the values of lacI and T7 promoter are held at 40 molecules based on the plasmid copy number. IPTG concentration can be varied but is initially set as 1 mM for the basic expression models of mRNA and protein.

Results:
SerA mRNA levels display a sharp peak at 1.74 µM 10 minutes after induction with 1000 µM IPTG, then a rapid decline as mRNA is translated into PHGDH protein and degraded.

Figure 3: serA mRNA expression in µM over time in minutes after induction with 1000 µM IPTG.

After induction with 1000 µM of IPTG, PHGDH protein expression rapidly increases from a baseline expression level of zero to saturation at approximately 0.84 µM by 35-40 minutes after induction, indicating steady-state protein expression at that point. A slow decline thereafter indicates PHGDH degradation. The protein levels in Figure 4 rise significantly later at around 20 minutes compared to the peak in Figure 3, consistent with the expected translation lag after transcription.

Figure 4: PHGDH protein expression in µM over time (minutes) after induction with 1000 µM IPTG

Varying IPTG concentrations affects serA mRNA and protein expression. As expected, at zero IPTG, no expression occurs due to tight repression of the T7 system without the inducer. All induced conditions display a sharp peak in serA mRNA at 10-15 minutes as transcription occurs; higher IPTG concentration increases the peak height, with approximate maximal expression reached after 1000 µM. For PHGDH protein expression, all curves show a rapid increase that plateaus, and maximum expression levels show minimal differences between induction with 1000 to 4000 µM IPTG. The maximum expression levels of the system are approximately 1.8 µM of mRNA and 0.84 µM of PHGDH. An implication of the variations shown in Figure 5 is that 500 µM IPTG is sufficient for near-maximal induction for protein expression, with only minimal increases in expression past 1000 µM IPTG.

Figure 5: PHGDH protein and serA mRNA expression in µM over time after induction with 0 to 4000 µM of IPTG in 500 µM increments.

Results

The key takeaways from the modelling included that amount of IPTG inducer for our constructs which utilize T7 promoter is no more than 0.5 mM IPTG (figure 5), and the maximum concentration of protein produced is 0.84 µM (Figures 4 & 5). A similar work was done with the second target, zwf, and the resulting maximum concentration was the same as that of serA. All the parameters were the same, except the degradation rate, where the difference was negligible. Hence, the maximum concentration rate of zwf is 0.84 µM. Through enzyme kinetics, the team then aimed to find upper values that could be plugged into the FBA to determine exact yields. Zwf was chosen as the final enzyme to increase yield as it was more effective compared to serA.

Enzyme Kinetics 

Overview

The research on enzyme kinetics first started to find FBA constraints, and the team decided that the Michaelis-Menten equation would be an appropriate model for illustrating the relationship between the substrate concentration and the rate of enzyme-catalyzed reaction. We conducted an analysis of kinetics on zwf, one of the enzymes involved in the FBA. 

Michaelis-Menten Model (Berg, Stryer, Tymoczko, & Gatto, 2015)

The Michaelis-Menten equation is a basic equation used to show the kinetics of different enzymes. It's written as:

V = $\frac{V_{\max}\lbrack S\rbrack}{K_{M}+ \lbrack S\rbrack}$

Where: 

To apply the Michaelis-Menten equation to our model, specifically zwf, we first conducted initial research on zwf's kinetics to find known values of $K_M$ or $V_{\max}$ in P. putida KT2440. A relevant paper found was Large-scale kinetic metabolic models of Pseudomonas putida KT2440 for consistent design of metabolic engineering strategies by Milenko Tokic and his group of researchers (Tokic, 2020), where they used an online database, BRENDA, for information regarding enzyme kinetics. This database cites a different literature, Quantifying NAD(P)H production in the upper Entner-Doudoroff pathway from Pseudomonas putida KT2440. (Olavarria, 2015). The paper includes the following table with kinetic parameters of zwf.

A theoretical $K_M$ value of 14 μm is obtained from the literature, with the conditions including pH 8.0, 30℃, recombinant enzyme, and NADP+. Another kinetic parameter included in the table is $k_{cat}$, the enzyme turnover number that represents the number of substrate molecules that are turned over or form a product in 1 second. For zwf in P.Putida, $k_{cat}$ is 102 $s^{-1}$. 

The total zwf concentration was previously calculated to be 0.84 μm, therefore

Now that the two unknown variables for Michaelis-Menten's equation are identified, this model can be plotted in MATLAB. (Figure 1.)

Another application of the Michaelis-Menten Curve is the Lineweaver-Burk Plot, where $\frac{1}{\lbrack S\rbrack}$ is plotted against $\frac{1}{V}$. The slope of the resulting linear plot is equivalent to $\frac{K_{M}}{V_{\max}}$. The derivation of the equation starts by taking the reciprocal of the Michaelis-Menten Model:

A typical equation for a linear function is y = mx+b:

Figure 2 is a representation of a typical Lineweaver-Burk plot, showing how the slope, y-intercept, and x-intercept relate to the kinetic parameters of an enzyme. Similarly, the Lineweaver-Burk plot can be plotted for zwf as seen in Figure 3. The Lineweaver-Burk plot is another method of representing enzyme kinetics. 

The Michaelis-Menten equation was utilized to find the $V_{max}$ value using known kinetic parameters of zwf from literature. The equation was plotted in the form of a Michaelis-Menten curve and a Lineweaver-Burk plot for visualization of kinetic parameters of zwf. The resulting Vmax value was used as an upper bound to constrain the fba model.

Conversion

FBA flux units are typically expressed in mmol/gDW/h. The maximum serA concentration was calculated to be 0.84 µM. It’s assumed that the max zwf concentration is also 0.84 µM, since the two are under the same genetic circuit with the same promoter, and the only difference would be the degradation rate, which had negligible difference.

\begin{align} 1 \mu M &= 1 mmol/1000 L \\ &= 0.001 mmol/L\\ 0.84 \mu M &= 0.00084 mmol/L \end{align} \begin{align} V_{max} &= k_{cat} \times [E] \\ &= 102 s^{-1} \times 0.00084 mmol/L \\ &= 0.08568 mmol/L/s \\ &= 0.08568 mmol/L/s \times 3600 \\ &= 308.448 mmol/L/h \end{align}

The dry cell weight concentration OD of 1 corresponds to 0.49 g/L of CDW (Verhoef, 2010). An OD of 0.8 is used in our design, so dry cell weight concentration is equal to 0.49 0.8 = 0.392 g/L. 

$308.448 mmol/L/h \div 0.392 gDW/h = 786.86 mmol/gDW/h$

This value corresponds to the upper value of the zwf that will be run in the FBA. For simplicity, the lower value was simply the negative of the upper value.

FBA-ZWF Results

Conclusion

In summary, the FBA has proven that zwf increase is effective in increasing NADPH yield, which is crucial for conversion from TPA to PHB. Enzyme kinetic modelling was used to determine upper bound values for zwf. With the exact max enzyme concentration and reaction rate, the team was able to determine exact yields from FBA. However, further research is needed to improve the Acetyl-CoA yield. Following these results, we confirmed our decision to utilize zwf in the experimentations, since it was also most effective in increasing NADPH and following PHB yield, a high value product produced from PET plastic.

References