Model

Overview

Building models that are closely related to the project and can guide each other with experiments is our core, therefore, this part is dedicated to simulate the execution process of circuit by constructing mathematical models under the guidance of experiments to improve the effect of circuit and optimize the experimental design. We hope that the experiments will provide first-hand data for the mathematical model to establish the basic model parameters, and in return, the mathematical model will provide guidance for efficient and optimized experiments and inspiration for experimental design to achieve a better interpretation of the project.

In our project, the main objective is to treat cardiovascular diseases caused by poor diet by reducing the concentration of TMA, therefore, the numerical model will enhance the working performance of the circuit from the perspective of optimizing the circuit design and point mutation to improve the degradation efficiency of TMADH. Meanwhile, to further enhance the therapeutic effect to expand the application scenarios, we try to design a lead compound of inhibitor from the perspective of inhibiting the production of TMA to further enhance the therapeutic effect and expand the application scenarios.

Degradation model of TMA

A model which represents multi-component and time-varying dynamic system is widely used in various biological problems. Among all of the models, differential equation is at the top of the list. To illustrate it, the regulatory network can be represented by an array of ordinary differential equations, where the interaction between a multitude of molecules (such as mRNA or protein) are quantified by the law of mass action. These equations designate the level of each protein or mRNA as a function of other components in the development of the system. These models usually include time-related variables, such as the concentration of protein and mRNA, other parameters, such as protein degradation parameters, are also included.

Therefore, clarifying the interaction and response-related parameters in the regulation network are required to establish gene expression models. Thanks to tons of work done by a large number of cell biologists, we have a deeper insight into the gene regulation networks. After revealing the interaction between genes, the model then can be constructed after finding the corresponding parameters.

In this model, we first need to address whether theophylline can activate the degradation system at a suitable concentration and without leakage in the absence of theophylline administration. Based on the experimental exploration of theophylline induction, we tried to fit the experimental data and find the threshold value for the induction concentration.

Figure 1. The simulation of theophylline induction system.

Based on the results of the simulations, we believe that the threshold of theophylline induction is around 0.883mM, which is much higher than the normal physiological concentration, and that the degradation system is not activated in vivo without theophylline administration. This may explain the fact that the experiments showed similar protein expression at 2 mM and 4 mM theophylline induction. This is because after reaching the threshold of induction concentration, changes in theophylline concentration would not greatly affect the expression of downstream genes.

On the other hand, whether TMA can be effectively degraded under the affect of the degradation system in a shorter period of time and with a lower impact on human body is also an issue we should focus on, so we need to construct models for TMA, DMA, PA and other substances to investigate the degradation efficiency.

Reaction formula

$$P_{T M A D H} \stackrel{\beta_{1}}{\longrightarrow} T M A D H $$

$$P_{D M A D H} \stackrel{\beta_{2}}{\longrightarrow} D M A D H $$

$$P_{P A D H} \stackrel{\beta_{3}}{\longrightarrow} P A D H $$

$$T M A+T M A D H \stackrel{K m_{1}, k_{1}}{\longrightarrow} D M A+PA $$

$$D M A+D M A D H \stackrel{K m_{2}, k_{2}}{\longrightarrow} M M A+P A $$

$$P A+P A D H \stackrel{K m_{3}, k_{3}}{\rightarrow} F A $$

$$T M A D H \stackrel{\alpha_{1}}{\rightarrow} \emptyset $$

$$D M A D H \stackrel{\alpha_{2}}{\rightarrow} \emptyset $$

$$P A D H \stackrel{\alpha_{3}}{\rightarrow} \emptyset $$

ODE

$$\frac{d[T M A D H]}{d t}=\beta_{1}-\alpha_{1} *[T M A D H] $$

$$ \frac{d[D M A D H]}{d t}=\beta_{2}-\alpha_{2} *[D M A D H] $$

$$\frac{d[P A D H]}{d t}=\beta_{3}-\alpha_{3} *[P A D H] $$

$$\frac{d[T M A]}{d t}=-\frac{k_{1} *[T M A] *[T M A D H]}{K m_{1}+[T M A]} $$

$$\frac{d[D M A]}{d t}=-\frac{k_{2} *[D M A] *[D M A D H]}{K m_{2}\left(1+\frac{[T M A]}{K i}\right)+[D M A]}+\frac{k_{1} *[T M A] *[T M A D H]}{K m_{1}+[T M A]} $$

