Comparison of Molecular Mechanics, Semi-Empirical Quantum

Comparison of Molecular Mechanics, Semi-Empirical Quantum Mechanical, and Density Functional Theory Methods for Scoring Protein–Ligand Interactions...
1 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Comparison of Molecular Mechanics, Semi-Empirical Quantum Mechanical, and Density Functional Theory Methods for Scoring Protein−Ligand Interactions Nusret Duygu Yilmazer and Martin Korth* Institute for Theoretical Chemistry, Ulm University, Albert-Einstein-Allee 11, 89069 Ulm, Germany ABSTRACT: Correctly ranking protein−ligand interactions with respect to overall free energy of binding is a grand challenge for virtual drug design. Here we compare the performance of various quantum chemical approaches for tackling this so-called “scoring” problem. Relying on systematically generated benchmark sets of large protein/ ligand model complexes based on the PDBbind database, we show that the performance depends first of all on the general level of theory. Comparing classical molecular mechanics (MM), semiempirical quantum mechanical (SQM), and density functional theory (DFT) based methods, we find that enhanced SQM approaches perform very similar to DFT methods and substantially different from MM potentials.

1. INTRODUCTION Accurate but fast computational modeling of complex systems is in the focus of researchers from several scientific fields. Computational approaches for the treatment of biomolecular systems are of especially high interest because of their importance for in silico drug design.1 A central challenge for such methods is the so-called “scoring” problem: Virtual drug design involves not only the “docking” of ligand (drug candidate) molecules into protein receptors but also their correct ranking in terms of overall free energy of binding.2 The overall goal is then that only the most promising candidates have to be tested in subsequent experiments, so that virtual screening can help to greatly reduce the experimental effort for developing new drugs, an approach that has now become a standard procedure for many pharmaceutical companies.3 Docking requires the prediction of the preferred orientation of the ligand to the receptor upon binding, but it is generally agreed on that common docking software can provide a reasonably good accuracy for identifying the most important binding poses.4,5 Scoring on the other hand is still a major problem because the accurate ranking of binding affinities turns out to be a very difficult task for the commonly used “scoring functions”.6,7 Such scoring functions treat the interactions of ligands and proteins at a rather low theoretical level, so that increasing the theoretical and computational effort seems to be a promising way to improve on the scoring problem; however, systematic studies of the benefits of higher-level methods for describing protein−ligand interactions are still very rare,8−10 though a number of studies on smaller host−guest model systems do exist.11,12 Fortunately, several research groups have made sets of protein−ligand complex structures and experimentally determined binding affinities, like for instance the PDBbind database,13,14 publicly available, so that the further © 2013 American Chemical Society

development and evaluation of scoring methods is possible on the ground of a large collection of experimental data.15 Important opportunities for the improvement of scoring functions include the more accurate treatment of polarization, solvation, and entropic effects, with the latter as probably the most pressing issue.4,5 A sometimes utter disregard of the first two problems is nevertheless quite premature because a really accurate assessment of entropic effects is only possible on the ground of correct enthalpies. Semiempirical quantum chemical (SQM) methods are very promising candidates for improving upon classical potentials with a more accurate quantummechanical treatment of polarization and solvation effects.16 Accordingly, these methods are now investigated in several research groups,8,17 with a large body of very successful studies especially by Hobza and co-workers,18 and were already successfully applied to systems up to tens of thousands of atoms using “linear scaling” approaches16as is of course also the case for linear scaling DFT schemes.19 We have recently contributed to this field with the development of empirical corrections for noncovalent interactions, which are of high importance for protein−ligand systems but usually treated quite badly at the semiempirical level.20 SQM methods augmented with our corrections reach the accuracy of dispersion-corrected density functional theory (DFT-D) approaches for a large number cases investigated, while still being about 3 orders of magnitude faster.21,22 In the following we would like to assess the performance of SQM and DFT electronic structure methods for the scoring of protein/ligand interactions. Our main focus is the comparison of SQM with DFT approaches Received: March 19, 2013 Revised: May 2, 2013 Published: June 12, 2013 8075

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

sophisticated Amber ff99sb/GAFF approach for further comparison.

