Identification of potent natural compounds in targeting Leishmania major CYP51 and GP63 proteins using a high-throughput computationally enhanced screening

Leishmaniasis is a disease caused by protozoan forms called Leishmania which infect animals and humans. The drugs have been in use since half a century due to which there have been mutations in the microbe-facilitating drug resistance. So this provides a reason for searching for effective drugs for the disease. In the current work, an effort has been to find such drugs that act on disease-relevant receptors by similarity indexing method, molecular docking, and dynamics studies. The study focused on the rapid expansion of potential anti-leishmanial compounds that could function as novel natural compound structures for future drug Similarity indexing of existing drugs with natural compounds using Tanimoto clustering resulted in 4 compounds with similarity index of greater than 0.7 (70% similarity). The molecular docking of the resulted compounds was carried out with therapeutic targets CYP51 and GP63 proteins. N-methyltyrosyl-N-methyltyrosyl-leucyl-alanine from Streptomyces griseus showed higher binding affinity in comparison to inhibitor and other selected natural compounds. Simulation studies revealed that the binding configuration of the compound with targets was highly stable all through 10 ns of simulation time with intact hydrogen bonding. The molecular docking and molecular dynamics studies for the selected natural bioactive compound N-methyltyrosyl-N-methyltyrosyl-leucyl-alanine from Streptomyces griseus showed better binding affinity with the selected therapeutics targets and can be further considered for in vitro and in vivo studies which may lead to a possible new drug for the treatment of Leishmaniasis.

Leishmaniasis is given in figure (Fig. 1) [2]. The severity of disease ranges from minor lesions that may heal in few months to years to highly severe lesions leading to death of the person if left untreated. The annual incidence is estimated to be 500,000 cases, with 50% occurring in India. The disease is endemic in the eastern States of Bihar, Jharkhand, Uttar Pradesh, and West Bengal. In these states, it is marginalized communities living in rural poverty who bear the greatest disease burden [3].
Leishmaniasis is transmitted by the bite of infected female Phlebotomine sand-flies. The sand-flies inject the infective stage (i.e., promastigotes) from their proboscis during blood meals. Promastigotes that reach the puncture wound are phagocytized by macrophages and other types of mononuclear phagocytic cells. These organisms inhibit the microbicidal activity of host macrophages, thus passing through the adaptive and innate immune responses of our body. Further, in the promastigote stage engulfed by macrophages with the help of C3b receptors, they resist proteolysis and degradation in the phagosome. These organisms have special surface glycoprotein called GP 63 which converts C3b into iC3b8 which allows phagocytic clearance.
The phagocytic protection includes protection from oxidative. By living inside the macrophages, it resists the humoral branch of the immune system [4].
Strategies to find a new compound from a natural source is carried out based on the tanimoto similarity indexing method and the most similar compound to the already established inhibitor of GP63 protein. Further, in silico studies provide an insight on the possible interaction of natural compounds with GP63 and CYP51 protein [5].

Protein structures and compound database
Protein sequence and structure was downloaded from Uniprot [6] and PDB database [7], respectively. The Uniprot id of Leishmanolysin protein is P23223 [6], the gene name is GP63, and the structure of the protein was modeled using RaptorX server [8] (Fig. 2a) and Lanosterol 14-alpha demethylase is C5WML [6], gene id is CYP51, pdb id of CYP51 is 3L4 D [9] (Fig. 2b), the and sequence and structure were downloaded from Uniprot website [6] and PDB website [7]. The inhibitors for GP63 protein were identified based on literature and other in silico methods.

Similarity indexing
Similarity indexing method using tanimoto clustering is one of the widely used method to calculate similarity coefficient based on atom fingerprints in compounds. The Similarity indexing was carried out using R programming [13] and ChemmineR package [14] using R studio interface.

