The Performance of Several Docking Programs at Reproducing Protein–Macrolide-Like Crystal Structures

1Many important drugs and hundreds of natural products with promising bioactivities are macrolides or have macrolide-related structures. Such compounds often show favorable pharmacokinetic profiles with regard to open-chain molecules, despite not always following Lipinski’s rules. Conformational search is important for identifying their bioactive conformers, which is a prerequisite for good results in a wide range of modeling methods. Indeed, adequate sampling is extremely important in docking calculations and small changes in the search parameters can have a significant effect on their performance.

2In this context, over the past two decades our research group has worked on the total synthesis of bioactive macrolide-like natural products, such as fluvirucins; amphidinolides X, Y, K, E and B2; and palmerolide A. We performed docking calculations and collaborated with biochemists, chiefly with the aims of elucidating the mechanisms of action of some cytotoxic macrolides. Since then, updated and new docking software has been released. This led us to wonder whether some explanations and predictions, included in PhD theses or articles in the past decade, would still stand in the present scenario. For example, an updated version of AutoDock 3.0.5 (henceforward AD 3.0, a popular program that we used years ago) is AutoDock 4.2.6 (henceforward AD 4.2). A faster program, AutoDock Vina 1.1.2 (henceforward AD Vina) was more recently launched by the same institution. It is expected that these programs will provide improved results and allow better explanations than those facilitated by AD 3.0. Last year, once most of our AD Vina results were at hand, PSOVina, a program that includes a particle swarm optimization algorithm, was released. PSOVina reduces the calculation time of AD Vina even further; since the accuracy is the same, only a few tests have been performed using it. All these packages are free. Other docking methods that we have examined are DOCK 6.5 (henceforward, DOCK) and Glide 6.6 (Glide). DOCK is also free. Glide, from Schrödinger, appears to be the most popular of the commercial docking programs, according to an exhaustive SciFinder search. Thus, in the present study, we restricted ourselves to these five programs. Comparative studies of different docking programs have been reported involving small molecules (not discussed here for the sake of space restrictions). However, we have found only two precedents of the interactions of macrolides and macrolactams with their receptors: the docking of dictyostatin into the taxol binding site by Glide, AutoDock 4.0 and RosettaLigand; and that of spinosyn precursors (22-membered macrolides) and their products by AD Vina and DOCK.

3Our objective was to evaluate how reliably these programs reproduce a set of protein–macrolide complexes for which crystallographic structures were reported with a good resolution (there are only eight disparate representative examples to date). Since this number is low, in our set we included some macrolactams, cyclodepsipeptides, and other heteroatom-containing macrocycles of similar size to our desired macrolides for which high-resolution crystallographic structures were also published. This gave us a total of 20 examples, which we ordered according to their PDB code (Figure 1). Six of these complexes (1NTI, 1UU3, 1W96, 2IYA, 2XBK, and 3EKS) have just been examined using other methods. We envisage expanding our test to other programs in the future and applying the most efficient to dock macrolide-like anti-tumor agents synthesized by our research group into cytoskeleton proteins.

4First of all, a self-docking study (cognate docking) was undertaken in which each ligand of Figure 1 was separated from each receptor (protein) and each pair was then subjected to a rigid docking calculation. We first used the standard parameters given for each program.

5The results obtained with the five programs mentioned in the introduction are shown in Table 1. All the methods gave very good results. In fact, the root mean square deviation values were quite low (mean RMSD ≤ 1.0 Å); the standard errors of these mean RMSDs (SEM) are also indicated. The standard deviations (SD) are shown in the penultimate row and the median values in the last row. There are other appropriate statistical methods to evaluate the self-docking performance, but the mean and median RMSD values are sufficiently representative.

6The lowest mean RMSD values were those provided by AD Vina (0.55) and DOCK (0.57). The speed of AD Vina was worth noting (just a few seconds per calculation), whereas the other programs required ca. 5 min per calculation. Figure 2 shows graphically which complexes were better or more poorly reproduced in silico.

7When only the macrolide subset (1, 5, 6, 8, 9, 11, 13, and 19) was taken into account, the average RMSD values turned out slightly better than the averages for all 20 compounds shown in Table 1:

8The database of conformers for each ligand was obtained by means of the MD/LLMod method. The values of the RMSD between the crystal and the predicted structures are widely used to confirm whether a close-match docked pose was predicted or not by the docking simulation. An RMSD value ≤ 2 Å is fairly good.

9Conformers lying 0–42 kJ·mol−1 (ca. 0–10 kcal/mol) above the lowest-energy minimum were retrieved. For some structures (e.g., for 17/3EKS) a set of as few as 17 conformers was obtained, whereas for larger macrocycles >250 conformations had to be calculated (see Table 2). For the cases 15/3DV1 and 16/3DV5, with more rotatable bonds, the numbers of conformers were 1106 and 1402, respectively.

10The RMSD values obtained for the lowest-energy poses predicted for each ligand are shown in Table 2 and the corresponding graphic representation is shown in Figure 3. The mean RMSD values for Glide, AD Vina, and DOCK were as low as 1.29/1.31/1.34 Å, respectively. The mean RMSD value for AD 4.2 turned out to be slightly higher (1.69 Å), while for AD 3.0 was 2.62 Å.

11From a statistics viewpoint, we have also included at the bottom of Table 2 the mean RMSDs after ruling out the 2–4 RMSD values drawn in red (outliers) in the corresponding columns. It is clear that all the methods, except AD 3.0, afford very similar results. The RMSD values for each conformer and a Table with statistical results are given as Supplementary Material.

12As mentioned above, the five docking programs were initially examined using their standard parameters. However, “flexible rings docking” with the program AD 4.2 (“AD 4.2 Flex”), which allows for the flexibility of the macrocycle, was also carried out. Since poorer values were obtained than those under the standard conditions, this option was not considered further. No induced fit approaches, taking into account the possible flexibility of the side chains of the receptors, were checked for any program, since we did not perform virtual screenings, but instead evaluated the capacity of the different docking programs to reproduce a set of complexes that contain poly-functional ligands. The ligands were relatively large, compared to most “small molecules” (drugs), and they contained a relatively rigid core (the macrocycle).

13The case of 3/1NM6 requires some comments. All the methods, with the exception of AD 4.2 and Glide, gave RMSD values that were too high. A look at the docking images indicated that there are similar hydrophobic moieties (aromatic rings) at both extremes of 3 and some docking programs indistinctly located the folded conformers in either arrangement; they could not distinguish one pose from the other.

14The fact is that, for each program, there were 2–4 complexes in which the ligands had too high RMSD values. If these two “worst points” were removed, the mean RMSD values decreased significantly. For instance, the mean RMSD values for AD Vina, DOCK, and Glide decreased to 0.90, 0.94, and 0.95, respectively. In other words, 80%–90% of the structures seemed to be well reproduced.

15As far as the calculation speed is concerned, some representative computing times with fourth-generation i7 processors are indicated. With 5/1PKF (412 conformers) the order was:

16In the case of 16/3DV5 (1402 conformers), the order was:

17Thus, AD Vina is not only among the three most accurate methods but it is quite rapid as well. AD 4.2 appears in an intermediate position (good accuracy but slow), while AD 3.0 offered the worst results while requiring more computer time than any other (the least efficient). This is understandable since AD 3.0 is an old version of AD 4.2. We needed to evaluate the difference, since, as commented in the introduction, we had often used the program more than a decade ago.

18DOCK turned out to be as accurate as AD Vina, as happened in the self-docking procedure, but it was slower than AD Vina. For molecules with a very large number of conformers, the time increased exponentially. The performance of Glide, taking into account that the overall time invested was usually very low, was excellent. However, for molecules with a very large number of conformers, the time grew exponentially, as happened with DOCK. All the calculations with Glide were carried out with the standard precision (SP) method, but see the next section for the use of Glide XP.

19When only the macrolide subset (1, 5, 6, 8, 9, 11, 13, and 19) was analyzed, the average RMSD values were slightly better than the averages for all 20 compounds shown in Table 2:

20With the macrolide subset, we estimated how much the new algorithm contained in PSOVina reduced the calculation time. As shown in Table 3, in six out of the eight instances there were time saving, operating with the same computer, equipped with an i7-4790 processor.