using realistic protein−ligand model structures, but in the course we also benchmark several proposed improvements and technical parameters at the respective SQM and DFT levels and make a first assessment of the differences to classical molecular mechanics (MM) results. Our choice of methods was motivated by the results of investigations into the performance of computational methods for smaller biomolecular model systems.23 From these studies, dispersion-corrected DFT and enhanced SQM methods (as well as self-consistent charge density functional tight binding24) emerged as the most promising candidates for going beyond the standard classical treatment of larger biomolecules. Errors with respect to highlevel (“gold-standard”) CCSD(T)/CBS data are on average below 1 kcal/mol for all of these methods, thereby improving greatly upon the accuracy found for state-of-the-art force field approaches.25 As DFT functionals, we chose BP86, for which the sophisticated COSMO-RS solvation model was originally parametrized (though we found only a very small dependency on the functional), PBE, the most commonly used GGA functional, and TPSS, which is among the most accurate functionals for biomolecular interactions when corrected for dispersion interactions (although differences between functionals are actually small for biomolecular interactions). As SQM methods we chose AM1, an older approach still widely used in the context of biomolecular interactions and PM6, which was found to be the most accurate approach apart from Thiel’s OMx methods,26 which are unfortunately not yet available for atoms beyond the first row. Motivated by our own contributions to the development of enhanced SQM methods, we are also interested in the importance of different dispersion (“-D”) and hydrogen-bond (“-H”) correction schemes. The dispersion term is usually set up similar to DFT-D ones and seems therefore less questionable, but the hydrogen-bond term clearly merits further investigation, which is why we have included two different variants of this term and an add-on term for halogen bonds in our tests: The DH2 correction uses socalled second-generation terms, which include all geometric information about the spatial arrangement of the hydrogen bond but are based on the distinction between donor and acceptor atom types, through the use of hydrogen-acceptor distances, which are not symmetrical with respect to the electronegative atoms involved. In the case of SQM methods, this dependence on the distance between the hydrogen and the acceptor atom is a major drawback, as it makes themunlike the common empirical dispersion correctionsan implicit bond-type term: If the acceptor atom changes, for instance in the case of proton transfer from the donor to the acceptor, the correction breaks down. The DH+ correction is a thirdgeneration approach which models hydrogen bonds as an interaction between two electronegative atoms (X, Y) that is smoothly switched on by the favorable placement of one (or more) hydrogen (H) atom(s) in between them. This avoids the mentioned problem for the application to SQM methods, connected to the necessity of a bond definition for the donor hydrogen in the second-generation approaches. The DH+ correction is thus more robust and transferable than DH2 but was able to achieve the same accuracy as DH2 with just one fit parameter per electronegative element involved in hydrogen bonding in benchmark set calculations. For the PM6-DH2X method, an additional H2-type correction term is applied also to halogen bonds.27 Concerning MM models we initially chose MMFF94, which we have found to be among the most accurate for biomolecular interactions, but then included also the more

2. GENERATION OF MODEL SYSTEM BENCHMARK SETS To acquire preferably realistic systems for the evaluation of QM and MM methods, we have set up benchmark sets of protein− ligand structures based on the PDBbind database and prepared them in a systematic way: Starting from the 1300 protein− ligand complexes in the refined set of the PDBbind2007 database13,14 cuts were made around the docked ligand at increasingly larger distances, and all systems with structural errors or uncertainties within the respective cutoff distance were discarded. Also, cases with ligand files with irregularities were directly skipped from the following. Problems which were encountered include atoms with an unreasonable number of bonds (assigned via comparison to sums of van der Waals radii), atoms missing from amino acids, and atoms too close to each other (judged by van der Waals radii). We did so to avoid uncertainties from fixing problems in the original PDBbind files. When generating our model systems, residues with at least one atom within the respective cutoff distance are kept completely, and resulting terminal structures are capped with hydrogen atoms. The resulting cut is therefore systematically larger (by about 3−5 Å) than the chosen cutoff distance. Protonation was necessary for histidine residues, which were all assumed to be neutral, and water molecules within the pocket were discarded. These are harsh approximations, which make comparison with experiment hardly possible at this stage but are valid for our purpose of comparing different computational methods. Overall charges were assigned according to automatic Lewis structure analysis (i.e., assigning charges to atoms according to the number of missing bonds found via comparison to sums of van der Waals radii) and doublechecked with the assignment by MOPAC. Cases in which our analysis and MOPAC did not agree were again skipped to avoid uncertainties from charge assignment. The PDBbind2007 refined set complexes range from about 1000 to almost 90 000 (on average around 7000) atoms in size, with ligands from about 20 to 140 (on average 50) atoms. The size of our model complexes is dependent on the applied cutoff distance, resulting in structures from about 80 to 820 atoms for 3.0 Å, up to about 500 to 2900 atoms for 10.0 Å. The 10.0 Å benchmark set was used to analyze the convergence of the ranking with respect to the cutoff distance, which was possible at the SQM level only. The performance of different methods, as well as technical parameters and possible improvements, was assessed with the 3.0 and 5.0 Å benchmark sets. 3. COMPUTATIONAL DETAILS MMFF9428 force field calculations were done with Open Babel.29 Amber ff99sb30 and GAFF31 force field calculations with Gasteiger charges were done with Amber 11,32 utilizing the AmberTools for automatic input preparation. Completely automatic preparation was possible for 352 entries, and all other cases were skipped to avoid uncertainties from fixing problems with manual atom type assignment. Semiempirical AM1,33 AM1-D,21,34 PM6,16 PM6-D,21 PM6-DH2,21 PM6DH2X,27 and PM6-DH+22 calculations were done with MOPAC2012,35 making use of the MOZYME linear-scaling algorithm and the COSMO36 solvation model. At the SQM and MMFF94 level, calculations on all generated model systems 8076

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