$$\frac{d[M M A]}{d t}=\frac{k_{2} *[D M A] *[D M A D H]}{K m_{2}\left(1+\frac{[T M A]}{K i}\right)+[D M A]} $$

$$ \frac{d[P A]}{d t}=\frac{k_{2} *[D M A] *[D M A D H]}{K m_{2}\left(1+\frac{[T M A]}{K i}\right)+[D M A]}+\frac{k_{1} *[T M A] *[T M A D H]}{K m_{1}+[T M A]}-\frac{k_{3} *[P A] *[P A D H]}{K m_{3}+[P A]} $$

$$\frac{d[F A]}{d t}=\frac{k_{3} *[P A] *[P A D H]}{K m_{3}+[P A]} $$

Figure 2. The degradation model of TMA.

Table 1. List of parameters involved in TMA degradation model.

Based on the above model results, it was found that PA and DMA accumulate during degradation, although the accumulation of DMA is small but due to the toxicity of DMA, this is a result we do not want to see in humans, so we decided to optimize the original genetic circuit.

From the results of simple simulations, it was found that different choices of RBS can have a large impact on the results of the degradation model of TMA, therefore, we need to explore the optimal experimental design under different RBS combinations. In order to ensure the therapeutic effect and applicable safety, the following three evaluation indexes were established.

(1)Degradation efficiency of TMA.

(2)Transformation efficiency of PA.

(3)Accumulation of DMA.

The above indicators were quantified and characterized by the integration of the area of the ordinary differential equation. The comprehensive evaluation model of Topsis [5] is a multi-objective decision analysis method applicable to a limited number of scenarios. We scored the TMA degradation process under different RBS combinations by the comprehensive evaluation model of Topsis based on the quantified results, combined with the characteristics of the data itself and subjective assumptions. Based on the scoring results, we established the optimal RBS circuit design as B0034-B0034-B0030.Then we plotted heat maps to characterize the results under each RBS combination, where the vertical axis represents the RBS corresponding to the upstream of tmd and the horizontal axis represents the RBS corresponding to the upstream of dmd and fdnH, respectively.

Figure 3. The evaluation scores of RBS combination.

In summary, we finally verified the feasibility of the model and improved the original experimental design to establish the optimized circuit.

Point Mutation Optimization Model for TMA Degrading Enzymes

Sensitivity Analysis

Based on the TMA degradation model constructed in Part1, we first performed a sensitivity analysis of the important parameters in it, as shown in Figure 4. According to the images, TMA degradation efficiency is significantly affected by the change of k₁ and Km₁, what's more, the interaction term is highly sensitive. As a result, the enzymatic catalytic efficiency of TMADH will greatly affect the effectiveness of our treatment, and it is therefore necessary to optimise the conformation of TMADH.

Figure 4. The Result of Sensitivity Analysis.

Table 2. The Result of the Sensitivity Coefficient.

Molecular Dynamics Simulation

TMADH is an oxidoreductase that removes a methyl group from TMA, while ETF is an important electron transfer flavoprotein endogenously expressed in E. coli. During the degradation of TMA ,TMADH binds to ETF to form a complex that facilitates electron transfer in the redox reaction. Besides, electron transfer is known to be the decisive step in the degradation of TMA[6].To simulate the combination mode of TMADH and ETF and the effect of point mutation on TMA degradation efficiency, molecular simulations of TMADH and ETF were carried out.

However, in TMADH, the active site consists of 6-S-cysteine FMN and 4Fe-4S centres, both covalently linked to protein residues, and this particular property introduces uncertainty into our molecular modelling work.

Figure 5. The Overall Structure of TMADH.

The protein sequence we chose was consistent with the wild type protein, number 1DJN, and then we selected it as a template on the website Swiss Model[7] to construct a possible protein conformation. To confirm the reliability of the homology modelling results, we uploaded the generated PDB file to Mol Probity[8] and analysed its conformation. The result is that there is only one unreasonable place.

Then it came to the next part of molecular dynamics -preparing the prmtop file and the inpcrd file which contains the force field parameters and the coordinate parameters of the molecule, respectively. We first calculate the relevant parameters by Multiwfn[9] and create the lib files for the non-standard units and proteins in the xleap program of AMBER[10] to obtain the two preparation files. The process is shown below.

Figure 6. Flow Chart of Generating Preparation Files.