21The RMSD values were similar, except for the complex of ligand 11/2IYA, for which it was higher with PSOVina (0.96) than with AD Vina (0.62). Owing only to this difference, the average RMSD value for PSOVina turned out to be slightly higher. Notwithstanding, the predicted binding affinities were practically identical.

22In several cases, the lowest-energy conformer predicted for each ligand was not the best pose (lowest RMSD). There was another conformer with a lower RMSD, the binding affinity of which was only around 0.0–0.5 kcal/mol above that of the lowest energy. In other words, there was often a conformer that overlapped with the crystalline structure better—we might call them the active conformer or that closest to the bioactive form, assuming that in many cases, but not always, the conformation of the ligand in the crystal structure may be the active conformation—while the score was only slightly different from that reported in Table 2. A summary of the results is shown in Table 4. At the bottom of Table 4, the average RMSDs together with the standard deviations are indicated for each method. Thus, for AD Vina the RMSD values are 0.83 ± 0.39 (they fall between 0.44 and 1.22), for AD 4.2 between 0.63 and 1.31, for DOCK between 0.45 and 1.41, and for Glide between 0.38 and 1.46. The medians are also included. It is observed that the mean and median RMSD values obtained with AD 3.0 are much higher than the rest. In short, all the methods (except AD 3.0) give nearly similar outcomes.

23Table 4 also shows the rank or pose order of the conformer with the optimal binding affinity; the number of red numbers in the pose-order columns is diagnostic (AD 3.0 is the worst in this sense).

24A comparison of Table 2 and Table 4 indicates that the performance of all the methods significantly improved, with the exception of two cases: (a) for 12/2VWC, Glide predicted an unacceptable RMSD value of 5.36; (b) for 14/2XX5, AD 3.0 gave an RMSD value of 4.81.

25Again, when only the macrolide subset (1, 5, 6, 8, 9, 11, 13, and 19) was considered, the average RMSD values changed slightly in relation to the general averages shown in Table 4:

26Re-docking of the main poses with computational methods that provide binding affinities closer to experimental values has a long tradition. Among the successful re-scorings reported, we will mention a few that have been published very recently. We realized that although AD Vina is a good compromise between speed and accuracy, we observed that the gaps between the predicted binding affinities for the best pose and for the next poses are minimal. In fact, many poses had very close scoring functions (often with differences of only 0.1 kal/mol).

27Re-scoring of the 20 molecules with X-Score software did not increase the small gaps. Re-scoring the AD Vina poses with AD 4.2, which obviously has advantages regarding the data compatibility between programs, improved the gaps between the poses, as shown in Figure 4 for the case of 3/1NM6. It should be noted that many conformations lying between −10.9 and −10.5 kcal/mol, only a few of which are shown in Figure 4 for the sake of simplicity, gave rise to larger energy differences after re-scoring: from −10.84 to −8.61 kcal/mol, in the examples shown. The conformer depicted in green (that of the lowest energy, −10.9 kcal/mol, according to AD Vina, see Table 2) moved up to −9.18 kcal/mol. The conformer depicted in magenta, which was not the lowest-energy one according to AD Vina (but it was the best pose), became the most favorable after the re-scoring.

28We call this re-scoring AutoDock 4.2//AutoDock Vina, or just AD 4.2//Vina, as done in quantum chemistry from the very beginning for single-point calculations of the energy at a high level of theory on equilibrium geometries obtained at a lower level. Table 5 summarizes the results for the AD 4.2//Vina re-scoring of the full series.

29Comparison with preceding Tables indicates how reliable this approach is. The mean RMSD and standard deviation does not significantly change (from 0.97 ± 0.34 to 1.05 ± 0.39), whereas the pose order of the lowest-energy predicted conformer moves to one (with only one exception, see 14). In other words, with AD 4.2//Vina, all the best poses for each ligand except one also appeared to be the lowest-energy conformers, in sharp contrast with the results in Table 4 for AD 4.2 (repeated in Table 5 for the sake of a direct comparison) and AD Vina. This may seem a fortuitous fact, to be interpreted by the authors of AutoDock, but it was well known that the scoring functions of AD 4.2 are more reliable than those of AD Vina, while AD Vina runs more rapidly than AD 4.2. Although the average RMSD value increased slightly (see Table 5), it is lower than the average values of 1.31 and 1.69 in Table 2 for AD Vina and AD 4.2, respectively.

