iGEMDOCK: a graphical environment of enhancing GEMDOCK using pharmacological interactions and post-screening analysis

1Structure-based drug design is widely used to identify lead compounds with the growing availability of protein structures. Many tools (e.g., GEMDOCK, DOCK, AutoDock, and GOLD ) have been developed for virtual screening (VS) and successfully identified lead compounds for some target proteins. However, the accuracy of these docking tools remained intensive because of the incomplete understandings of ligand binding mechanisms. In addition, most of scoring functions are lack of pharmacological interactions that are essential for ligand binding or biological functions. Recently, some approaches have been proposed to derive pharmacological interactions from known compounds. These approaches apparently increase hit rates for identifying the active compounds which are often similar to the known compounds. In addition, these approaches are often unable to be applied for new targets, which have no known active compounds.

2Generally, a VS procedure consists of four main steps: preparations of the target protein and the compound library, docking and post-screening analysis (e.g., clustering compounds and pharmacological interactions). Most docking programs (e.g. DOCK and AutoDock) only provide docked poses or loosely coupled these steps. They often provided limit ability for post-screening analysis. Therefore, a VS framework, providing an easy-to-use graphic and integrated environment, is an emergent task for drug discovery.

3To address these issues, we have developed a structure-based VS framework, named iGEMDOCK, from preparations through to post-screening analysis. iGEMDOCK is an integrated environment, which integrates the heavily modified and enhanced in-house tool GEMDOCK, protein-ligand profiles, pharmacological interactions, and compound clusters. GEMDOCK was comparative to several docking tools (e.g. DOCK and GOLD) and has been successfully applied to identify new inhibitors and new binding sites for some targets. Notably, iGEMDOCK derives the pharmacological interactions from screening compounds without using a set of known active compounds. The pharmacological interactions, which often form binding pockets with specific physico-chemical properties of the target protein, can represent conserved interactions between the interacting residues and the screening compounds. We initially validated the pharmacological interactions on three therapeutic protein targets, including estrogen receptor α for antagonists and agonists and thymidine kinase. Our experimental results show that the derived pharmacological interactions are often essential for the ligand binding or maintaining biological functions for these targets. In addition, iGEMDOCK provided a post-screening analysis module, which is convenient for clustering compounds and visualizing the pharmacological interactions by interaction profiles. We believe that iGEMDOCK is useful for drug discovery and identifying essential residues and interactions for understanding the binding mechanisms.

4To initially validate the pharmacological interactions, we selected three therapeutic protein targets, including estrogen receptor α for agonists (ERA, PDB code 1gwr), estrogen receptor α for antagonists (ER, PDB code 3ert), and thymidine kinase (TK, PDB code 1kim) because these proteins were well studied. The catalytic mechanisms, biological functions, key functional residues, and active compounds of the three targets were available in the literatures. Estrogen receptor is an important therapeutic target for osteoporosis and breast cancer, and TK is a drug target for the therapy of herpes simplex virus type-1. Moreover, we also evaluate the docking and screening accuracy of iGEMDOCK. For docking, a highly diverse dataset comprising 305 protein-compound complexes (i.e., CCDC/Astex set) was selected; for screening, we prepared 10 known active compounds and 990 compounds were randomly selected from Available Chemical Directory (ACD) proposed by Bissantz et al. for each therapeutic protein target.

5iGEMDOCK is an integrated VS environment from preparations through post-screening analysis with pharmacological interactions (Fig. 1A). First, iGEMDOCK provides interactive interfaces to prepare both the binding site of the target protein and the screening compound library (Figs. 1B and 1C). Each compound in the library is then docked into the binding site by using the in-house docking tool GEMDOCK. Subsequently, iGEMDOCK generates protein-compound interaction profiles of electrostatic (E), hydrogen-bonding (H), and van der Waals (V) interactions. Based on these profiles and compound structures, iGEMDOCK infers the pharmacological interactions and clusters the screening compounds for the post-screening analysis (Figs. 1D and 1E). Finally, iGEMDOCK ranks and visualizes the screening compounds by combining the pharmacological interactions and energy-based scoring function of GEMDOCK.

6iGEMDOCK mines the pharmacological interactions based on protein-compound interaction profiles (Fig. 2). The size of each profile is N×2K where N and K are the numbers of screening compounds and interacting residues of the target protein, respectively. Here, an interacting residue is divided into two interacting groups: main and side chains. A profile P(I) with type I (E, H, or V) is given as (Fig. 2A):