comparison to our DFT reference values. We will accordingly denote an R value of 1.00 in combination with a τ value of 0.94 or higher as perfect agreement with the reference. Figure 1 shows PM6-DH+ binding energies for all four benchmark sets and allows us to investigate the convergence of binding energies with increasing cutoff distance. As a result, a reasonable cutoff distance can be chosen for further comparison. Figure 2 presents values for different SQM methods, consisting of a basic model and a number of different correction schemes for noncovalent interactions. Evaluation of these data allows us to pick the best suited SQM methods. Figure 3 shows data for different DFT functionals and basis sets and is complemented by Figure 4, which illustrates the effect of going from the simple COSMO to the more advanced COSMO-RS solvation model, as well as Figure 5, which presents similar data for different dispersion correction schemes. Figure 6 finally gives a comparison of the best suited MM, SQM, and DFT methods before we take a closer look at the problems of SQM methods in Figure 7.Tables 2 and 3 list additional data for different cut sizes and different MM methods, and Table 4 summarizes the error statistics discussed in paragraph E (see below for details). A. Cutoff Distances. Figure 1 shows that the ranking of the calculated binding energies converges quite quickly with increasing cutoff distance from R = 0.83 at 3.0 Å over R = 0.96 at 5.0 Å to R = 1.00 at 7.0 Å. Although there are still sizable contributions beyond a distance of 7−10 Å, these contributions seem to be rather uniform and contribute similarly to all cases. This is reasonable because even the smallest protein is much larger than the largest ligand, and all proteins therefore offer a rather uniformly large number of interaction partners. Additionally, the almost perfect convergence at about 7.0 Å (corresponding to an effective distance of about 10−12 Å) is in good agreement with what would be expected from studies on small model systems.18 Figure 1 therefore indicates that quantum mechanical calculations can be restricted to smaller model systems without losing predictive capability. Looking at Kendall τ values in Table 1, we find the same overall trend as for the Pearson R values. Nevertheless, rank correlation between the 7.0 and 10.0 sets is less perfectly converged than suggested by the linear correlation measure. Nevertheless, even for the 7.0 set, the differences are similar to those between different DFT functionals. Though it would have been desirable to continue with the 10.0 Å set for all of the following, we had to restrict ourselves to the smaller sets for the DFT calculations, for which a very similar cutoff distance convergence behavior is found. Real scoring application would of course require the larger cutoff, which is fortunately well within reach of the SQM methods. B. SQM Methods. Using Figure 2 we are able to make some conclusions about which SQM methods are best suited for our purpose. Studies on small model systems suggest that PM6-DH+ should be among the most accurate models,21,22 and this view is indeed again supported by the comparison with DFT-D data (see below, paragraph D). We therefore use PM6DH+ as our SQM reference and look at the correlation of AM1, PM6-D, PM6-DH2, and PM6-DH2X data with our reference. From this comparison, we can deduce that there seems to be an advantage of PM6 over AM1, a significant importance of dispersion corrections, but only a negligible contribution from other empirical corrections for hydrogen and halogen bonding. Looking at Kendall τ values in Table 1 we find the same trends, with an increased emphasis on the differences between methods. This is in good agreement with recent findings by

were possible, giving 736 data points at 3.0 Å, 733 at 5.0 Å, 725 at 7.0 Å, and 714 at 10.0 Å cutoff distance. In some cases, enlarging the cutoff distance resulted in trapping structural errors, so that additional cases had to be discarded for the larger cutoff distances. BP86,37,38 PBE,39 and TPSS40 DFT calculations with empirical dispersion corrections of D2,41 D3,42 and D3 plus three-body-dispersion (named D33 in the following) type were done with the “huge” variant of Turbomole 6.443,44 using TZVP and TZVPP Gaussian AO basis sets,45 the RI approximation for two-electron integrals,46,47 and COSMO as well as COSMO-RS (via COSMOtherm) solvation models.36 At the DFT level, not all generated model systems could be successfully treated because of software and RAM limitations as well as SCF convergence problems. As a result, we finally arrived at 695 BP86/TZVP and 487 BP86/TZVPP data points for the 3.0 Å, as well as 539 BP86/TZVP, 539 PBE/TZVP, and 513 TPSS/TZVP data points at 5.0 Å. Each of these data points corresponds to an interaction energy computed from singlepoint energy calculations on the complex, pocket, and ligand structures in the unrelaxed PDBbind geometry E interaction = Ecomplex − Epocket − E ligand

4. RESULTS AND DISCUSSION We have calculated binding energies for four benchmark sets, named after the used cutoff distance 3.0, 5.0, 7.0, and 10.0 Å, with about 700 systems each, making use of MM, SQM, and DFT methods and several possible enhancements. We aim at a better ranking of protein−ligand interactions and accordingly focus on the correlation of the data points, which is presented in graphical form for easy assessment and comparison. A common numerical measure for correlation is the Pearson R value, which assesses the linear correlation of two sets of data. For ranking purposes, it is more instructive to look at nonlinear correlation measures, like Spearman’s ρ or Kendall’s τ, which are a better measure of how well the actual rankings are correlated. We therefore give Kendall τ values in comparison to Pearson R values in Table 1 for further analysis. For the comparison of different DFT-D functionals we find R values of 1.00 and τ values of 0.94 (see Table 1), and these values are thus limits for what we can claim to be an improvement in Table 1. Kendall τ Values in Comparison to Pearson R Values for All Data Presented in Figures 1 to 7 entry

figure

R

τ