30Taking the 20 poses of lowest energy in each case, we also tested the re-scoring of the four poorest DOCK results in Table 4 with GBSA score, which uses the pairwise GB solvation model with special parameters, and with Amber Score, which includes an implicit solvation and a short MD simulation that takes into account some flexibility of the receptor. Whereas no improvement was observed for 3/1NM6, poses were found for complexes 7/1UU3, 15/3DV1, and 19/3UYK which showed the optimal RMSD values and optimal relative energies, especially with the Amber Score (see Table 6).

31We also performed Glide XP calculations (eXtra Precision, as a matter of fact a docking procedure with more precise scoring functions) and a re-scoring with the popular PRIME MM-GBSA module for the six poorest results in Table 4, which were clearly improved with both programs in all cases (Table 6).

32Finally, we examined whether these poses, which are in principle more reliable, overlap with the experimental pose (crystalline structure) or not. The drawings and RMSD for the superimpositions (RMSD overlap) are given as Supplementary Material, but the key points are: With AD 4.2//Vina, 9 out of the 20 ligands overlap almost perfectly with their respective poses in the crystals (ligands 1, 3, 7, 8, 9, 10, 12, 17, and 18, RMSD overlap = 0.32–0.67 Å), 6 ligands can be superimposed moderately well with their crystalline structures, at least their cyclic moieties (RMSD overlap = 0.79–0.93), and there are four calculated structures that do not match sufficiently well (5, 11, 13, and 14, RMSD overlap = 1.04–1.10). The worst result is for ligand 19 (RMSD overlap = 2.14); we have no explanation for it, but the crystalline structure of the complex was poorly reproduced by most methods (see Table 2, Table 3, Table 4, Table 5 and Table 6).Three poses that improved after re-scoring of DOCK with Amber (Table 6) showed values of RMSD overlap for 7, 15, and 19 equal to 0.87, 1.12, and 0.75, respectively (two match partially and one not sufficiently well, also by visual inspection).For the six Glide poses re-scored with MM-GBSA (Table 6), the RMSD-overlap values were 1.13, 1.36, 0.69, 0.64, 0.73, and 0.67. In other words, three ligands (12, 13, and 19) were almost superimposable with their respective experimental structures, but there were two that did not match well. It is worth noting that 19/3UYK is better described by Glide (docking and re-scoring) than by AD 4.2 (docking and AD 4.2//Vina re-scoring, respectively).

33With AD 4.2//Vina, 9 out of the 20 ligands overlap almost perfectly with their respective poses in the crystals (ligands 1, 3, 7, 8, 9, 10, 12, 17, and 18, RMSD overlap = 0.32–0.67 Å), 6 ligands can be superimposed moderately well with their crystalline structures, at least their cyclic moieties (RMSD overlap = 0.79–0.93), and there are four calculated structures that do not match sufficiently well (5, 11, 13, and 14, RMSD overlap = 1.04–1.10). The worst result is for ligand 19 (RMSD overlap = 2.14); we have no explanation for it, but the crystalline structure of the complex was poorly reproduced by most methods (see Table 2, Table 3, Table 4, Table 5 and Table 6).

34Three poses that improved after re-scoring of DOCK with Amber (Table 6) showed values of RMSD overlap for 7, 15, and 19 equal to 0.87, 1.12, and 0.75, respectively (two match partially and one not sufficiently well, also by visual inspection).

