Molecular docking studies for the identification of novel melatoninergic inhibitors for acetylserotonin-O-methyltransferase using different docking routines

Background N-Acetylserotonin O-methyltransferase (ASMT) is an enzyme which by converting nor-melatonin to melatonin catalyzes the final reaction in melatonin biosynthesis in tryptophan metabolism pathway. High Expression of ASMT gene is evident in PPTs. The presence of abnormally high levels of ASMT in pineal gland could serve as an indication of the existence of pineal parenchymal tumors (PPTs) in the brain (J Neuropathol Exp Neurol 65: 675–684, 2006). Different levels of melatonin are used as a trait marker for prescribing the mood disorders e.g. Seasonal affective disorder, bipolar disorder, or major depressive disorder. In addition, melatonin levels can also be used to calculate the severity of a patient’s illness at a given point in time. Methods Seventy three melatoninergic inhibitors were docked with acetylserotonin-O-methyltransferase in order to identify the potent inhibitor against the enzyme. The chemical nature of the protein and ligands greatly influence the performance of docking routines. Keeping this fact in view, critical evaluation of the performance of four different commonly used docking routines: AutoDock/Vina, GOLD, FlexX and FRED were performed. An evaluation criterion was based on the binding affinities/docking scores and experimental bioactivities. Results and conclusion Results indicated that both hydrogen bonding and hydrophobic interactions contributed significantly for its ligand binding and the compound selected as potent inhibitor is having minimum binding affinity, maximum GoldScore and minimum FlexX energy. The correlation value of r2 = 0. 66 may be useful in the selection of correct docked complexes based on the energy without having prior knowledge of the active site. This may lead to further understanding of structures, their reliability and Biomolecular activity especially in connection with bipolar disorders.

tryptophan metabolism) [1]. ASMT lies under three critically important sub-classes of enzyme, transferases, one-carbon group transferases, and methyltransferases. ASMT protein is found in both prokaryotes and eukaryotes. In humans, ASMT gene encoded this enzyme. There are two identical copies of ASMT (one on the X arm and another on the Y arm of chromosome) [2,3]. High Expression of ASMT gene is evident in pineal parenchymal tumors (PPTs). The presence of abnormally high levels of ASMT in pineal gland could serve as an indication of the existence of PPTs in the brain [4]. Different levels of melatonin are used as a trait marker for prescribing the mood disorders e.g. seasonal affective disorder, bipolar disorder, or major depressive disorder. In addition, melatonin levels can also be used to calculate the severity of a patient's illness at a given point in time. As studies have revealed a direct correlation between the amount of ASMT in the pineal gland and the melatonin level, additional knowledge of ASMT could provide valuable insight into the nature and onset of these impairing disorders [5]. In studying the various properties associated with protein-ligand interactions, docking is a powerful tool. Since molecules in nature have a tendency to be found in their lowest energy form, the final configuration should also be of low energy [6]. Understanding these properties is crucial in the rational design of potent inhibitors.
Molecular Docking is an effective and competent tool for in silico screening. It is playing an important and ever increasing role in rational drug design [7,8]. Docking is a computational procedure of searching for an appropriate ligand that fits both energetically and geometrically the protein's binding site. In other words, it is a study of how two or more molecules e.g. ligand and protein, fit together. The problem is like solving a 3D puzzle [9]. During the past decade, for understanding the formation of intermolecular complexes, the application of computational methods in this arena has been subjected to intensive research. It is commonly known that molecular binding of one molecule (the ligand) to the pocket of another molecule (the receptor), which is commonly a protein, is responsible for accurate drug activity. Molecular docking has been proved very efficient tool for novel drug discovery for targeting protein. Among different types of docking, protein-ligand docking is of special interest, because of its application in medicine industry [10]. Proteinligand docking refers to search for the accurate ligand conformations within a targeted protein when the structure of proteins is known [11].
Docking procedures are basically the combination of search algorithms and scoring function. The largest number of search algorithms and scoring functions are available. Search algorithms predict the ligand binding orientation and conformations commonly referred to as posing [11]. Some common search algorithms are [9]: Monte Carlo methods, Genetic algorithms, Fragment-based methods, Point complementary methods, Distance geometry methods, Tabu searcher and Systematic searches. In order to differentiate between the active and random compounds, the scoring functions are employed. The scoring functions predict binding free energies in ligand-protein docking generally in 7-10 kJ/mol [12]. Numbers of molecular docking software are employed in drug research industry [9]. The most popular and commonly used softwares for molecular docking are AutoDock [13][14][15], AutoDock/Vina [16], GOLD [17,18], FlexX [19], FRED [20], DOCK [21] and ICM [22]. For docking purpose, Auto-Dock/Vina employed Broyden-Fletcher-Goldfarb-Shanno algorithm and it significantly improves the average accuracy of the binding mode predictions compared to AutoDock 4 [16]. FlexX employed an IC algorithm. IC algorithm attempts to reconstruct the bound ligand by first placing a rigid anchor in the binding site and later using a greedy algorithm to add fragments and complete the ligand structure. GOLD considers the degree of freedom in the binding site that corresponds to reorientation of hydrogen bond donor and acceptor groups. This degree of freedom represents only a very small fraction of the total conformational space that is available but should account for a significant difference in binding energy values [23].
In connection with efforts rendered in searching for novel inhibitors of ASMT, we perform a comparative docking study with four extensively used programs: AutoDock/Vina, GOLD, FlexX and FRED. The docking accuracy and scoring reliability of the selected docking approaches were evaluated by docking seventy three melatoninergic ligands with ASMT and correlating the predicted binding affinities with the experimental values.