PM6-DH+ 3.0 vs 10.0 PM6-DH+ 5.0 vs 10.0 PM6-DH+ 7.0 vs 10.0 AM1 vs PM6-DH+ PM6 vs PM6-DH+t PM6-D vs PM6-DH+ PM6-DH2 vs PM6-DH+ PM6-DH2X vs PM6-DH+ PBE-D/TZVP vs BP86-D/TZVP TPSS-D/TZVP vs BP86-D/TZVP BP86: TZVPP vs TZVP BP86: COSMO vs COMSO-RS BP86: D2 vs D3 BP86: D33 vs D3 PM6-DH+ vs BP86-D/TZVP MMFF94 vs BP86-D/TZVP

1a 1b 1c 2a 2b 2c 2d 2e 3a 3b 3c 4 5a 5b 6a 6b

0.83 0.96 1.00 0.80 0.89 0.99 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.93

0.76 0.90 0.96 0.56 0.67 0.90 0.94 0.94 0.94 0.94 0.99 0.95 0.94 0.99 0.93 0.76 8077

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

Figure 1. Correlation between PM6-DH+ data for benchmark sets generated with different cutoff distances. From left to right: (a) 3.0, (b) 5.0, and (c) 7.0 in comparison to 10.0 Å values. All computations with COSMO solvation (COSMO label skipped from axis labeling for better readability).

Figure 2. Correlation between different SQM approaches. From top left to bottom right: (a) AM1, (b) PM6, (c) PM6-D, (d) PM6-DH2, and (e) PM6-DH2X in comparison to PM6-DH+ values. All computations with COSMO solvation (COSMO label skipped from axis labeling for better readability).

Ryde an co-workers.8 Since the computational cost is very little in comparison to the PM6 calculation and since very substantial improvements from hydrogen-bond corrections were proven for small systems, we nevertheless recommend to use PM6DH2 or PM6-DH+ when investigating biomolecular systems with SQM methods. The reason is the following: while hydrogen-bond corrections might not play a decisive role for protein−ligand interaction energies, they are very likely to do so for protein geometries. Geometry optimization of hydrogenbonded systems with SQM methods was found to be much improved with hydrogen-bond corrections18 because without