Molecular docking
Molecular docking studies for Leishmanolysin (GP63) and Lanosterol 14-alpha demethylase (CYP51) with validated inhibitors and bioactive compounds of obtained from similarity indexing studies were carried out using extra precision glide docking (Glide XP) [15]. The 3D structures of compounds downloaded from PubChem and UNPD database were imported to Maestro Interface of Schrödinger and LigPrep [16] was used to understand ligand-protein interaction. Molecular docking studies were performed for the downloaded compounds. The prepared ligands were docked against protein GP63 and CYP51 (PDB ID: 3L4D) using Glide (grid-based ligand docking) program incorporated in the Schrodinger molecular modeling package by Maestro 11.0 [17]. GP63 and CYP51 are important in therapeutic impact of Leishmaniasis treatment. Benzyloxycarbonyl-Tyr-Leu-NHOH is a highly potent inhibitor that has been validated to inhibit the activity of GP63 thereby helping in therapeutic action in treatment of cancer. The crystal structure of GP63 was modeled using raptorX server using sequence downloaded from Uniprot and the structure of CYP51 was obtained from RCSB Protein Data Bank having X-ray diffraction resolution of 2.75 Å.
The 3D structures of the ligand for docking were downloaded in sdf format from PubChem and UNPD. The protein crystal structures were further optimized and minimized using protein preparation wizard using default settings for rectifying PDB structure for docking process. The ligands were optimized using the LigPrep which generate low energy 3D structure using the OPLS-2005 force field. Thereafter, molecular docking was performed with extra precision glide docking (Glide XP) to evaluate the binding affinity of these ligands towards the protein.

Molecular dynamics (MD)
MD simulations for protein-compound complex were conducted using Desmond's explicit solvent MD package [18] with built-in, optimized fluid simulation potential (OPLS 2005) force field. Protein Preparation Wizard (macro models), LigPrep (chemical molecules), and Epik (ligand protonation states) ensured the correctness of the chemical structures provided to Desmond. Further protein was subjected to restrained minimization for hydrogen atoms only by converging heavy atoms to RMSD value of 0.30 Å [19][20][21]. The system was set up for simulation using a predefined water model (simple point charge, SPC) as a solvent in a cubic box with periodic boundary conditions specifying box shape and size as 10 Å × 10 Å × 10 Å. Using the system-built option, the desirable electrically neutral simulation system was built in 10 Å buffer with 0.15 M NaCl (physiological concentration of monovalent ions). The system's relaxation was accomplished through the hybrid implementation of Steepest Descent and Broyden-Fletcher-Goldfarb-Shanno limited-memory algorithms [22]. Molecular dynamics for the system developed was carried out for 10 ns and the condition for the system was set at temperature 310 K with thermostat model and pressure of 1.01325 bar with ensemble class of NPT and energy of 1.2 at 100 ps recording interval.

Similarity indexing
The results from similarity indexing were obtained in terms of similarity coefficient and based on the similarity coefficient values, and 4 compounds from natural compound databases were selected having similarity more than 70%, i.e., similarity coefficient more than 0.7 ( Table 1). The source of natural compounds was identified to be from bacteria, fungi, and plant family based on the literature, and the 2D structures of the compounds are given in figures (Fig. 4a-d).

Molecular docking results
Molecular docking of inhibitor and bioactive compounds obtained from similarity indexing method was carried out using GLIDE tool in Schrodinger interface. The results from the docking study are represented by a score called docking score/glide score [27] which denotes the binding affinity. The lower the value the better the binding affinity of the compound with the target. In this study, the docking score of UNPD86509 from Streptomyces griseus showed better results in comparison with inhibitor 46173719 (Table 2 and Figs. 5 and 6). The interaction of compounds with the target GP63 can be seen in Fig. 5. Benzyloxycarbonyl-Tyr-Leu-NHOH (46173719) showed hydrogen bonding with target at ALA-43 and THR-34 position, and UNPD86509 showed hydrogen bonding with target at THR-34.