ASMT and melatoninergic inhibitors
The protein used in the docking study was obtained through homology modeling by Azam et al., [24]. Dogsite web server was employed to detect the binding pocket of ASMT (Table 1) [25]. Seventy three structurally diverse ASMT inhibitors (Additional file 1) with representative good biological activity were selected from the literature [26][27][28][29][30][31]. The 2D structures of the melatoninergic inhibitors were drawn using chemical  [32]. The conformational energies of inhibitors were minimized by using UCSF Chimera [33]. The minimized structures were then subjected to docking studies.

Docking protocol
Molecular docking protocols are widely used for predicting the binding affinities for a number of ligands. In current work, our aim was to examine the possibility of an existing relationship between the experimental bioactivities of the inhibitors under study and the docking scores. In order to get accurate results, all the docking experiments were performed with the default parameters. The time to dock one ligand was approximately 1-2 min. Docking with AutoDock/Vina, GOLD and FRED was performed on a Linux workstation (openSUSE11.4) with an Intel Pentium D processor (3.0 GHz) and 1 GB of RAM where as FlexX was run on windows 7 equipped with an Intel® Atom™ processor (1.67 GHz) and 1GB of RAM.

Docking using AutoDock/Vina
Intermediary steps, such as pdbqt files for protein and ligands preparation and grid box creation were completed using Graphical User Interface program AutoDock Tools (ADT). ADT assigned polar hydrogens, united atom Kollman charges, solvation parameters and fragmental volumes to the protein. AutoDock saved the prepared file in PDBQT format. AutoGrid was used for the preparation of the grid map using a grid box. The grid size was set to 60 × 60 × 60 xyz points with grid spacing of 0.375 Å and grid center was designated at dimensions (x, y, and z): -1.095, -1.554 and 3.894. A scoring grid is calculated from the ligand structure to minimize the computation time. AutoDock/Vina was employed for docking using protein and ligand information along with grid box properties in the configuration file. AutoDock/Vina employs iterated local search global optimizer [34,35]. During the docking procedure, both the protein and ligands are considered as rigid. The results less than 1.0 Å in positional root-meansquare deviation (RMSD) was clustered together and represented by the result with the most favorable free energy of binding. The pose with lowest energy of binding or binding affinity was extracted and aligned with receptor structure for further analysis.