Getting the prmtop and inpcrd files, we imported the ff99SB protein force field, the GAFF organic small molecule force field and the TIP3P water model force field into the AMBER tleap program. Afterwards,the input files for solvent minimization, solvent molecular dynamics, system minimization, system molecular dynamics, side chain and solvent molecular dynamics were generated and ran the following process in the force field.

Figure 7. Flow Chart of Molecular Dynamics.

After obtaining the molecular dynamics results, we used capptraj in AMBER to generate the trajectory file without water due to the performance of trajectory analysis and calculate the RMSD of the protein. To gain the representative conformation generated during the simulation, we used cpptraj to get the average structure based on all the different structures observed in the simulation and calculated the RMSD of various conformations based on the average one. This intended to find the best approximation of the average structure,which was the representative.

The next part of our work was optimisationing the conformation we got. It was known that Rosetta, a tool for structure prediction, design and remodelling of proteins and nucleic acids, can be used to model macromolecular structures. One of the commonly used optimisation tools in Rosetta is RELAX, which has a multi-step iterative simulated annealing algorithm, consisting of three steps including weight escalation, amino acid side chain rearrangement and energy minimisation [11]. We used this tool to further optimise the result and correct unreasonable residues, bringing the parameters of each residue within the permissible range.

Figure 8. Ramachandran Plot of the Molecular Dynamics Results (left) and the Rosetta Optimized Results (right).

Molecular Docking

TMADH is a homodimer with each subunit containing an unusual 6-S-cysteine FMN cofactor covalently linked to Cys-30, a ferric oxide-reducing protein type 4Fe-4S centre covalently linked to four cysteines (Cys-345, Cys-348, Cys-351, Cys-364) and an ADP. Moreover, the 4Fe- 4S centre is located within the concave surface of TMADH and is involved in electron transfer in redox reactions.

Figure 9. Several Important Ligands of TMADH.

In order to render TMADH a higher TMA degradation efficiency, we tried to optimize it by point mutation. As is known, residues Tyr-442 and Val-344 are located in the centre of the concave surface of TMADH and the concave region is postulated to form an ETF docking site.Tyr-442 is the amino acid closest to the iron-sulphur centre and is in van der Waals contact with Val-344. Val-344 is sequentially adjacent to Cys-345, which is a ligand for one of the iron atoms. Up to now, a possible electron transfer pathway through Cys-345, Glu-439 and Tyr-442 has now been identified by pathway algorithm [6]. And it is known from experimental data in the literature [12] that mutation of V344C gives a better TMA degradation efficiency. Therefore, we decided to target Tyr-442 and Val-344 for mutation to cysteine.

It is also known through the literature that there is a possible interaction between Tyr442 in TMADH and Argα237 in ETF [13], so both of two were used as docking sites and uploaded to HDOCK[14] for rigid macromolecular docking.

We attempted to dock TMADH with ETF obtained by molecular dynamics, but the docking results were not satisfactory. We speculated that this was due to the conformational change in the molecule during molecular docking and what'more, during rigid docking, the conformation might not match the one simulated by molecular dynamics. Then we used the Rosetta-optimised TMADH to dock with ETF in HDOCK and found that they docked successfully with the highest confidence score (1.0000) and that the distance between Tyr442 (TMADH) and Argα237 (ETF) was within the electron transfer-limited distance (14Å).

Figure 10. TMADH and ETF Docking Results.

Figure 11. Interface Between TMADH and ETF.

We uploaded the docking result to Prodigy[15] to measure its binding free energy and dissociation constant but found that the measured dissociation constants for V344C and wild type were equal due to the low precision of the website measurement method. Surpringly, the results were in good agreement with the literature finding of reduced degradation efficiency for the Y442C mutant and slightly increased degradation efficiency for the V344C mutant. In this way, we believe that this proves the potential of V344C mutant to improve the degradation efficiency of TMA to some extent. Subsequent experiment also verified this potential.

Table 3. Binding Free Energy and Dissociation Constants of Wild Type and Mutant After Docking.

Rigid docking has some limitations and may not be able to simulate the diverse conformation of molecular docking. Based on a review of the literature, we found that ETFs will undergo conformational changes upon binding to TMADH and that ETFs from methylotrophic bacteria are similar in conformation to human ETFs in the activated state[16].Thus we decided to use the human ETF as a template for homology modeling to get the conformation of the ETF that we needed in the activated state.