them the delicate balance between the uniform dispersion and the directional hydrogen-bond interactions is not correctly modeled, so that hydrogen-bond binding motifs are lost. C. DFT Methods. Figure 3 indicates that the choice of the DFT functional is of little importance when looking at rankings only. Similarly, increasing the basis set beyond the TZVP level seems to be of little use here, as the observed gain in interaction energy is rather uniform for all investigated cases. Another allclear signal is given by Figure 4, which shows that the differences between the simple COSMO solvent model and the more sophisticated COSMO-RS approach (including several 8078

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

Figure 3. Correlation between different DFT functionals and basis sets. From left to right: (a) PBE-D2/TZVP, (b) TPSS-D2/TZVP, and (c) BP86D2/TZVPP in comparison to BP86-D2/TZVP values. All computations with COSMO solvation (COSMO label skipped from axis labeling for better readability).

additional terms, for instance, for hydrogen-bond interactions with the solvent) do not matter too much for our purpose. Also different dispersion correction schemes do not seem to be a decisive component when just ranking protein−ligand interactions as shown in Figure 5. Neither the latest (D3) improvement of the widely used Grimme scheme (D2) nor the additional inclusion of three-body dispersion effects (termed D33 here) has a large impact on the ranking. Also the Kendall τ values in Table 1 support these conclusions. This is of course no counter-argument for recent claims about the importance of COSMO-RS and D3/D33 in the literature when calculating absolute values for noncovalently bound systems for direct comparison with experiment,12 only that the case seems to be different when looking at rankings on a large scale only (but see below for an analysis at a smaller scale). D. Comparison of Methods. In Figure 6 we compare methods at the MM, SQM, and DFT level with each other. Judging from the DFT values, SQM methods seem to be a clear improvement of MM methods. Looking at the Kendall τ values in Table 1 we find the same overall trend, with an even stronger emphasis on the differences especially between SQM and MM methods. (We are presenting MMFF94 data here, but it was

Figure 4. Correlation between the solvation contributions from COSMO and COSMO-RS for BP86/TZVP calculations.

Figure 5. Correlation between different dispersion correction schemes. From left to right: (a) D3 against D2 and (b) D33 against D3 contributions for BP86 calculations. 8079

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

Figure 6. Correlation between methods at different theoretical levels. From left to right: (a) MMFF94 and (b) PM6-DH+ each in comparison to BP86-D2/TZVP values.

Figure 7. More detailed comparison between PM6-DH+ and BP86-D2/TZVP. From top left to bottom right: (a) Overall interaction energy, (b) energy component scaled over all interaction energy (see text), (c) dispersion contribution, (d) COSMO solvation contribution, and (e) electronic contribution. All computations with COSMO solvation (COSMO label skipped from axis labeling for better readability).

have nevertheless evaluated the more sophisticated Amber ff99sb/GAFF approach for the subset of structures for which automatic atom-type assignment was possible; see below for the

recently shown that interaction energies from most MM methods are highly correlated,48 so that we do not expect qualitatively different results for other MM approaches. We 8080

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

Table 2. Pearson R and Kendall τ Values for All Data Presented in Figure 7 in Comparison to the Corresponding Values for the 3.0 Set of Structures 3.0

contributions to the respective DFT values, multiplied with the slope of the linear regression of the added up contributions to the overall DFT values. Table 2 shows that by making use of the same scaling factors a similar improvement is found also for the smaller 3.0 cuts. We have argued above that because MM interaction energies are highly correlated48 we do not expect qualitatively different results for different MM approaches. To check whether this is really true for more sophisticated approaches like the combined Amber ff99sb/GAFF approach, we have made additional comparisons for all cases where automatic atom type assignment was possible (roughly half of the cases; see above). The results, given in Table 3, clearly support our initial assumption. The performance of MMFF94 and ff99sb/ GAFF is found to be similar in comparison to DFT values judging from both Pearson R and Kendall τ values: R/τ are 0.93/0.76 for the former vs 0.96/0.73 for the latter method, while R/τ values are 1.00/0.93 for the SQM approach. A comparison based on energies including solvation effects is even more biased in favor of the QM methods because of the comparably large differences between QM and MM solvation models: R/τ are 0.60/0.51 for ff99sb/GAFF/GBSA in comparison to 0.92/0.74 for the SQM method. Using more sophisticated charge and solvation models would very likely improve these results, but given the very small impact of the rather big change from MMFF94 to ff99sb/GAFF, we feel justified in not expecting qualitatively different results. E. Accuracy on a Smaller Scale. Up to here we have talked about correlation values only and looked at them on a rather large scale. Ryde and co-workers have shown that absolute binding energies calculated with implicit solvent models are anyhow problematic,50 so that using scaled or shifted QM contributions seems unavoidable in any case, but apart from prescreening applications one might still be interested to look at the differences on a smaller scale. A good correlation for a set of very diverse cases with energy differences of up to 400 kcal/mol can still mean substantially different results when we compare a number of cases within a much smaller energy range of say 10−20 kcal/mol. This can be seen from looking at the error statistics corresponding to the given correlation values, given in Table 4. As we are interested in rankings first of all, we nevertheless want to ignore systematic shifts or scaling factors between values. We therefore give error measures also for data adjusted by scaling and shifting in the following. Scaling is done by multiplying all data points with one over the slope of the corresponding linear regression, shifting by subtracting minus one time the offset of the linear regression from all data points. We find very large mean deviations (MDs) and mean absolute deviations (MADs) between nonenhanced and enhanced SQM methods (−54 and 54 kcal/mol for AM1 vs

5.0

entry

figure

R

τ

R

τ

overall overall scaled solvation dispersion electronic

7a 7b 7c 7d 7e

0.92 0.95 1.00 0.99 0.93

0.74 0.79 0.95 0.92 0.81

0.93 0.95 0.99 0.98 0.99

0.77 0.79 0.95 0.89 0.84

results.) In Figure 7 we take a closer look at the DFT and SQM values. It is first of all noteworthy that the agreement between DFT and SQM gets worse (from R = 0.99 to 0.93, see Figure 7a), when solvation effects are included at the respective level of theory. Also the comparison between dispersion contributions indicates significant differences between the SQM and DFT approaches. These are important findings because they show in what respect SQM methods need to be further improved (and that there is quite some room for improvement): Better polarization is needed for better implicit solvent model results, and the balance of intermolecular interactions is quite badly reproduced at SQM in comparison to the DFT level, very likely also because common SQM methods (and especially PM6) already model some noncovalent effects (but without the proper long-range behavior), so that adding dispersion corrections is a somewhat dodgy process. Recent development of more physically sound SQM methods by for instance Truhlar and co-workers49 could therefore be of high importance for quantum chemistry based scoring functions, only that also the balance of intermolecular interactions should be taken into account right from the start. We find no correlation of the deviations with minimum or average interaction distances, also when excluding hydrogen bonds (for which rather small distances could still be fine) or hydrogen atoms in general (for which the effects should be less pronounced than for heavier atoms). The very high correlation of the different components of the SQM and DFT values (see Figure 7c−e, R = 0.99/0.98/0.99 for the dispersion/solvation/ electronic components) suggests a quick test: Scaling the SQM components before adding them up results in an improved correlation of the SQM with the DFT values (see Figure 7b, R = 0.95). Such an “energy component scaled” (ECS) PM6-DH+ variant, with scaling factors of 0.96, 1.7, and 0.91 for the electronic, dispersion, and solvation energy contributions, could thus even serve as a preliminary solution for the problem. E ECS = a ·E electronic + b·E D + c·E COSMO

The parameters a = 0.96, b = 1.7, and c = 0.91 were here derived from the slopes of the linear regressions for the SQM

Table 3. Comparison of Amber(ff99sb/GAFF), MMFF94, SQM, and DFT Data ff99sb/GAFF vs MMF94 ff99sb/GAFF vs BP86-D2/TZVP ff99sb/GAFF/GBSA vs BP86-D2/TZVP/COSMO MMFF94 vs BP86-D2/TZVP PM6-DH+ vs BP86-D2/TZVP PM6-HD+/COSMO vs BP86-D2/TZVP/COSMO a

R

τ

0.95 0.96 0.60 0.93 1.00 0.92(0.95a)

0.75 0.73 0.51 0.76 0.93 0.74(0.79a)

Energy component scaled variant; see text for details. 8081

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

Table 4. Mean Deviations (MDs) and Mean Absolute Deviations (MADs) for All Data Presented in Figures 2 to 6

a

entry

figure

MD

MAD

MDa

MADa

AM1 vs PM6-DH+ PM6 vs PM6-DH+ PM6-D vs PM6-DH+ PM6-DH2 vs PM6-DH+ PM6-DH2X vs PM6-DH+ PBE-D/TZVP vs BP86-D2/TZVP TPSS-D/TZVP vs BP86-D2/TZVP BP86: TZVPP vs TZVP BP86: COSMO vs COMSO-RS BP86: D2 vs D3 BP86: D33 vs D3 PM6-DH+ vs BP86-D2/TZVP MMFF94 vs BP86-D2/TZVP

2a 2b 2c 2d 2e 3a 3b 3c 4 5a 5b 6a 6b

−54 −41 −7 2 2 −2.6 −1.2 0.5 2.9 3.2 −1.2 −12 −55

54 41 7 2 2 3.0 2.4 0.8 6.8 3.4 1.2 14 57

−5 −2 0 0 0 0.0 0.0 0.0 −0.4 0.0 0.0 −0.1 0.4

20 14 4 2 2 2.3 2.3 0.4 4.4 1.9 0.4 7.1 30.9

Values scaled and shifted before error assessment; see text for details.

PM6-DH+, −41 and 41 kcal/mol for PM6 vs PM6-DH+) and still sizable differences when we correct for systematic shifts and scaling (−5 and 20 kcal/mol for AM1, −2 and 14 kcal/mol for PM6). MD and MAD values are much smaller in between enhanced SQM methods (−7 and 7 kcal/mol for PM6-D, 2 and 2 kcal/mol for PM6-DH2 and PM6-DH2X) and become very small when systematic shifts and scaling are taken into account (0 and 4 kcal/mol for PM6-D, 0 and 2 kcal/mol for PM6-DH2 and PM6-DH2X) MD/MAD values for PBE and TPSS versus BP86 with TZVP basis sets are −2.6/3.0 and −1.2/2.4 kcal/mol before shifting and scaling and 0.0/2.3 and 0.0/2.3 kcal/mol after shifting and scaling. MD/MAD values for BP86/TZVPP versus BP86/TZVP values are 0.5/0.8 kcal/mol before and 0.0/0.4 kcal/mol after adjustment, with the TZVPP results shifted by −0.7 kcal/mol and scaled by a factor of 0.97. MD/MAD values for COSMO versus COSMO-RS are 2.9/6.8 kcal/mol and change to −0.4/4.4 kcal/mol after adjustments, with the COMSO-RS results shifted by −8.4 kcal/mol and scaled by a factor of 1.05. MD/MAD values for D2 versus D3 and D33 versus D3 are 3.2/3.4 and −1.2/1.2 kcal/mol and change to 0.0/1.9 and 0.0/0.4 kcal/mol after adjustments with a shift of −1.2 and −0.1 kcal/mol from D3 and D33 and scaling factors of 1.03 and 0.98 from D3 and D33. Looking at the different theoretical levels, we find large MD/ MAD values of −12 and 14 kcal/mol for the comparison of PM6-DH+ versus BP86/TZVP but much larger MD/MAD values of −55 and 57 kcal/mol for the comparison of MMFF94 with BP86/TZVP. After shifting the SQM and MM values by 0.8 and 56 kcal/mol and scaling by 0.91 and 1.00, MD and MAD values drop to −0.1 and 7.1 kcal/mol for PM6-DH+ and 0.4 and 30.9 kcal/mol for MMFF94. The observed MADs of about 7 kcal/mol for the SQM and about 30 kcal/mol for the MMFF94 method illustrate the advantages of the former especially when looking at smaller energy scales and/or less diverse systems. Typical deviations from DFT values are therefore at the size of about 1% for SQM and 3% for MM methods on the electronic scale (with an energy range of about 1000 kcal/mol) and 5% for SQM and 15% for MM methods taking also solvation into account (with a resulting energy range of about 200 kcal/mol, excluding additional errors from the respective solvent models). For our ECS-PM6-DH+ variant, we arrive at MD/MAD values of 0.0/8.7 kcal/mol with respect to the DFT reference values after adjusting for a systematic shift of 18 kcal/mol

