Screening metallo-beta-lactamase inhibitors effectively is of great significance and necessity. However, a good start can greatly help the experimental group narrow down the selection of inhibitors to reduce workload as well as expenses and help the experimental group to find effective inhibitor molecules in a limited amount of time. So, we built a multistep virtual screening model based on molecular docking to discover the therapeutically effective components, which were used to predict potential inhibitors. Firstly, the structure of NDM-23 is constructed by homology modeling with a template of NDM-1 which has the best match and a high crystal structure resolution. After being examined for conformation and energy, the protein is used for molecular docking. Next, four pre-treated databases are docked with NDM-23 through Surflex-Dock method. Then, according to the docking scoring functions and analysis results, we obtain a potential hit candidate list for further testing by the experimental group. Furthermore, molecular docking methods were used to predict the possible ligand-receptor binding model of the selected effective inhibitors in the experimental group.
Part 1 Introduction
As the development of Computational Chemistry, Molecular biology, Medicinal Chemistry and computer technology, CADMD (Computer Aided Drug Molecular Design) has become one of the main methods to design drugs. CADMD allows the computational strategies to permeate all aspects of drug discovery today, such as the virtual screening (VS) techniques for hit identification and methods for lead optimization. Compared with traditional experimental high-throughput screening (HTS), VS is a more direct and rational drug discovery approach and has the advantage of low cost and effective screening. The basic methods of VS are shown in fig.1.
Methods of CADMD
The ligand structure-based drug design(SBDD) methods, such as quantitative structure activity relationship (QSAR), are based on the principle of similarity to find ligands that potentially act in the same way as previous inhibitors. This method has the characteristics of high speed but large error, especially for those molecules that have a lot of rotatable bonds and are flexible. Because our molecular databases have many natural compound molecules with high molecular weight and many rotatable bonds, we select the molecular docking which based on matching criteria as the method of virtual screening. In addition to predicting the affinity and docking scores that are very important in virtual screening, the molecular docking approach can also be used to model the interaction between a small molecule and a protein at the atomic level, which allow us to characterize the behavior of small molecules in the binding site of target proteins as well as to clarify fundamental biochemical processes. The docking process involves two basic steps: prediction of the ligand conformation as well as its position and orientation within these sites (usually referred to as pose) and assessment of the binding affinity. These two steps are related to sampling methods and scoring functions, respectively.
In molecule docking, both the structure of ligands and a protein is indispensable. Since the structure of NDM-23 is unknown, we resorted to homology modelling to obtain ingredients for molecular docking. Homologous modeling, also known as comparative modeling of protein, is a method of converting amino acid sequence into 3D structures based on a known 3D structure of a related homologous protein. This approach relies on the similarity between the target protein and the template protein. Generally, a model established with more than 50% sequence similarity is accurate and sufficient for drug discovery. And when the sequence similarity is less than 30%, the homologous modeling is usually not accurate enough.
Based on the above ideas, the general steps of our model are as follows.
Part 2 Homology construction
Homologous modeling has two fundamental principles. The first principle is that the spatial structures of proteins are determined by its amino acid sequence. So, the primary structure itself is enough to obtain its secondary, tertiary structure, theoretically. The second basic theory is that in evolution, the conservation of tertiary structures is much higher than that of primary structures. Thus, proteins that diverged appreciably in sequence but have detectable similarity will also share common structural characteristics. By aligning and comparing with the homologous protein, we can get the structure of the target protein.
In our project, we used Swiss-model and Modeller for homologous modeling, both of which are free for academic users. For the modeller, we use it on UCSF Chimera because UCSF Chimera provides a graphical interface for modeler.
STEP 1 Blast search for templates
The protein to be modeled is the target, and a related known structure used for modeling is a template. In this step, the goal is to find the template that matches the target best. We use swiss-model and Chimera's Blast Protein tool to search the Protein Data Bank (PDB) for sequences similar to the target. The hits are listed from best to worst.
In the table, higher blast score and identity, lower E value indicate better results. According to the scores, 4ey2A and 4exsA have the highest ratings. Next, by comparing the resolutions of A and B, we chose A as the template for homologous mode bonds.
STEP 2 Verifying the Alignment
Comparative modeling requires a template structure and a target-template sequence alignment. The sequence alignment is important; it controls which residues in the template are used to model which residues in the target, and any inaccuracies in the alignment will result in the application of incorrect constraints during 3D modeling. Regardless of how the sequence alignment was obtained, it should be examined and adjusted as needed before initiating the more computationally intensive 3D modeling calculations. In the previous step, we found a high degree of similarity between the template and the target, so we chose the single-template modeling approach. Next, we view the Blast alignment to see if it looks reasonable.
Black outline box encloses a residue that is in the sequence but not in the associated structure.
According to the above results, there is only one amino acid difference between the template and the target. Also, the mutation did not occur in the loop region.
STEP 3 Run Modeller
Launch the calculation with default settings on the platform of Modeller. Five comparative models will be generated. Keep the ligands and the water molecules. The Modeller run may take several minutes. When the five models have been generated, they will be opened in Chimera and their evaluation scores shown in a Model List dialog. According to the user's guide, GA341 is model score derived from statistical potentials. Normally, a value > 0.7 generally indicates a reliable model, >95% probability of having the correct fold. zDOPE is normalized Discrete Optimized Protein Energy (DOPE), an atomic distance-dependent statistical score, negative values indicate better models. Taken these two factors into consideration, we choose model 001 as our result for this part. Below is Ramachandran Plot of our result. It is clear that the structure built is reasonable.
STEP 4 Model Evaluation
After getting the structure of our target protein, it is necessary to check whether its conformation and energy are reasonable. In this step, we use Procheck and Errat to evaluate the model. Procheck mainly detects whether the Angle between residues is reasonable, and the generated Ramachandran plot serves as the evaluation basis. In general, a model is likely to be of high quality when the residual number in the core is greater than 90%. ERRAT is used to analyze the statistics of non-bonded interactions between different atom types. And When the ERROT SCORE is greater than 85%, it means the model has a good quality. Both of these tests can be performed on the SAVES server of UCLA-DOE. PROSA can be used to evaluate the energy of the model. In the result diagram of PROSA, the Z-score of all experimentally determined protein chains in the PDB database will be displayed. If the Z-score of the target protein is within this range, it indicates that the structure of the target protein is acceptable.
The three tests above shows the reliability of our model in terms of conformation and energy. In addition, the only amino acid difference between our target and the template is the mutation on the alpha helix, which does not involve the change of the ring region, so our model is considered to be reasonable and reliable..
Part 3 Molecule docking
The basic principles of molecular docking are mainly derived from lock-and-key model and induced fit model. In the lock and key model, the spatial shapes of the protein active sites and ligands are complementary, successfully explaining enzyme specificit. But this theory has a big limitation, because it believes that the enzyme and protein are rigid in the interaction. So then, koshland came up with the idea of induced-fit, which is that both the enzyme and the substrate are flexible, and during the docking of the molecules, the ligand and the protein adjust their conformation to get the best fit. Molecular docking technology focuses on the simulation and identification of molecular interactions, including hydrophobic interactions, hydrogen bonds, van der Waals forces, electrostatic interactions, and solvation effects. And the score function was used to optimize the binding position and conformation of ligands and calculate affinity energy and binding free energy to assist virtual screening. In molecular docking, the ligand is often considered flexible while the protein keeps rigid due to the balance between calculation accuracy and calculation speed. By docking the molecules in the library with the receptor protein, we can find potential drug molecules that bind well to the protein according to the affinity score function.
3.1 Preparation of protein target for virtual screening
Under most circumstances, the crystal structure solved by experiments are not entirely accurate, for the possibility of containing non-protein structure, such as crystal water and eutectic molecules. Also, the hydrogens of the structure are normally left out. Therefore, before the actual docking, the preparation of the protein is vital. Firstly, we gather the chemical structure formula in two methods. One is download directly from the website available, another is to call the SKETCH module to construct ligands. In order to perform a batch docking, then we built multiple molecules in the same folder which can be recognized by Sybyl. After getting these basic works done, we moved on to the optimization of structures.
We call the Biopolymer module to prepare the receptor, details are as follow:
- Delete unwanted amino acids and waters.
- Fix the main chain and side chains of the protein.
- Add charges to the terminals of the main chain.
- Add hydrogens to the whole protein.
- Set protonation types for residues.
- Set atom type to AMBER7FF99
- Add Gasteiger-Huckel charge
- Set termination to Gradient 0.001
- Max iterations is set to 1000
3.2 Validation of docking method
To maintain the reliability of the docking protocol, we first extract compound from the crystal structure and perform the docking to ensure that the docking algorithm can reproduce the observed binding mode as present in the parent molecule. To be specific, validation is carried out in order to ensure that ligands bind to the binding pocket in the correct conformation. The structure of NDM-23 beta‐lactamase was retrieved through homology modeling. Ligand from the NDM-23 crystal structure was extracted and docked back into the active site of NDM-23.
First, the inhibitors were extracted and redocked in its parent molecule. Redocked inhibitor was found to interact with the same amino acids of the active site as was in the crystal structure. The score of similarity score of these two conformations was 0.902, indicating that the protocol set for screening the compound library is accurate.
3.3 Ligand molecule preparation
As for ligands, we call Ligand Preparation module to do following adjustments:
- Add Gasteiger-Huckel charge
- Set force field to Tripos
- Set termination to Gradient 0.001
- Max iterations is set to 1000
After minimized the energy of ligands and protein, we obtained an optimized version of them. In order to be more accurate, we then tried to get the global optimal with algorithm of simulated annealing.
Parameters are set as follow:
- Set cycles as 20
- Heat molecules at 1000K for 1000fs
- Anneal molecules at 300K for 1500fs
- None minimization before annealing
3.4 Molecular Docking（Model construction）
Active pocket selection
According to NCBI, we assign the residues of receptor binding cavity as follow: Leu65, Met67, Trp93, His120, His122, Gln123, Asp124, Glu152, His189, Cys208, Lys211, Leu218, Gly219, Asn220, His250. Water molecules in the binding domain were retained, because water-bridge and two zinc ions served as the important electron transfer medium in the process of breaking beta-lactam rings. We chose to retain H2O 598.
Docking and result selection
Surflex-Dock is particularly successful at eliminating false positive results and can, therefore, be used to narrow down the screening pool significantly, while still retaining a large number of active compounds. Surflex-Dock uses an empirically derived scoring function that is based on the binding affinities of protein-ligand complexes and on their X-ray structures. The Surflex-Dock scoring function is a weighted sum of non-linear functions involving van der Waals surface distances between the appropriate pairs of exposed protein and ligand atoms. We used the Total_score to evaluate the docking results.
3.5 Docking results and conclusions
We selected top 300 molecules from nature products and top 1000 chemistry compounds as our predicted results. The experiment group later on tested 1,300 molecules and confirmed 3 natural products and 5 chemistry compounds are positive inhibitors.
|NAME||FORMULA||DOCKING RESULTS TOTAL_SCORE||EXPERIMENT RESULTS INHIBITON RATE|
|Ginkgoic Acid C15:1||11.4099||99.87%|
3.6 Possible receptor-ligand binding model
After getting the inhibitors screened by the experimental group, we redocked these molecules to NDM-23 to explore their possible receptor-ligang binding models. In this step, we comprehensively considered C-score and whether the ligands bind into the recipient cavity. The docking results of natural compounds and some small molecules are shown below. In these results，it was widely observed that electron dense parts, such as hydroxyl, ketone, ester, sulfonamide and carboxyl, interact with zinc ions and amino acid residues. In addition, hydrophobic effect also plays an important role in maintaining the stability between inhibitor molecules and receptor proteins.
Figure A shows the binding model of Tannic Acid to NDM23. Tannins are polyphenols. In the docking results, Tannic Acid extensively interacts with amino acid residues near the active pocket of NDM-23 through its benzene ring and phenolic groups on the benzene ring, including His122, Asp124, Asn220, GLY222, Glu152, Lys211, Asp223, etc. All these residues play a key role in the activation of the enzyme. Also, the oxygen atoms of the hydroxyl group on the oxygen six-membered ring can directly interact with the Zn2+ ions in the enzyme activity center. The complex formation of hydrogen bonded networks contributes greatly to the high affinity of Tannins, which may explain why it has the lowest IC50 in experiments. In addition, NDM-23 has a large and open active pocket, which may facilitate the binding of natural compounds with larger carbon skeletons, such as Tannins.
Figure B shows the binding model between (-)-Epigallocatechin gallate and NDM-23. In combined mode, similar to Tannins, hydroxyl groups on the benzene ring are also involved in the interaction of Zn2+ ions in the active pocket. This shows the potential activity of polyphenols as NDM-23 inhibitors, which can inhibit the enzyme activity by interacting with Zn2+. Hydrogen bonding of Asp223 and Glu152 also assisted the binding of (-)-Epigallocatechin gallate as well, reflecting indicating the important roles of Asp223 and Glu152 play in similar inhibitors. However, unlike Tannins, (-)-Epigallocatechin gallate has can play an important hydrophobic interactions with Leu65 and Trp93, Met67 in the loop 3 of enzymatic activity center. These hydrophobic interactions are vital and instructive for the reason that they can be used as targets in lead compound optimization.
Figure C shows the binding model of Ginkgolic Acid C15:1 with NDM-23. Shape matching and carboxyl group are important reasons why Ginkgo Acids are effective. As shown in the diagram, the carboxyl group directly interact with two positively charged Zn2+ ions and Lys211 through charge attraction to the active center. Meanwhile, Asn220 is also influential to the localization of carboxyl group.
Figure D shows the potential binding pattern of Corilagin within the active site. Unlike other natural compounds, Corilagin has showed a low docking score (3.5939), which was consistent with the mechanism of non-competitive inhibition measured by the experimental group. This reflects reveals that it may not act in the active pocket of the enzyme, but at other sites of the enzyme. Since metallo-beta lactamases exhibit have significant sequence variability outside the active site, this may be used to explain why Corilagin is not an effective does not exhibit inhibitor totion of other kinds of metallo-beta lactamases.
Figure A shows the docking result between Adaparin and NDM-23. Adaparin's adamantine alkyl matched perfectly with the enzyme's binding pocket. The proper use of adamantine alkyl is expected to improve the specific interactions of drugs. Also, the carboxyl group had a relatively close relationship to Zn2+ in the docking results. The significant interaction between carboxyl and Zn ions are missed in the figure. It might be ascribed to the flexibility of zincs of the NDM family in the active pocket, but the flexibility of Zn2+ is not considered in molecular docking process.
Figure B shows the possible binding mode of Silver Sulfadiazine. In the docking mode, one sulfadiazine oxygen atom interacts with Zn2+ while the other oxygen atoms are stabilized by Asn220. The hydrophobic action of His122 is also crucial in the docking.
Figure C is a schematic diagram of the docking mode of Dicloxacillin Sodium, showing the typical binding mode of beta-lactam antibiotics to the NDM family. Figure D shows the binding mode of Sertaconazole Nitrate with NDM-23. The imidazole group indicates the complexation of Zn2+. In addition, hydrophobic effects of Leu65 and Trp93 also assist in molecular binding.
Analysis of the docking results showed that most of these inhibitors have binding effects on zinc ions in the active pocket center. In addition, several inhibitors interact with amino acids near the active pocket through their electron dense parts such as hydroxyl groups and carboxyl groups. These inhibitors can be characterized by polyphenols. Glu152 Asn220, Gly222 Asp223, Lys211, Asp124 are often participate in the hydrogen bonding interaction with polyphenols. For example, Tannic Acid, showing a high inhibitory effect on metal lactamase, can be used as the direction of future skeleton research.
In addition to phenolic groups, sulfonamides, imidazoles, carboxyl groups, amantadyl groups, ester groups, etc., have also shows significant the potential as active groups of drug molecules. It is noteworthy that the hydrophobic interactions caused by amino acid residues such as Leu65, Trp93 and His122 may play a significanpredominantt role in the binding of small molecules to receptors. To conclude in this step, we discussed the possible binding mode and exact inhibition mechanism of the screened inhibitors, which laid a foundation for further structure-based optimization in the future.
 Schwede T, Kopp J, Guex N, et al. SWISS-MODEL: An automated protein homology-modeling server. Nucl. Acids Res. 2003, 31(32):3381~3385.
 Sali A and Blundell T L. Comparative protein modelling by satisfaction of spatial restraints. Mol. Biol. 1993, 234(14):779~815.
 Covalent Inhibition of New Delhi Metallo‐β‐Lactamase‐1 (NDM‐1) by Cefaclor[J]. Pei W. Thomas,Michael Cammarata,Jennifer S. Brodbelt,Walter Fast. ChemBioChem. 2014 (17).
 Laskowski R. MacArthur M. Moss D. et al. PROCHECK: a program to check the stereochemical quality of protein structures. Appl. Cryst.1993. 26:283~291.
 Chris C. Todd O. Verification of protein structures: Patterns of nonbonded atomic interactions. Prorein Science, 1993. 2:1511~1519.
 Markus W, Manfred J. ProSA-web: interactive web service for the reco- gnition of errors in three-dimensional structures of proteins. Nucleic Acids Research 2007. 35:407~410.
 Bernice Wright,Kimberly A. Watson,Liam J. McGuffin,Julie A. Lovegrove,Jonathan M. Gibbins. GRID and docking analyses reveal a molecular basis for flavonoid inhibition of Src family kinase activity[J]. The Journal of Nutritional Biochemistry,2015,26(11).
 Homology modeling and virtual screening of inhibitors against TEM‐ and SHV‐type‐resistant mutants: A multilayer filtering approach[J] . Mohammad H. Baig,Vishal M. Balaramnavar,Gulshan Wadhwa,Asad U. Khan. Biotechnology and Applied Biochemistry. 2015 (5).
 Zhang Yanhui. Research on homology modelling and molecular docking of P-glycoprotein[D]. Tianjin University.2010.
 Lu yulin. Virtual screening on inhibitors of New Delhi metallo-beta-lactamase 1[D]. Huazhong Agricultural University,2016.
 Rotondo Caitlyn M,Wright Gerard D. Inhibitors of metallo-β-lactamases.[J]. Current opinion in microbiology,2017,39.
 Linciano Pasquale,Cendron Laura,Gianquinto Eleonora,Spyrakis Francesca,Tondi Donatella. Ten Years with New Delhi Metallo-β-lactamase-1 (NDM-1): From Structural Insights to Inhibitor Design.[J]. ACS infectious diseases,2019,5(1).
 Ju Lin-Cheng,Cheng Zishuo,Fast Walter,Bonomo Robert A,Crowder Michael W. The Continuing Challenge of Metallo-β-Lactamase Inhibition: Mechanism Matters.[J]. Trends in pharmacological sciences,2018.