The development and implementation of a range of molecular docking algorithms, based on different search methods (Taylor et al., 2002, Halperin et al., 2002) was observed in the last few years. This approach has had several recent successes in drug discovery (Sechi et al., 2005; Liu et al., 2005).

In the field of molecular modeling, docking is a method which predicts the preferred orientation of one molecule to a second when bound to each other to form a stable complex (Lengauer and Rarey, 1996). Knowledge of the preferred orientation in turn may be used to predict the strength of association or binding affinity between two molecules using for example scoring functions.

Docking is frequently used to predict the binding orientation of small molecule drug candidates to their protein targets in order to in turn predict the affinity and activity of the small molecule. Hence docking plays an important role in the rational design of drugs (Kitchen et al., 2004). Given the biological and pharmaceutical significance of molecular docking, considerable efforts have been directed towards improving the methods used to predict docking.

Evaluation of existing docking algorithms can assist in the choice of the must suitable docking programs for any particular study. Effectively, several studies estimating and comparing the accuracies of protein-ligand programs like Dock, ICM, Gold have been reported (Perola et al., 2004 ; Bursulaya et al., 2003).

The goal of this study was to evaluate the ability of ArgusLab, a freely distributed molecular modeling package in which molecular docking is implemented, to reproduce crystallographic binding orientations and to compare its accuracy with that of the widely well established docking packages, Autodock and FlexX.

ArgusLab4.0 has fast become a favorite introductory molecular modeling package with academics mainly because of its user-friendly interface and intuitive calculation menus (Thompson, 2004). The ArgusDock docking engine, implemented in ArgusLab, approximates an exhaustive search method. Flexible ligand docking is possible with ArgusLab, where the ligand is described as a torsion tree and grids are constructed that overlay the binding site. Ligand’s root node (group of bonded atoms that do not have rotatable bonds) is placed on a search point in the binding site and a set of diverse and energetically favorable rotations is created. For each rotation, torsions in breadth-first order are constructed and those poses that survive the torsion search are scored. The N-lowest energy poses are retained and the final set of poses undergoes coarse minimization, re-clustering and ranking.

AutoDock3.0 explores the conformational space of the ligand using the Lamarkian genetic algorithm (LGA), which is a hybrid of a genetic algorithm (GA) with an adaptive local search (LS) method (Morris et al., 1998). In this approach, the ligand’s state is represented as a chromosome, which is composed of a string of real-valued genes describing the ligand location (three coordinates), orientation (four quaternions) and conformation (one value for each torsion). The simulation is started by creating a random population of individuals. It is followed by a specified number of generation cycles, each consisting of the following steps: mapping and fitness evaluation, selection, crossover, mutation and elitist selection. Each generation cycle is followed by a local search. The solutions are scored using an energy-based scoring function, which includes terms accounting for short-ranged Van Der Waals and electrostatic interactions, loss of entropy upon ligand binding, hydrogen bonding and solvation.

FlexX1.11 (Rarey et al., 1996; Kramer et al., 1999) employs an incremental reconstruction algorithm. In this algorithm rigid base fragments are identified first. At the next step, the selected fragment is placed into the active site of the receptor using a hashing technique. The complete ligand is constructed by adding the remaining components step by step. At each step of reconstruction a specified number of optimal partial solutions are selected for the next extension step. The scoring is done using a modified Böhm scoring function, which includes the following terms: entropic, which accounts for loss of entropy upon ligand binding; hydrogen bonding; ionic, accounting for electrostatic interactions; aromatic, which accounts for interactions between aromatic groups; and lipophilic, which accounts for hydrophobic interactions. All terms, except the entropic term, are scaled by a corresponding heuristic distance and an angle dependent penalizing function.

In all algorithms studied here, the receptor is treated as a rigid body and a grid potential is used to evaluate the scoring functions. This simplification allows one to perform docking more efficiently, which is especially crucial in database screening. Arguslab requires a PDB format file for both ligand and receptor. The binding site was defined from the coordinates of the ligand in the original PDB file. Argusdock exhaustive search docking engine was used, with grid resolution of 0.40 Å. Docking precision was set to ‘high precision’ and ‘flexible ligand docking’ mode was employed for each docking run.

AutoDock requires the receptor and ligand coordinates in MOL2 format. Nonpolar hydrogen atoms were removed from the receptor file and their partial charges were added to the corresponding carbon atoms. The program Mol2topdbqs was used to transform the receptor MOL2 file into the PDBQS format file containing the receptor atom coordinates, partial charges and solvation parameters. The program AutoTors was used to transform the ligand MOL2 file into a PDBQ file, merge nonpolar hydrogen atoms and define torsions. The grid calculations were set up with the utility Mkgpf3 and maps were calculated with the program AutoGrid. The grid maps were centered on the ligand's binding site and were of dimension 61× 61 × 61 points. The grid spacing was 0.375 Å yielding a receptor model that included atoms within 22.9 Å of the reference binding site center. The default parameter settings generated by the program Mkdpf3 were used for docking. For each complex 10 dockings were performed. The initial population was set to 50 individuals; maximum number of energy evaluations was 2.5×105; maximum number of generations was 27,000. The other parameters provided by the default setting were the same as in the followed reference (Morris et al., 1998).