7where pi,j is a binary value (0 or 1) for the compound i interacting to the residue group j. In the E and H profiles, the pi,j is set to 1 (green) if hydrogen-bonding or electrostatic interactions are yielded between the compound i and the residue j (energy ≤ -2.5 kcal/mol); otherwise, pi,j=0 (black). For the V profile, pi,j = 1 if the interacting energy is less than -4 kcal/mol (Fig. 2A).

8After the generations of the profiles, we identified the pharmacological interactions. For each interacting residue group, the z-score value is used to measure the interaction conservation between the interacting groups and the screening compounds. For computing the z-scores of interacting groups in a profile, 1000 randomly shuffled profiles are utilized to obtain the standard deviation (σ) and mean (μ). The z-score of the interacting residue group j is defined as, where fj is given as, where N is the number of screening compounds. Finally, we normalize the z-score value as follows:

9where Wj is the interaction conservation of the residue group j related to the largest z-score (zmax) among all of interacting groups in the binding site. Here, an interaction conservation is viewed as a pharmacological preference and an interaction is considered as the pharmacological interaction if Wj ≥0.4. For example, for the hydrogen profile of the target ERA, the pharmacological preferences of E353 and R394 are 0.64 and 0.80, respectively; for the V profile, the preferences of L387, L391, and F404 are 1.00, 0.61, and 0.90, respectively (Fig. 2B). In this case, over 300 (>30%) screening compounds form hydrogen bonds with the residues E353 or R394 by polar moieties (e.g., hydroxyl group (27%), carboxyl group (20%), sulfuric acid monoester (9%), ketone (8%), and phosphoric acid monoester (6%)). Moreover, the aromatic rings of the screening compounds are often sandwiched by vdW interacting residues L387, L391, and F404 (Fig. 2D).

10Based on the pharmacological interactions, we developed a pharmacological scoring function for identifying the active compounds from thousands of screening compounds. The pharmacological scoring function is given as

11where EGEMDOCK is the docked energy of GEMDOCK and E(E)pharma, E(H)pharma, and E(V)pharma are the pharmacological scores of electrostatics, hydrogen-bonding, and vdW interactions, respectively. The E(I)pharma with interaction type I (i.e., E, H, or V) is defined as

12where ej is the energy obtained by the GEMDOCK scoring function for the residue group j. Finally, iGEMDOCK provides the ranks of energy-based and pharmacological scoring functions for all screening compounds.

13iGEMDOCK is an easy-to-use VS environment and includes three main modules (Fig. 1): docking and virtual screening tool (GEMDOCK); post-screening analysis methods; and visualization tools (RasMol and Java Treeview). We employed ERA as an example to present these modules, procedures and features of iGEMDOCK.

14For protein-ligand docking/screening module, iGEMDOCK provides an interactive interface for the preparations of the binding site and compound library; setting docking parameters; and monitoring progress status (Fig. 1B). For most docking tools, users usually need to prepare the binding site structure and compound library through complicated steps (e.g., add hydrogen atoms and generate the grids of the protein). Here, iGEMDOCK provides a straightforward method to derive the binding site from the bounded ligand. For example, the binding site of ERA (PDB code 1gwr) was obtained from the estradiol (Fig. 1C). iGEMDOCK is able to automatically consider the effects of hydrogen atoms when preparing the binding site and the compound library. In addition, iGEMDOCK allows users to visualize and refine the binding site of the target protein. Additionally, iGEMDOCK offers the similar way to prepare screening compounds and docking parameters (e.g., the population size and the number of generations).

15After the screening process, iGEMDOCK utilizes the post-screening analysis module to infer pharmacological interactions and cluster screening compounds based on protein-ligand complexes and compound structures (Fig. 1D). First, iGEMDOCK generates interaction profiles and calculates the pharmacological preference (Wj) of each interacting group for deriving the pharmacological interactions (Fig. 2). These pharmacological preferences and interactions are shown in an interactive window (Fig. 2G); furthermore, RasMol displays the pharmacological interactions with conserved interacting residues and functional groups of compounds (Fig. 2H). Additionally, iGEMDOCK supports a hierarchical clustering method to cluster screening compounds according to interaction profiles and the atomic composition (Fig. 1E). The atomic composition, which is similar to the amino acid composition of a protein sequence, is useful for measuring compound similarity. iGEMDOCK provides an interactive interface for visualizing compound similarity with a hierarchical tree by Java Treeview. Finally, iGEMDOCK ranks and visualizes the screening compounds by combining the pharmacological interactions and the energy-based scoring function.