between DFT and SQM values, which is only a minor improvement over the unscaled version with MD/MAD values of −0.7/9.3 kcal/mol (now including the errors from the respective solvation model). F. Remaining Issues. Though we think that we have shown that enhanced SQM methods and model system approaches are a promising path to better scoring functions, it is perfectly clear to us that we are still quite far away from reaching the final goal. The problem is nicely illustrated by comparison of the DFT-D values with experimental data, where essentially no correlation at all is found. We believe the main reasons for this to be our simple preparation of the model systems (relying on a single binding mode, ignoring protonation problems, cavity water, and relaxation effects, etc.) and probably even more important the neglect of entropic effects, but would again like to emphasize that the assessment of entropic effects will profit greatly from working with accurate enthalpic contributions. Otherwise, besides actually accounting for entropy changes, entropic terms have to also “magically” correct enthalpic errors. To summarize our findings: •quantum mechanical calculations can be restricted to smaller model systems without losing predictive capability, and the resulting model systems are well within the reach of SQM methods but still challenging for DFT approaches •SQM methods perform significantly different, and based also on former studies we recommend using SQM-DH+ methods •when ranking very diverse cases over a comparably large energy scale, we find almost no differences between different dispersion corrected DFT approaches and implicit solvent models, but based again also on former studies we recommend using DFT-D3/TZVP/COSMO(-RS) methods •when comparing methods at different theoretical levels, we find that enhanced SQM approaches like SQM-DH+ perform very similar to DFT-D and improve substantially upon classical potentials in this respect •looking at differences on a smaller energy scale, the peculiarities of these approaches become more important, with typical deviations to the DFT values of 5% for SQM and 15% for MM methods after adjustment for systematic differences