35For the six Glide poses re-scored with MM-GBSA (Table 6), the RMSD-overlap values were 1.13, 1.36, 0.69, 0.64, 0.73, and 0.67. In other words, three ligands (12, 13, and 19) were almost superimposable with their respective experimental structures, but there were two that did not match well. It is worth noting that 19/3UYK is better described by Glide (docking and re-scoring) than by AD 4.2 (docking and AD 4.2//Vina re-scoring, respectively).

36The proteins selected for docking analysis satisfied the following criteria: they do not form covalent bonds with their respective ligands; the main rings have ≥11 members; and the structures were refined at resolutions <2.5 Å (Figure 1).

37Ligand coordinate files were extracted from the corresponding PDB files and used as reference structures for root mean square deviation (RMSD) calculations. The ligand input files were represented from their one-dimensional notations, SMILEs, with the aim of avoiding any bias from the bioactive conformation in conformer calculations. The 3D structures were generated with Corina and were refined with LigPrep.

38Macrocycle from the MacroModel 10.0 package was used for the conformational searches. OPLS2005 was the force field selected. The solvent effect was described by means of the GB/SA continuum solvent model. From the different methods of conformational search, we chose molecular dynamics/large-scale low-mode sampling (MD/LLMod). This MD/LLMod method is based on a two-stage search strategy. The first stage is high temperature MD-based simulated annealing, while the second stage is a large-scale low-mode search (LLMod). The MD/LLMod parameters were obtained after 5000 simulated annealing cycles and then 5000 energy-minimization iterations, with a maximum of 42 kJ·mol−1 (ca. 10 kcal/mol) above the lowest-energy conformer.

39The receptor input files were prepared as usual. First, we removed all water molecules, ligand atoms, and those ions that do not belong to the active site of the receptor from the corresponding PDB file. Next, hydrogen atoms were added, as well as protons, and amino acid side-chain partial charges (in accordance with the protonation state at the physiological pH) were assigned according to the requirements of each program. This was followed by a local minimization to relieve potential bad contacts. The minimization was performed in the presence of restraints to maintain the protein conformation very close to that observed in the experimental model. For those proteins with one cytochrome molecule close to the binding site (1PKF, 2C6H, and 2XBK), the Chimera USCF program was used; the corresponding hydrogen atoms and the AM1-Bcc partial charges of each atom were added and the cytochrome molecule was then included.

40AutoDock programs use a stochastic search algorithm (Lamarckian Genetic Algorithm, LGA, for AutoDock 3.0 and 4.2.6, and Iterated Local Search, ILS, for AutoDock Vina and PSO Vina). DOCK and Glide use a systematic search algorithm.

41Docking results were assessed with RMSD (root-mean-square deviation of the heavy atoms) of each predicted pose versus the crystal structure. The RMSD values were calculated with the rmsd.py script downloaded from the Script Center of the Schrödinger Resources website.

42We evaluated the performance of some popular docking programs at reproducing X-ray crystallographic structures that contain macrolides and related poly-functionalized ligands. Thus, the conclusions that follow regarding the performance of each method refer only and exclusively to these macrocyclic substrates when bound to rigid receptors; any generalization to other types of complexes will or would require an independent study. To our initial surprise, the outcome was very satisfactory, as most of the programs arrived at optimal poses in spite of the complexity of the structures. In general, for our set of 20 complexes, Glide 6.6 (Glide), AutoDock Vina (AD Vina), and DOCK 6.5 (DOCK) were efficient docking procedures for macrolides and analogues of intermediate size (Table 2). For large structures with many conformational isomers, the computing times of the programs Glide and especially of DOCK grew exponentially. If only true macrolides are considered (eight examples available), the preceding sentences are still valid. AD 4.2 outperforms AD 3.0, as expected for an updated and improved version, but we needed to confirm that this was the case for macrolide-like ligands. Consequently, our results obtained with AD 3.0 more than 10 years ago should be revised and corrected. When the best poses were compared (Table 4), AD Vina, Glide, DOCK, and AD 4.2 gave mean RMSD values below 1.0 Å. Since the energies given by AD Vina for the different conformers are very close, re-scoring with AD 4.2 (here so-called AD 4.2//Vina) was performed. This increased the energy gaps between the different poses and the agreement between the computational results and crystalline structures significantly improved. In general, re-scorings of DOCK and Glide also improved the agreement, which could have been expected but it had to be demonstrated for macrolide-like compounds. We hope to extend this study to other docking programs and to take advantage of the present findings in connection with the design and synthesis of analogs of anti-tumor macrolides.