We used HDOCK to dock the specified residues, and the resulting superior conformation was used as the natural conformation for Rosetta local docking. After obtaining the highly scored conformation, we relaxed it in Rosetta and calculated the binding free energy as well as the distance between the FAD and the iron-sulphur centre. Although the docking results are not ideal, this provides an idea for future exploration.

Virtual screening of CutC inhibitor lead compounds

In the communication with our partnership, they suggested that in order to achieve a better quality therapeutic effect, we could try to reduce the concentration of TMA in two directions: TMA degradation and TMA production. Therefore, in this part We intended to improve the therapeutic efficiency by designing enzyme inhibitors and the project’s application prospect in different scenarios was proposed.

In human intestine, CutC is the main lytic enzyme that converts choline to TMA. We hope to design an easily synthesized and harmless inhibitor to inhibit the activity of CutC by competing with choline for the active site. However, when discussing the biosynthesis with the experimental group, we learned that chemical small molecules are too difficult to biosynthesize, while peptides can be more easily biosynthesized which contributes to our project promotion. As a result, We determined to design a macrocyclic peptide inhibitor to inhibit the activity of CutC enzyme.

First, we were supposed to design the pharmacophore part of the macrocyclic peptide. By reviewing several literature, we found the best competitive inhibitor of chemical small molecule[17], the superior effect of which has been proved through experiments. On the basis of the molecular structure, we used RDkit to screen all the dipeptide structures composed of non-essential amino acids based on the molecular fingerprint, and obtained the following six structurally similar dipeptides for the subsequent validation, as shown in Figure 12.

Figure 12. Molecular similarity between chemical inhibitors and dipeptides.

And then, to verify the competitive inhibition effect of the above dipeptides, we performed 50 molecular docking using SaliVina[18][19] 6 dipeptide molecules and choline with CutC, respectively, and the results of the formation of hydrogen bonds between molecules were counted as shown in the following figure 13 and figure 14. Kruskal-Wallis test was performed on the results, and it was obtained that the docking results of His-Met and Met-His were significantly different from choline, while His-Val and Val-His, His-Leu and Leu-His were similar. Thus we tentatively considered the latter group of dipeptide as possible competitive inhibitors.

Figure 13. The hydrogen bonds interaction of Choline and peptide with CutC.

Figure 14. The comparison of hydrogen bonds between Choline and peptide with CutC.

In order to reflect the inhibition effect of small molecules more visually, we subjected choline and dipeptide to competitive multi-substrate docking with CutC. Using the results of choline docking with CutC as a control, we analyzed it through PLIP[20] which was shown in the figure below with choline indicated in orange. Choline has hydrogen bonding interactions with 784Thr and 785Ser in CutC, Π-Cation interactions with 677Phe, and Salt Bridge interactions with 498ASP in most docking results, and hydrogen bonding interactions with 788Tyr and 976Ser in very few docking results, and Salt Bridge interactions with 776Glu.

Figure 15. The interaction of choline and CutC.

The results of competitive docking of His-Leu, His-Val, Leu-His and Val-His with choline are shown in the table 4 below. We found that His-Leu and His-Val had competitive inhibition effect, for they formed hydrogen bonds or hydrophobic interactions with residues present in some controls and significantly affected the binding of choline to CutC. However, while Leu-His and Val-His were less effective, and the binding of choline to CutC was largely unaffected.

Table 4. Results of multi-substrate molecular docking.

Figure 16. Docking result of His-Leu and choline with CutC.

Figure 17. Docking result of His-Val and choline with CutC.

Therefore, we drew a conclusion that His-Leu and His-Val may have an inhibitory effect on CutC, and their docking results are shown in Figure 16 and 17. It is worth trying to design macrocyclic peptide inhibitors based on them, but unfortunately due to time limitation, we failed to complete the design of the subsequent macrocyclic peptides. However, we thought that such an attempt do offer the future iGEMers of treatment track a feasible idea and contribute to the promotion of our project in more application scenarios.

Reference

[1]Gasteiger E, Hoogland C, Gattiker A, et al. Protein identification and analysis tools on the ExPASy server[J]. The proteomics protocols handbook, 2005: 571-607.