5. CONCLUSIONS We have presented a comparison of different SQM and DFT approaches based on a large number of realistic protein−ligand model complexes. Comparably small differences are found for 8082

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

Known Three-Dimensional Structures. J. Med. Chem. 2004, 47, 2977− 2980. (14) Wang, R.; Fang, X.; Lu, Y.; Yang, C. -Y.; Wang, S. The PDBbind Database: Methodologies and Updates. J. Med. Chem. 2005, 48, 4111− 4119. (15) Nicola, G.; Liu, T.; Gilson, M. K. Public Domain Databases for Medicinal Chemistry. J. Med. Chem. 2012, 55, 6987−7002. (16) Stewart, J. P. Optimization of Parameters for Semiempirical Methods V: Modification of NDDO Approximations and Application to 70 Elements. J. Mol. Model. 2007, 13, 1173−1213. (17) Raha, K.; Merz, K. M. Large-Scale Validation of a Quantum Mechanics Based Scoring Function: Predicting the Binding Affinity and the Binding Mode of a Diverse Set of Protein−Ligand Complexes. J. Med. Chem. 2005, 48, 4558−4575. (18) Hobza, P. Calculations on Noncovalent Interactions and Databases of Benchmark Interaction Energies. Acc. Chem. Res. 2012, 45, 663−672. (19) Cole, D. J.; Skylaris, C. -K.; Rajendra, E.; Venkitaraman, A. R.; Payne, M. C. Protein−Protein Interactions from Linear-Scaling FirstPrinciples Quantum-Mechanical Calculations. Europhys. Lett. 2010, 91, 37004−37009. (20) Korth, M. Empirical Hydrogen-Bond Potential Functions-An Old Hat Reconditioned. Chem. Phys. Chem. 2011, 12, 3131−3142. (21) Korth, M.; Pitonak, M.; Rezac, J.; Hobza, P. A Transferable HBonding Correction for Semiempirical Quantum-Chemical Methods. J. Chem. Theory Comput. 2010, 6, 344−352. (22) Korth, M. Third-Generation Hydrogen-Bonding Corrections for Semiempirical QM Methods and Force Fields. J. Chem. Theory Comput. 2010, 6, 3808−3816. (23) Hobza, P.; Muller-Dethlefs, K. Non-covalent Interactions. Theory and Experiment; RSC Publishing: London, 2010. (24) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260−7268. (25) Paton, R. S.; Goodman, J. M. Hydrogen Bonding and ?-Stacking: How Reliable are Force Fields? A Critical Evaluation of Force Field Descriptions of Nonbonded Interactions. J. Chem. Inf. Model. 2009, 49, 944−955. (26) Weber, W.; Thiel, W. Orthogonalization Corrections for Semiempirical Methods. Theor. Chem. Acc. 2000, 103, 495−506. (27) Rezac, J.; Hobza, P. A Halogen-Bonding Correction for the Semiempirical PM6 Method. Chem. Phys. Lett. 2011, 506, 286−289. (28) Halgren, T. A. Merck Molecular Force Field. I. Basis, Form, Scope, Parameterization, And Performance of MMFF94. J. Comput. Chem. 1996, 17, 490−519. (29) O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchinson, G. R. Open Babel: An Open Chemical Toolbox. J. Chem. Inf. 2011, 3, 33−47. (30) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (31) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (32) Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossvazary, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M. -J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 11; University of California: San Francisco, 2010. (33) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. P. AM1: A New General Purpose Quantum Mechanical Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (34) Mc Namara, J. P.; Hillier, I. H. Semi-Empirical Molecular Orbital Methods Including Dispersion Corrections for the Accurate

different DFT methods when screening diverse structures over a wide energy range, as well as a very similar performance of SQM and DFT methods as opposed to MM potentials. Looking at a smaller scale we find typical deviations with respect to DFT-D values of about 5% for SQM and 15% for MM methods after adjusting for systematic differences. Further analysis has shown that there is still room for improvement on the SQM side, also because the balance of electronic, dispersion, and solvation interaction contributions is very different from the one found for the presumably more accurate DFT methods. Getting the delicate balance of intermolecular forces right should therefore be a major goal for further SQM method development.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support from the Barbara Mez-Starck Foundation is gratefully acknowledged. We gratefully thank the bwGRiD project for computational resources.