Docking using GOLD (Genetic Optimization for Ligand Docking)
GOLD utilizes genetic algorithm to explore the rotational flexibility of receptor hydrogens and ligand conformational flexibility [18]. In GOLD docking was carried out using the wizard with default parameters population size (100); selection-pressure (1.1); number of operations (10,000); number of islands (1); niche size (2); and operator weights for migrate (0), mutate (100), and crossover (100) were applied. The active site with a 10 Å radius sphere was defined by selecting an active site residue of protein. Default Genetic Algorithm settings were used for all calculations and a set of 10 solutions were saved for each ligand. GOLD was used by a GoldScore fitness function. GoldScore is a molecular mechanism like function and has been optimized for the calculation of binding positions of ligand. It takes into account four terms: Where S hb_ext is the protein-ligand hydrogen bonding and S vdw_ext are the van der waals interactions between protein and ligand. S hb_int are the intramolecular hydrophobic interactions whereas S vdw_int is the contribution due to intramolecular strain in the ligand.

Docking using FlexX
FlexX (which is now a part of LeadIT) is a flexible docking method that uses an Incremental Construction (IC) algorithm and a pure empirical scoring function similar to the one developed by Böhm and coworkers [36] to place ligands into the active site. IC algorithms first dissect each molecule into a set of rigid fragments according to rotatable bonds, and then incrementally assemble the fragments around the binding pocket [19]. For docking studies, the pdb files of ligands were transformed into a SYBYL mol2 file format and a ligands library was generated. A receptor description file was prepared through the FlexX graphic interface. An active site was defined by selecting the residue of the protein. The active site includes protein residues around 10 Å radius sphere centered on the center of mass of the ligand. Based on energy values, top ten ranked poses for each ligand in data set were selected for further analysis.
The free binding energy ΔG of the protein-ligand complex is given by: Here, f (ΔR, Δα) is a scaling function penalizing deviations from the ideal geometry and N rot is the number of free rotatable bonds that are immobilized in the complex. The terms ΔG hb , ΔG io, ΔG ar and ΔG 0 are adjustable parameters. ΔG lipo is lipophilic contact energy (Rarey et al., [19]).

Docking with FRED (Fast Rigid Exhaustive Docking)
FRED uses multi-conformer docking algorithm which separately generates a set of lowenergy conformers, and then do rigid docking for each conformer [37]. In order to carry out correct docking, FRED required accurately prepared receptor file as well as a ligand conformer library. The receptor file was prepared by using make-receptor file provided in FRED whereas ligand conformer library was created in Omega 2.3.2 (Open-Eye Scientific Software) with default settings. The volume of the docking box centered on the receptor was expanded in all directions until it was approximately 31671 Å 3 . The dimensions of the box were: 28.10 Å × 32.91 Å × 34.25 Å. FRED with a Gaussian type fitting scoring function Chemgauss4 was used to dock ASMT with ligands conformer library in order to obtain a potent inhibitor against ASMT. Chemgauss4 uses the potentials between the chemically matched positions around the ligand docked pose. Those chemical positions are complementary to the nearby specific groups in the receptor. Generally, the interactions are either hydrogen bond donors or acceptors and a favorable hydrogen bond score is obtained when a polar hydrogen position on one molecule overlaps a lone pair position on another molecule. The interactions which can be scored by Chemgauss functions are: steric, acceptor, donors, coordinating groups, metals, lone pairs, polar hydrogens and chelator coordinating groups [38].