16The pharmacological interactions derived by iGEMDOCK are often involved in biological reactions or essential for ligand binding. We examined the pharmacological interactions on three selected target proteins (ERA, ER, and TK). First, we compared the pharmacological interactions, derived from 1000 screening compounds, to the consensus interactions, derived from 10 active compounds (Table 1 and Fig. 3). Here, the residue i is considered as "hot spot" if the consensus interaction ratio ≥0.5. Among 10 predicted pharmacological interactions (residues) for ERA, 9 pharmacological interactions (9 of 9 residues) agree with hot spots except the L387 with the hydrogen-bonding interaction. For TK, 8 of 14 pharmacological interactions (7 of 9 residues) are the hot spots. These results indicate the pharmacological interactions (residues) from screening compounds are often essential for the ligand binding. For example, 10 active compounds of TK form stacking interactions with the residue Y172 (vdW preference is 1.0 defined in Equation (1)) that stabilizes the binding of thymine or purine moieties.

17We also examined the pharmacological interactions by their biological functions or binding mechanisms. For estrogen receptor α, H524 (hydrogen-bonding preferences are 1.0 and 0.42 for ERA and ER, respectively) is involved in a hydrogen-bonding network; similarly, E353 and R394 (hydrogen-bonding preferences ≥ 0.5 for both ERA and ER) interact the structural water to form the hydrogen bonding network (Table 1 and Fig. 3). These two hydrogen bonding networks are essential for estrogen receptor modulators to trigger the responses of estrogen receptor α. For ER and ERA, hydrophobic interacting residues, L346, L387, F404, and L525 with high vdW interaction preferences, contact with the sterols or flavones scaffolds of the active compounds. These residues contribute the major vdW interactions for the ligand binding of estrogen receptor α.

18For TK, R222 and R163 play major roles for inhibitor and substrate binding, and their hydrogen-bonding preferences are 1.0 and 0.99, respectively (Table 1). Our method identified the electrostatic interactions of R222 and R163 (preferences are 1.0 and 0.4, respectively), which help to transfer phosphate in the substrate phosphorylation. However, these two electrostatic interactions are not observed from 10 active compounds (Fig. 3). For the residue Q125 (H preference 0.40), the TK activity was decreased over 90% if Q125 mutated into Asp, Glu, or Asn. The residues M128, Y172, H58, R163, and Y88 constitute a pocket to fix the substrate, and their vdW preferences are 0.58, 1.00, 0.68, 0.56, and 0.87, respectively (Table 1). For the substrate binding, M128 and Y172 sandwich the thymine moiety and W88 is a part of the quasi-helical motif. These results demonstrated that the pharmacological interactions derived by iGEMDOCK are often involved in the biological functions and the ligand binding.

19To initially evaluate the utility of iGEMDOCK for docking and virtual screening, we selected a highly diverse dataset with 305 protein-ligand complexes (i.e., CCDC/Astex set ) and ERA, ER, and TK with 1000 compounds as test sets. Please note that the docking and screening tool of iGEMDOCK is GEMDOCK which was well-studied for VS and some applications. In order to compare with previous works, we followed the docking procedures and performance indices proposed by Nissink, et al. A docked result was considered as a success solution if the root-mean-square derivation (RMSD) ≤2.0 Å between the docked solutions and X-ray crystal structures. For these 305 complexes, the success rates of iGEMDOCK and GOLD are 78% and 68%, respectively (Table S1 in additional file 1).

20The pharmacological scoring function was then applied to identify the active compounds from the 1000 compounds of ERA, ER, and TK. Furthermore, we compared the screening results with those of using the energy-based scoring function of GEMDOCK. These two approaches were tested on the same datasets. The true hits of the three testing sets were used to access the screening accuracy of the two approaches (Fig. S1 in additional file 1). The screening accuracy was generally improved when iGEMDOCK considered the pharmacological interactions.