FlexX requires a MOL2 format file for the ligand and a PDB format file for the receptor. The default settings as provided with the FlexX package were used for flexible docking and database screening. The conformational flexibility of the ligand is modeled by a discrete set of preferred torsional angles for acyclic single bonds. The rings were considered rigid, since the program CORINA for treating multiple conformations of the rings was not included in the distribution. The active site and the interaction surface of the receptor were defined by using a reference ligand and a 6.5 Å cutoff distance. Base fragments were selected automatically. The maximum number of base fragments was 4. The base fragment was placed into the active site using two algorithms. The first one superimposes triples of interaction centers of a base fragment with triples of compatible interactions in the active site. The second algorithm, called matching, is used when the base fragment had fewer than three interaction centers. The sampling was done with 400 solutions per partial solution at each iteration of incremental construction.

The best ranking poses predicted by the three programs Arguslab, Autodock and FlexX are shown in the Figure 1 and their root mean square deviation (RMSD) values from the original crystallographic pose determined. It can be observed that both Autodock and FlexX outerform Arguslab. For RMSD interval < 2Å, the difference in docking accuracies between the three programs is so important but decrease significantly in RMSD interval < 3Å.

Figure 1: Best pose with reference to crystallographic pose

Figure 2 shows the evaluation of the docking algorithms for their sampling accuracy. The percentage of poses with RMSD within 2Å from the experimental structure was 65% for Autodock, 55% for FlexX and only 30% for ArgusLab. This confirms the results reported by earlier studies, Autodock seems to be highly efficient in terms of sampling (Badry et al., 2003). However, under less rigorous conditions, the performance of ArgusLab is hugely improved with 64% of the top ten poses falling within 3Å of the crystallographic pose. This signifies that ArgusLab still gives some biological results and can be used in educational demonstrations.

Figure 2: Top ten poses with reference to crystallographic pose

The effect of a ligand parameter on docking accuracy is another kind of analysis we have carried out (figure 3). It is a well-known fact that as the number of rotatable bonds of the ligand increases, the docking accuracy falls since a much larger conformational space has to be sampled. The complexes in the present study were divided into three groups, ligands with 1 to <10 rotatable bonds, ligands with 11 to < 15 rotatable bonds and those with > 15 rotatable bonds. The results confirm earlier works. Indeed for all algorithms, the docking accuracy decreases when the number of rotatable bonds increases. Also in all cases, accuracy of both Autodock and FlexX is approximately double that of ArgusLab. This decrease is very pronounced when the number of rotatable bonds exceed 15. Though, an essential remark is that docking time in both ArgusLab and FlexX is typically much shorter than that of Autodock.

Figure 3: Ligand rotatable bonds in relation to docking accuracy.

To further evaluate these docking programs, another test we have conducted is to study the chemical nature of their protein-ligand interactions and then to check the success rate of each scoring function (Figure 4). The classification is aided by using X-Score. For any given proteinligand complex, if the contribution of the H-bond term in XScore is 50% larger than the hydrophobic term, it is classified as the "hydrophilic" type. If the contribution of the hydrophobic term is 50% larger than the H-bond term, it is classified as the "hydrophobic" type. Otherwise, the complex is considered to have mixed hydrophilic and hydrophobic factors in the protein-ligand interaction and thus is classified as the "mixed" type. We have used X-Score for this classification process because it is the only one with open source codes, so we can analyze the hydrophobic and the hydrophilic terms conveniently. Best results for hydrogen bond driven complexes are given by FlexX (73%) and for hydrophobic-burial driven ones are given by Autodock (67%). There is no perceptible change in the docking accuracy of A rgusLab with degree of hydrogen bonding. Studies for determination of IC50 and MIC, in specialized laboratory, are needed to confirm these in silico results.

Figure 4: % of hydrogen bonding in terms of docking accuracy.

Our results suggest that the two docking programs Autodock and FlexX do a reasonable job in docking and should aid significantly the drug discovery process. However, Autodock outperforms the two other programs and its use for molecular docking seems to be most advantageous. This study shows that commercial packages surpass the freely available docking program in all parameters experienced. The study also revealed that, in less rigorous conditions, ArgusLab can be used for demonstration of molecular docking method to novices in this area owing to its easiness to use graphical user interface. Moreover, some future advances can be made in this program at the expense of the docking time.

Understanding the ruling principles whereby protein receptors recognize, interact, and associate with molecular substrates and inhibitors is of paramount importance in drug discovery efforts. Protein-ligand docking aims to predict and rank the structures arising from the association between a given ligand and a target protein of known 3D structure. Despite the breathtaking advances in the field over the last decades and the widespread application of docking methods, several downsides still exist, in particular, protein flexibility. Indeed, a critical aspect for a thorough understanding of the principles that guide ligand binding in proteins is a major hurdle in current protein-ligand docking efforts and needs to be more efficiently accounted for. In the future the key concepts of protein-ligand docking methods will be outlined, with major emphasis being given to the general strengths and weaknesses that presently characterize this methodology.

The authors would like to thank Pr David Perahia for the AutoDock 3.0 academic license.