Results and discussion
In order to recognize an accurate docking routine for carrying out molecular docking studies of ASMT protein as well as to identify the potent inhibitors against that protein, four widely used docking routines (AutoDock/Vina, GOLD, FlexX and FRED) are compared in this work. Each docking routine returned top ten ranked docked poses for each ligand. The number and categories in which the dock poses fall are summarized in Table 2. Among all the studied docking routines, AutoDock/Vina was found to be the best for carrying out blind docking and in generating poses that bind best deep inside the 5 Å of the binding pocket. It generated 70% good poses. Gold and FRED also performed well with 65% and 45% of good poses respectively. FlexX returned no significant results in generating accurate poses. This variation might be because of the algorithms employed by the routines, grid box specification and active site residue specification. Overall results showed relatively poor performance ranging from 41 to 70%. These percentages were surprisingly low, indicating the docking programs often failed to find the correct binding mode. The percentages of docked poses for ASMT melatoninergic inhibitors docked complexes obtained by the different docking routines are shown in Figure 2.
Evaluation of docking accuracy of docking programs requires the programs run at approximately comparable speeds. The average time required for carrying out the docking calculation using different scoring functions and docking routines of a single ligand is shown in Table 3. According to Table 3 in the context of the average time required for docking a single ligand, FRED is the fastest algorithm with average 1.4 sec, followed by Gold and AutoDock/Vina with an average time of 1.66 and 2-3 sec respectively. FlexX seemed to be the more time-consuming approach.  The correlation between the experimental bioactivities (PIC50) and binding affinities and the calculated scores generated by each docking routine for ASMT docked complexes are shown in Figure 3. The best squared correlation coefficient of r 2 = 0. 66 was observed between binding affinities (AutoDock/Vina) and experimental values (Figure 3a). Low values of correlation for the Gold scores (r 2 = 0.37) and FlexX scores (r 2 = 0. 28) was obtained. FRED (Chemgauss4), due to its rigid-body approach, presented the lowest value of the squared correlation coefficient among the tested routines ( Figure 3d).
The 2D view of protein-ligand interactions of the best poses generated by all the four studied routines are shown in Figure 4. As clearly depicted in Figure all molecules exhibit the same binding mode. Interestingly, important interactions can be found between these atoms and the residues Ser213, Ser98, Val97, Thr100, Val211, Ser227, Arg210, Arg280, Phe212, Leu198, Ile198, Ser104, Thr195, Leu160, Tyr327, Tyr108, Trp117, Leu326, Phe29, Phe19, Gln334, Asn330, Ala159, Lys107, Met105 which directly participate in the catalytic mechanism of this enzyme. The ligand-enzyme complex is stabilized mainly by hydrogen bonds and hydrophobic interactions. All the top docked poses generated by each docking routine exhibited well-established bonds with   For B22 the binding affinity is found to be −9.2 Kcal/mol. The orientation and hydrogen bonding, ionic interactions of B22 within ASMT active site are shown in Figure 5a.
Different sets of hydrogen bonding interactions with polar side chain residues of Arg280, Ser104, Thr100, Tyr108 and with phenol of Tyr327 are observed at distances within 4 Å. An ionic interaction with the side chain residue of Arg280 is also observed. The indole ring and other side chain carbons of B22 formed strong hydrophobic interactions with nonpolar residues Leu160, Tyr327, Tyr108, Leu326, Trp117 and Phe26 and are suggested to increase binding affinity (Table 4). Figure 5b shows the binding mode of top pose B22ASMT complex generated by GOLD with GoldScore of 64.88. Compound B22 was mediating hydrogen bond interactions with the side chain residues of Arg280, Ala159, Tyr327 and Tyr108. The phenol and indole ring of B22 formed favorable hydrophobic contacts with Tyr108, Tyr327, Tyr366, Trp117 and Leu326 and ionic interaction with Arg280.
The FlexX generated 10 solutions for B22. The highest-ranking solution has a binding energy of −25.45 kJ/mol. Hydrogen bonds with a backbone and side chain residues of Arg210, Lys223 and Thr207, back bone residues of Val211, back bone residues of Phe212 and Glu224 are observed. The strong hydrophobic interactions were with Thr195, Lys223, Leu228, Arg169, pyrrole ring of His209, Val211 and aromatic ring of Phe212. Ionic interactions were with Arg210 and Glu224 (Figure 5c).
The top ranked pose of ASMT-A3 generated by FRED has Chemgauss4 score of −9.87. Phenolic side chain of Tyr327, Tyr108, amide side chain of Asn330 and side chain residues of Lys107 formed strong hydrogen bonding interactions with electronegative atoms of ligand. Phenol group of Tyr108, and the aromatic ring of Trp117 and Tyr 327 were highlighted as major contributor of hydrophobic interactions. A3 also resulted in favorable ionic interactions through the active site with the back bone residues of Asn 330 and    side chain residues of Lys107 ( Figure 5d) ( Table 4). The lowest binding affinity, high Gold-Score, a low binding energy and number of observed hydrogen bonding interactions between the backbone residues and the various important residues at the entrance of the pocket enables B22 to be a strongly binding inhibitor of ASMT. As shown in Table 5, among the provided data set AutoDock/Vina docked all the seventy three inhibitors with binding affinities in a convincing range. Like AutoDock/Vina, GOLD also docked all the provided inhibitors however GOLD tends to be limited in terms of ligand flexibility. Kellenberger reported that out of 100 protein-ligand complexes, GOLD successfully docked 80% of the ligands at an RMSD value of 2.0 Å [37]. Likewise Perole [39] while doing virtual screening of HIV-1 protease established that when looking at only the top 10% of the rankings, GOLD will accurately bind 60% of the ligands. Additionally GOLD depends more on the binding site and it generates worst results when hydrophobic interactions plays an important role in binding [39]. The performance of FRED with respect to ranking the inhibitors in the list was poor. FRED ranked inhibitor B22 on the 11 th place whereas all the remaining routines ranked B22 in 1 st place. FRED docked fifty three inhibitors out of seventy three provided inhibitors which is a low number as compared to other routines. Kellenberger et al., concluded in their work that FRED had problems with small, polar and buried ligands [37]. FlexX docked sixty eight inhibitors out of seventy three. The most severe restriction of FlexX is that it treats receptor as a rigid entity. Leach in 1994 reported an approach handling receptor flexibility [40].
In the light of the above analysis, the B22ASMT docked pose generated by Auto-Dock/Vina produced the best results. It forms hydrogen bonds, hydrophobic and ionic interactions with the important residues of the binding pocket of ASMT thus stabilizing the structure of target receptor. The dock pose with least binding energy has the highest affinity and hence is the best docked conformation.