21We compared iGEMDOCK (pharmacological scoring function) with three screening methods (DOCK, GOLD, and FlexX) on the ER and TK sets (Table S2 in additional file 1). To compare with previous works, we followed the experiment design and performance indices used by Bissantz et al. When true-positive rate is 80%, the false positive rates were 2.3% (iGEMDOCK), 13.3% (DOCK), 57.8% (FlexX), and 5.3% (GOLD), for ER. The false positive rates were 7.8% (iGEMDOCK), 23.4% (DOCK), 8.8% (FlexX), and 8.3% (GOLD) for TK.

22To identify leads from vast amount of docked poses generated during the virtual screening procedure is the key step for the drug discovery. In addition, the top-ranked compounds based on the scoring functions are not advisable since these compounds may be similar in structures or physico-chemical properties. For these two issues, iGEMDOCK provides a post-screening analysis module to cluster compounds based on the interactions profiles and the atomic compositions. Selecting representative compounds from each cluster is able to maintain compound diversity and then reduces the false positives. Further, when active compounds are available, users can choose the similar compounds in the same cluster of the actives based on hierarchical trees (Fig. 1E). Notably, iGEMDOCK visualizes the interaction profiles of the compounds, and thereby the top-ranked compounds with pharmacological interactions can be selected according to the interaction table (Fig. 2G).

23The post-screening analysis module of iGEMDOCK is useful for clustering and selecting compounds based on interaction profiles. We selected a set of compounds, including 10 ERA active compounds and top-ranked 100 compounds based on both the pharmacological and energy-based scoring functions. The hierarchical clustering dendrogram and the profile revealed that the protein-ligand interactions derived from the pharmacological scoring function are significantly different from those derived from the energy-based scoring function (Figs. 4A and B). The compounds with the high pharmacological scores and the active compounds consistently keep the pharmacological interactions; conversely, the compounds derived from the energy-based scoring function are often lack of the pharmacological interactions (Fig. 4B). This result indicates the pharmacological interactions are useful for identifying active compounds. For example, the two active compounds, ERA_03 and ERA_06, were ranked as 187 and 160 using the energy-based scoring function, respectively. When the pharmacological interactions were considered, the ranks of ERA_03 and ERA_06 were 91 and 87, respectively (Fig. 4C).

24Some compounds having high pharmacological scores are structurally and chemically similar to the active compounds (Fig. 4C). For example, MFCD00003670 (Tetrahydrocortisol) and MFCD00012742 (Tetrahydropapaveroline) were analogues of the ERA active compounds, and both of them were clustered into the same cluster. In addition, the pharmacological scoring function can reduce the ill-effect of most energy-based scoring functions which often favor high molecular weight or highly-polar compounds. For instance, the ranks of two high molecular weight and polar compounds, MFCD00009772 and MFCD00016940, are 267 and 274, respectively (Fig. 4C). To further examine the pharmacological scoring function, we analyzed the relationship between the molecular weights of the active compounds and the rank improvement using the pharmacological scoring function (Fig. 5). The pharmacological scoring function generally improves the screening accuracy when the molecular weights of the active compounds are less than 400.

25In summary, iGEMDOCK can mine the pharmacological interactions from the screening compounds without known active compounds. Therefore, iGEMDOCK can provide a good starting point for deriving pharmacological interactions (residues) and identifying new potential active compounds for a new protein structure. In addition, iGEMDOCK offers the visualization of the interaction profiles, pharmacological interactions, and the hierarchical clustering dendrogram. Users are able to easily observe and select compounds for post-screening analysis to enrich accuracies.

26We have developed a structure-based VS framework iGEMDOCK from the preparations through to the post-screening analysis. iGEMDOCK is an integrated and easy-to-use environment which is especially useful for post-screening analysis and inferring pharmacological interactions from screening compounds. The friendly user interface is helpful to biologically oriented nonexperts. The experimental results show that the pharmacological interactions are often essential for the binding of the active compounds and involved in biological mechanisms. The pharmacological interactions can reduce the ill effects of energy-based scoring functions to enrich the hit rates in VS. We believe iGEMDOCK is useful for drug discovery and understanding protein-ligand mechanisms.

27The authors declare that they have no competing interests.

28KCH, YFC, and JMY conceived and designed the experiments. YFC, SRL, and JMY implemented the program. KCH, YFC, and JMY performed the experiments and analyzed the data. KCH, YFC, and JMY wrote the paper.