[2]Basran J, Chohan K K, Sutcliffe M J, et al. Differential coupling through Val-344 and Tyr-442 of trimethylamine dehydrogenase in electron transfer reactions with ferricenium ions and electron transferring flavoprotein[J]. Biochemistry, 2000, 39(31): 9188-9200.

[3]Meiberg J B M, Harder W. Dimethylamine dehydrogenase from Hyphomicrobium X: purification and some properties of a new enzyme that oxidizes secondary amines[J]. Microbiology, 1979, 115(1): 49-58.

[4]Ashraf R, Rashid N, Basheer S, et al. Glutathione-dependent formaldehyde dehydrogenase homolog from Bacillus subtilis strain R5 is a propanol-preferring alcohol dehydrogenase[J]. Biochemistry (Moscow), 2017, 82(1): 13-23.

[5]Olson D L. Comparison of weights in TOPSIS models[J]. Mathematical and Computer Modelling, 2004, 40(7-8): 721-727.

[6]Wilson E K, Huang L, Sutcliffe M J, et al. An exposed tyrosine on the surface of trimethylamine dehydrogenase facilitates electron transfer to electron transferring flavoprotein: kinetics of transfer in wild-type and mutant complexes[J]. Biochemistry, 1997, 36(1): 41-48.

[7]Waterhouse A, Bertoni M, Bienert S, et al. SWISS-MODEL: homology modelling of protein structures and complexes[J]. Nucleic acids research, 2018, 46(W1): W296-W303.

[8]Chen V B, Arendall W B, Headd J J, et al. MolProbity: all-atom structure validation for macromolecular crystallography[J]. Acta Crystallographica Section D: Biological Crystallography, 2010, 66(1): 12-21.

[9]Lu T, Chen F. Multiwfn: a multifunctional wavefunction analyzer[J]. Journal of computational chemistry, 2012, 33(5): 580-592.

[10] Case D A, Ben-Shalom I Y, Brozell S R, et al. AMBER 2018; 2018[J]. University of California, San Francisco, 2018.

[11]Tyka M D, Keedy D A, André I, et al. Alternate states of proteins revealed by detailed energy landscape mapping[J]. Journal of molecular biology, 2011, 405(2): 607-618.

[12]Loechel C, Basran A, Basran J, et al. Using trimethylamine dehydrogenase in an enzyme linked amperometric electrode Part 2. Rational design engineering of a 'wired' mutant[J]. Analyst, 2003, 128(7): 889-898.

[13]Burgess S G, Messiha H L, Katona G, et al. Probing the Dynamic Interface between Trimethylamine Dehydrogenase (TMADH) and Electron Transferring Flavoprotein (ETF) in the TMADH− 2ETF Complex: Role of the Arg-α237 (ETF) and Tyr-442 (TMADH) Residue Pair[J]. Biochemistry, 2008, 47(18): 5168-5181.

[14]Remmert M, Biegert A, Hauser A, et al. HHblits: lightning-fast iterative protein sequence searching by HMM-HMM alignment[J]. Nature methods, 2012, 9(2): 173-175.

[15]Honorato R V, Koukos P I, Jiménez-García B, et al. Structural biology in the clouds: the WeNMR-EOSC ecosystem[J]. Frontiers in Molecular Biosciences, 2021: 708.

[16]Leys D, Basran J, Talfournier F, et al. Extensive conformational sampling in a ternary electron transfer complex[J]. Nature Structural & Molecular Biology, 2003, 10(3): 219-225.

[17]Gabr M, Świderek K. Discovery of a Histidine‐Based Scaffold as an Inhibitor of Gut Microbial Choline Trimethylamine‐Lyase[J]. ChemMedChem, 2020, 15(23): 2273-2279.

[18]Trott O, Olson A J. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading[J]. Journal of computational chemistry, 2010, 31(2): 455-461.

[19] O'Boyle N M, Banck M, James C A, et al. Open Babel: An open chemical toolbox[J]. Journal of cheminformatics, 2011, 3(1): 1-14.

[20]Adasme M F, Linnemann K L, Bolz S N, et al. PLIP 2021: expanding the scope of the protein–ligand interaction profiler to DNA and RNA[J]. Nucleic acids research, 2021, 49(W1): W530-W534.

[21]Mulligan V K, Workman S, Sun T, et al. Computationally designed peptide macrocycle inhibitors of New Delhi metallo-β-lactamase 1[J]. Proceedings of the National Academy of Sciences, 2021, 118(12): e2012800118.