Conclusion
Docking and scoring have evolved significantly over the past years. It has become a valuable tool in drug discovery process. Our goal of this study was to explore the feasibility of four different docking approaches: (AutoDock/Vina, GOLD, FRED and FlexX) for our target ASMT and to find out the lead compound. We compared the predictive power of each docking and scoring function. Our results suggest that all docking programs studied here do a reasonable job in docking and should aid significantly the drug discovery process. However, AutoDock/Vina consistently outperformed as compared to other programs and was found to be relatively more useful in blind docking pose prediction. Moreover, analysis of the docked ligands with the protein brought into focus some important interactions operating at the molecular level. The results of the ligand docking showed that the binding pocket involves the amino acid residues Ser213, Ser98, Val97, Thr100, Val211, Ser227, Arg210, Arg280, Phe212, Leu198, Ile198, Ser104, Thr195, Leu160, Tyr327, Tyr108, Trp117, Leu326, Phe29, Phe19, Gln334, Asn330, Ala159, Lys107, Met105. The important hydrogen bond forming amino acid residues was Arg280, Thr100, Tyr108, Asn330, Trp117 and Tyr327. In conclusion we have discovered a highly potent lead compound which will be useful for the design of novel less toxic and highly efficient drug for the treatment of bipolar disorders and PPTs.

Additional file
Additional file 1: 2D view of melatoninergic inhibitors.