■

REFERENCES

(1) Jorgensen, W. L. Efficient Drug Lead Discovery and Optimization. Acc. Chem. Res. 2009, 42, 724−733. (2) Leach, A. R.; Shoichet, B. K.; Peishoff, C. E. Prediction of Protein−Ligand Interactions. Docking and Scoring: Successes and Gaps. J. Med. Chem. 2006, 49, 5851−5855. (3) Klebe, G. Virtual Ligand Screening: Strategies, Perspectives and Limitations. Drug Discovery Today 2006, 11, 580−594. (4) Gilson, M. K.; Zhou, H.-X. Calculation of Protein−Ligand Binding Affinities. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 21−42. (5) Gallicchio, E.; Levy, R. M. Recent Theoretical and Computational Advances for Modeling Protein−Ligand Binding Affinities. Advances in Protein Chemistry and Structural Biology. Adv. Protein Chem. Struct. Biol. 2011, 85, 27−80. (6) Cheng, T.; Li, X.; Li, Y.; Liu, Z.; Wang, R. Comparative Assessment of Scoring Functions on a Diverse Test Set. J. Chem. Inf. Model. 2009, 49, 1079−1093. (7) Plewczynski, D.; Lazniewski, M.; Augustyniak, R.; Ginalski, K. Can We Trust Docking Results? Evaluation of Seven Commonly Used Programs on PDBbind Database. J. Comput. Chem. 2011, 32, 742− 755. (8) Mikulskis, P.; Genheden, S.; Wichmann, K.; Ryde, U. A Semiempirical Approach to Ligand-Binding Affinities: Dependence on the Hamiltonian and Corrections. J. Comput. Chem. 2012, 33, 1179−1189. (9) Antony, J.; Grimme, S.; Liakos, D. G.; Neese, F. Protein−Ligand Interaction Energies with Dispersion Corrected Density Functional Theory and High-Level Wave Function Based Methods. J. Phys. Chem. A 2011, 115, 11210−11220. (10) Antony, J.; Grimme, S. Fully Ab Initio Protein−Ligand Interaction Energies with Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2012, 33, 1730−1739. (11) Muddana, H. S.; Gilson, M. K. Calculation of Host−Guest Binding Affinities Using a Quantum-Mechanical Energy Model. J. Chem. Theory Comput. 2012, 8, 2023−2033. (12) Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem.Eur. J. 2012, 18, 9955−9964. (13) Wang, R.; Fang, X.; Lu, Y.; Wang, S. The PDBbind Database: Collection of Binding Affinities for Protein−Ligand Complexes with 8083

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084

The Journal of Physical Chemistry B

Article

Prediction of the Full Range of Intermolecular Interactions in Biomolecules. Phys. Chem. Chem. Phys. 2007, 9, 2362−2370. (35) OPENMOPAC. www.openmopac.net (accessed Mar 13, 2013). (36) Klamt, A. The COSMO and COSMO-RS Solvation Models. Wire Comput. Mol. Sci 2011, 1, 699−709. (37) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100. (38) Perdew, J. P. Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B 1986, 33, 8822−8824. (39) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (40) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401−146405. (41) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (42) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104−154123. (43) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic Structure Calculations on Workstation Computers: The Program System Turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (44) TURBOMOLE V6.4 2012, a development of University of Karlsruhe and Forschungszen-trum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007; available from http://www. turbomole.com. (45) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets of Triple Zeta Valence Quality for Atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (46) Eichkorn, K.; Treutler, O.; Ö hm, H.; Häser, M.; Ahlrichs, R. Auxiliary Basis Sets to Approximate Coulomb Potentials. Chem. Phys. Lett. 1995, 242, 652−660. (47) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary Basis Sets for Main Row Atoms and Transition Metals and Their Use to Approximate Coulomb Potentials. Theor. Chem. Acc. 1997, 97, 119−124. (48) Englebienne, P.; Moitessier, N. Docking Ligands into Flexible and Solvated Macromolecules. 5. Force-Field-Based Prediction of Binding Affinities of Ligands to Proteins. J. Chem. Inf. Model. 2009, 49, 2564−2571. (49) Isegawa, M.; Fiedler, L.; Leverentz, H. R.; Wang, Y.; Nachimuthu, S.; Gao, J.; Truhlar, D. G. Polarized Molecular Orbital Model Chemistry 3. The PMO Method Extended to Organic Chemistry. J. Chem. Theory Comput. 2013, 9, 33−45. (50) Genheden, S.; Mikulskis, P.; Hu, L.; Kongsted, J.; Soderhjelm, P.; Ryde, U. Accurate Predictions of Nonpolar Solvation Free Energies Require Explicit Consideration of Binding-Site Hydration. J. Am. Chem. Soc. 2011, 133, 13081−13092.

8084

dx.doi.org/10.1021/jp402719k | J. Phys. Chem. B 2013, 117, 8075−8084