Molecular dynamics and simulations
The molecular dynamics for the GP63-UNPD86509 complex and CYP51-UNPD86509 complex carried out using Desmond showed a very stable configuration, and the complex showed a degree of affinity in terms of hydrogen bonding throughout the simulation time.

Simulation of GP63-UNPD86509 (natural compound) complex
The simulation showed very stable RMSD for both calpha (Protein RMSD Å) and Lig-fit (Ligand RMSD Å). The protein RMSD stayed in the range of 9.5-1.1 Å and the ligand RMSD stayed in the range of 7.5-10 Å. The protein ligand contacts showed a high degree of affinity with ALA-35 and ILE-39 by holding hydrogen bonding for most duration of 10 ns (Fig. 7). The orientation of ligand in the binding site of GP63 protein (Fig. 8a) and alignment of ligand molecules based on their states at the beginning of the simulation and after simulation are shown in Fig. 8b. Alignment of the ligand UNP86509 with pre-and post-simulation states showed that 36% of the change in orientation of movable bonds and this orientation was calculated using Discovery Studio Visualizer 2019 academic free version [28].

Simulation of CYP51-UNPD86509 (natural compound) complex
The simulation showed very stable RMSD for both calpha (Protein RMSD Å) and Lig-fit (Ligand RMSD Å). The protein RMSD stayed in the range of 2.1-2.8 Å and the ligand RMSD stayed in the range of 1.8-2.4 Å. The protein ligand contacts showed high degree of affinity  with TYR-102,TYR-115,MET-359, andMET-459by holding hydrogen bonding for complete duration of 10 ns (Fig. 9). The orientation of ligand in binding site of CYP51 protein (Fig. 10a) and alignment of ligand molecules based on their states at the beginning of simulation and after simulation are shown in Fig. 10b. Alignment of the ligand UNP86509 with pre-and postsimulation states showed that 38% of the change in orientation of movable bonds and this orientation was calculated using Discovery Studio Visualizer 2019 academic free version [28].

Discussion
GP63 and CYP51 proteins from Leishmania major are very important in the growth of the organism, and these targets have been used as one of the major therapeutic targets in the treatment of Leishmaniasis. Many compounds have been used to inhibit these targets but have failed to provide the needed impact in controlling the disease. Natural compounds have always been the source of drugs from a long time, and it is very important that we consider these natural compounds as very important source. Similarity indexing method provides fingerprint-based values to consider compounds having similarity with already established inhibitors. The detailed study of the compounds obtained from similarity indexing can help in finding a potential long-term cure for the Leishmaniasis. In silico analysis conducted in this study should serve as the basis for further wet lab validation of inhibition of GP63 and CYP51.

Conclusion
Natural compounds from different sources are always known to have answers to drug-like compounds; many compounds have been known to have anti-Leishmaniasis properties. Similarity indexing has given some of the compounds that can be drug-like molecules. The in silico study conducted here gives natural compounds based on fingerprints and further compares the in silico interaction of GP63 and CYP51 with a validated inhibitor and the bioactive natural compounds obtained from similarity indexing. The study showed very promising results in comparison of docking score for inhibitor benzyloxycarbonyl-Tyr-Leu-NHOH (46173719) with bioactive compound UNPD86509 from Streptomyces griseus. The docking score of interaction for CYP51 with 46173719 showed − 10.08 and docking score of CYP51with UNPD86509 showed − 11. 65. The binding energy of interaction for GP63 with 46173719 showed − 8.5 and docking score of GP63 with UNPD86509 showed − 7.8. Molecular dynamics studies also showed very stable configuration for CYP51-UNPD86509 and GP63-UNPD86509 complex in comparison with the inhibitor complex. From this study, we can conclude that the bioactive compounds obtained from similarity index can be a potential inhibitor for the CYP51 and GP63 proteins in Leishmania major which has great therapeutic role in the treatment of Leishmaniasis. The findings in this study should also encourage in vitro and in vivo studies evaluations of the same with the identified natural compounds.