Team:UESTC-China/Model
Overview
Our project is to construct a three-enzyme system to degrade environmental pollutant 1,2,3-trichloropropane (1,2,3-TCP) into glycerol. The system consists of three enzymes DhaA31, HheC and EchA, which catalyze five-step conversions [1]. The efficiency of such multi-enzyme processes can be enhanced by optimization of the stoichiometry of the biocatalysts. Therefore we built a mathematical model, and the enzyme stoichiometry of each step in the system was optimized. Mathematical modeling together with one-pot multi-enzyme experiments were performed to provide detailed insights into system dynamics, enabled us to select a suitable engineered enzyme, and optimal ratios of enzymes. Our model was validated by using sets of experimental data collected from the recent published article. The results showed that all these data after being adjusted could be fitted well using our model, indicating that the validated model can be used to optimize our biodegradation system. We used the model to guide the design of our system, and to improve the efficiency of our system as well.
Figure 1. The biodegradation pathway of 1,2,3-TCP
Model Foundation
∞ Basic Model
To explore the kinetics of a multi-enzyme system, it is difficult to manipulate the process by taking into account of all related factors, such as the reversibility of each enzymatic reaction, the inhibition of the various products on the enzyme, and the affinity of the enzyme for different substrates, etc. It is difficult to study at the same time. Thus, the system needs to be simplified.
Scheme (1)
We first look at each reaction as independent kinetics; we suppose that the reaction catalyzed by all three enzymes is consistent with the Michaelis-Menten equation. Under the above simplifying assumptions, the each step of the whole system can be described as the reaction scheme (1) (E is the enzyme, S, is the substrate, P is the product ), where complex represents for the enzyme and substrate complex (k1,k2 means reaction rate constant). The reaction rate versus substrate concentration at a constant concentration of the enzyme was described in equation (1). So the conversion of the prochiral 1,2,3-TCP into 2,3-DCP by DhaA31 was described using equation (2).
The formation of 2,3-DCP and the subsequent conversion into ECH were described using equations (3) and (4). The consecutive conversions of CPD and GDL by HheC and EchA were described using equations (5) and (6), respectively. The final reaction step, i.e. the conversion of GDL into the final product GLY by EchA, was described using Michaelis-Menten equations (7).
∞Improved model
1. The competitive effect of two substrates on the same enzyme
Although we obtained the GDL expression by the relationship between the Michaelis-Menten equation and the chain reaction, the model only considers the case of single enzyme reaction, and the affinity of the enzyme to different substrates is too complicated to analyze. But the reaction process has two enzymes, that were used twice, we hope that through this aspect to improve the model. During this process, we assumed. We assume that there is no inhibitory effect of the substrate during the reaction
Scheme (2)
As the Scheme 2 shows, we can use HheC and EchA as the same reaction process to establish the model, the two substrates together to share the same enzyme can be approximated as competitive inhibition of the kinetic equation. (E is the enzyme, S, I are the two substrates, ES, EI are the intermediate complex, P, Q are the product)
In Eqs. (8), E is Total enzyme concentration, after the reaction began, divided into [ES], [EI], and the free enzyme [Ef].
Eq. (10) is constructed based on the 'steady state' theory and the law of quality, combine with Eqs. (9), we can get the Eqs. (11) – (14), then the changes in the concentrations of substrates (i.e 2,3-DCP), product 1,2,3-TCP, and intermediates (i.e ECH, CPD, GDL), can be calculated by using Eqs. (3) - (7) and Eqs. (11) - (14).
2. Activity differences toward two enantiomers of 2,3-DCP
HheC have the high enantioselectivity (ER≥100), which made the accumulation of S-enantiomer of 2,3-DCP. (Figure 1) Nonselective DhaA converts the prochiral 1,2,3-TCP into both enantiomers of 2,3-DCP in an almost equimolar ratio. Because HheC prefers (R)-2,3-DCP, (S)-2,3-DCP tends to accumulate. The conversion of the prochiral 1,2,3-TCP into either (R)-2,3-DCP or (S)-2,3-DCP by DhaA31 was described using equation (15).
Figure 2. 1,2,3-trichloropropane (1,2,3-TCP) is converted to both (R) and (S) enantiomers of 2,3-dichloropropane-1-ol (2,3-DCP)
The rate constants kcat,TCP, (R)-DCP and kcat,TCP,(S)-DCP were determined from the overall kcat for the appropriate DhaA variant with 1,2,3-TCP and the ratio of (R)- and (S)-2,3-DCP produced from the prochiral 1,2,3-TCP by DhaA31 and DhaA. The formation of (R)-2,3-DCP and (S)-2,3-DCP and the subsequent conversion of both enantiomers into (R, S)-ECH was described using equations (16) and (17). The remaining reaction steps were assumed to be non-selective. The consecutive conversions of (R,S)-ECH by HheC and EchA was described using equations (17).
3. Inhibition of products
Because the product has an inhibitory effect on the reaction, The GDL concentration progress curve (Eq. 20) was fitted to equation (19), where P is the concentration of the reaction product and Ki is a product inhibition constant[2].
Model fitting with experimental data
Recently, several article related to the studies on the same biodegradation pathway of 1,2,3-TCP have been published. We validated our models by adopting the data sets from these articles[2].
The conversion of glycidol (GDL) into GLY was described by a Michaelis–Menten equation, with one inhibition constants defining the inhibitory effects of GLY (product inhibition constant Ki=1.00 mM).
The model simulation curves based on first setting data did not match well with the experimental data. especially for the time required for reaching the equilibrium. It was not surprising that experimental data and simulation curves based on the first set of parameters did not fit well because kinetic parameters of an individual enzyme may not be the same as its apparent catalytic activity in the mixture of multiple enzymes. The first set of parameters need be adjusted for the best fitting. Taking into account other studies on this pathway[1], and by referring to our own experimental kinetic parameters, we adjust the kcat of the DhaA31 and HheC. Among a number of parameters tested, it was found out that adjusting three kinetic parameters.（i.e. increasing kDhA31cat,TCP,(R)-DCP from 0.58 to 1.33 s-1, increasing kDhA31cat,TCP,(S)-DCP from 0.47 to 1.31 s-1 increasing kDhA31cat,(S)-DCP from 0.08 to 0.33 s-1 ,Simulationcurves based on the second set of data exhibited the best fitting (Fig 3). Then we used the fitting model to optimize our system.
Figure 3. (A) simulation results of the synthesis of GLY from 1,2,3-TCP using experimental data collected from Dvorak[2]. Experimental conditions: 10mM Tris-SO4 buffer (pH 8.5) containing 2 mM 1,2,3-TCP, 2.4mg DhaA31, 2.4mg HheC, 2.4mg EchA at 37℃(B) Under the best fitting, the curves of 1,2,3-TCP and various intermediates, glycerol concentration over time.
Model guiding system optimization
∞ Optimal enzyme loading ratio
The important function of modeling is the prediction of the optimal reaction conditions (e.g. enzyme ratio and catalytic constant of the enzyme). Predictive simulation of the kinetic model was conducted based on the second set of parameters under conditions: in the 10mM Tris-SO4 (pH=8.5) containing 2 mM 1,2,3-TCP, DhaA31, HheC, EchA at 37℃.
The enzymes have similar molecular weights (34.1, 29.3, and 36.5 kDa, respectively), so mass ratio roughly equals molar ratio. Degradation has relationship with these two factors (i) expression profile, (ii) viability. So our kinetic model was employed for rational selection of suitable plasmid combinations and balancing the stoichiometry of the enzyme DhaA:HheC:EchA. The pathway was optimized toward faster removal of toxic metabolites and higher production of GLY.
Significant efforts have been invested in the past few years in engineering the first enzyme of the 1,2,3-TCP pathway: haloalkane dehalogenase (DhaA)[2]. Constructed mutants(DhaA31) could possibly help overcome the previously mentioned limitations of the 1,2,3-TCP pathway in the last studies[2]. So we chose the DhaA31 to construct plasmid. We then evaluated the effects of enzyme stoichiometry on efficiency. An algorithm was used to minimize the total enzyme load without sacrificing productivity (we define the sum of the three enzymes in contrast to 1.0)[3]. We set the step to 0.01, through the algorithm traversal of all the ratio of the situation. The modeled time courses for individual reactions were similar in each case, but the optimized ratios and total enzyme loadings differed significantly between pathways (Fig 4). We find that the optimal enzyme loading ratio is DhaA31: HheC: EchA=0.45:0.45:0.10.
Experiment at this enzyme ratio was conducted to validate the model prediction. As shown in Fig.4, the simulation curves of TCP, intermediates and the glycerol fit well with experimental data, indicating the validity of the model.
As shown in Fig 5, the simulation curves of 1,2,3-TCP and glycerol was consistent with our prediction of the optimum enzyme ratio, which indicated that DhaA31 and HheC had greater effect on degration of TCP.
Figure 4. The four - dimensional relationship between the product and the relative enzyme ratio. Color showing how biocatalyst concentration ratio affects productivity in the three-enzyme conversion of 1,2,3-trichloropropane into glycerol.
Figure 5. The simulation curves of 1,2,3-TCP and GLY in different enzyme loading ratio
∞ Improvement of project
According to the optimal enzyme loading ratio, we found that the ratio of DhaA31:HheC≈1:1. So we design the same promoter to drive the enzymes. 2A Peptides were used to control the expression of the polygene on the same vector is consistent with its protein activity and expression efficiency.
There was epoxide hydrolase in the plant. Besides, we considered the introduction of more foreign genes had exogenous inhibition of gene expression. Based on the optimum enzyme ratio previously obtained, loading quality of EchA had little effect on glycerol yield. We constructed a new plasmid which only had the first two enzymes to improve the productivity of glycerol.
Because the rate of degradation is related to expression profile and viability. However, overexpression of exogenous genes is harmful to the plant. So we paid attention to the viability of enzyme. According to our team in iGEM2013, the activity of HheC-W249P is wild type ten times. So we constructed our plasmid by this mutant.
∞ Modelling guided enzyme mutation
In our work, we found that the viability of HheC-W249P can also be improved, so we focused on obtaining a reasonable and ideal structure by molecular docking. Through the screening of PDB database, we found X-ray crystal structure of HheC with serial number of 1pwz. (Fig 6)
Figure 6. HheC structure diagram (PDB serial number 1pwz)
First of all we had to do some simple processing, the use of Warren Lyford DeLano developed by the open source visualization of the PyMOL software 1pwz selected as the monomer WT, and the system of water molecules removed. We used alanine scanning technology, the active pocket around the 4A within the amino acid ALA replacement. The predicted structure of the mutant was obtained by SWISSMODEL homologous modeling. The WT and mutant models were subjected to Autodock removal of water, plus polar hydrogen, plus Kollman charge treatment. When setted the docking simulation parameters, it was reasonable to consider whether the substrate and the enzyme S132 and Y145 could form hydrogen bonds. Therefore, W139A, F186A, Y187A may be a potential mutation hot spots. (table 2)
Table 2. The docking results of the distance analysis results
| NAME | s132 | y145 |
|---|---|---|
| WT | 9.2 | 7.6 |
| F82A | 8.9 | 7.5 |
| E85A | 8.9 | 7.4 |
| F86A | 8.9 | 7.5 |
| W139A | 3.0 | 3.0 |
| P184A | 9.1 | 7.5 |
| Y185A | 8.8 | 7.4 |
| F186A | 3.0 | 2.9 |
| Y187A | 2.9 | 2.9 |
| P188A | 9.9 | 7.4 |
| W192A | 8.9 | 7.4 |
In the figure 7, red is the active center of the enzyme, and the yellow circle shows the tryptophan at position 139 in WT. As the volume of the 139W side chain is too large, the substrate is prevented from approaching the active site. On the right, the substrate can be docked smoothly when tryptophan is mutated to a smaller side chain.
Figure 7. (A) WT docking result (B) W139A docking result
Next, we will target these mutations and saturation mutations, screening for 2,3-DCP has a highly effective activity of mutants. We constructed plasmids through different mutants, and our kinetic model was employed to improve the 1,2,3-TCP degradation pathway using engineered enzyme variants and balanced enzyme ratios. And model will show excellent correspondence with experimental data. An approximate estimation of expression profiles obtained from sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) analysis of cell free extract (CFE) will show the relative ratio of enzymes which will be constructed by different mutants and promoter. Combining viability of mutants with optimal ratios of enzymes for the maximal production of glycerol, we predict improvements required for further pathway optimization.
References
- Dvorak, P., et al., Immobilized synthetic pathway for biodegradation of toxic recalcitrant pollutant 1,2,3-trichloropropane. Environ Sci Technol, 2014. 48(12): p. 6859-66.
- Dvorak, P., et al., Maximizing the efficiency of multienzyme process by stoichiometry optimization. Chembiochem, 2014. 15(13): p. 1891-5.
- Kurumbang, N.P., et al., Computer-assisted engineering of the synthetic pathway for biodegradation of a toxic persistent pollutant. ACS Synth Biol, 2014. 3(3): p. 172-81.
