Superior Performance of the SQM/COSMO Scoring Functions in Native Pose Recognition of Diverse Protein-Ligand Complexes in Cognate Docking.

General and reliable description of structures and energetics in protein-ligand (PL) binding using the docking/scoring methodology has until now been elusive. We address this urgent deficiency of scoring functions (SFs) by the systematic development of corrected semiempirical quantum mechanical (SQM) methods, which correctly describe all types of noncovalent interactions and are fast enough to treat systems of thousands of atoms. Two most accurate SQM methods, PM6-D3H4X and SCC-DFTB3-D3H4X, are coupled with the conductor-like screening model (COSMO) implicit solvation model in so-called "SQM/COSMO" SFs and have shown unique recognition of native ligand poses in cognate docking in four challenging PL systems, including metalloprotein. Here, we apply the two SQM/COSMO SFs to 17 diverse PL complexes and compare their performance with four widely used classical SFs (Glide XP, AutoDock4, AutoDock Vina, and UCSF Dock). We observe superior performance of the SQM/COSMO SFs and identify challenging systems. This method, due to its generality, comparability across the chemical space, and lack of need for any system-specific parameters, gives promise of becoming, after comprehensive large-scale testing in the near future, a useful computational tool in structure-based drug design and serving as a reference method for the development of other SFs.

Introduction

In structure-based drug design, docking/scoring is a prime and well-established computational tool. Molecular docking generates ligand geometries bound to the protein (poses), whereas scoring using scoring functions (SFs) ranks them by the predicted affinity (score). Owing to the approximations embodied in docking/scoring methods for the sake of their acceleration, their accuracy has often been compromised.[1] Nevertheless, recent methodological advances made docking/scoring methods an indispensable tool in discovering new protein ligands.[2]

The “docking power” or “sampling power”[3,4] of a docking/scoring method is assessed by its ability to identify the native ligand pose (root-mean-square deviation (RMSD) from the crystal pose <2 Å) in protein–ligand (PL) complexes. Comprehensive testing across diverse PL complexes has shown that in up to 80% of PL complexes this task can be accomplished.[4−8] However, classical SFs had troubles with the identification of the native binding mode as the best-scoring pose (especially in the case of metalloproteins, halogenated ligands, inorganic ligands, etc.).[4] Thus, reliable identification of native PL poses within a diverse set of PL complexes using a single SF remains a challenging task.[3,4,9]

The four major approaches toward scoring are empirical,[10−12] knowledge-based,[8,13,14] statistical/machine learning,[15,16] and physics-based.[17,18] The first three approaches require a training set, and by use of parametrization and statistics, useful models can be obtained.[19] However, because these approaches are dependent on the training set, their predictive power is limited. In contrast, physics-based methods rely on a generally valid description of PL interactions. Traditionally, such approaches were limited to molecular mechanics (MM) methods and simplified variants thereof. Thus, these approaches were inherently limited by the underlying approximations, most importantly the implicit treatment of electrons.

A general solution to the problem of accurately calculating noncovalent interactions in PL systems is the use of quantum mechanics (QM).[20] With QM methods, phenomena of quantum origin, such as charge transfer, are described without further ad hoc parametrization. This is important for systems involving halogen bonding,[21,22] metalloprotein binding,[23] inorganic ligands,[24−26] or covalent bond formation.[27] But because of the high computational demands, QM calculations of sufficient quality (e.g., DFT-D3 level with a triple-ζ basis set) are limited to a few hundred atoms. This limitation can be overcome by use of fragmentation[28,29] or a QM/MM approach.[30−34] Another route is the use of semiempirical QM (SQM). The first QM-based SF was introduced by the Merz group.[35] They combined the Austin model 1 (AM1) SQM method with empirical dispersion (D) and implicit solvation (Poisson–Boltzmann (PB) model). Its validation on a large dataset of PL complexes showed its superior performance, especially for metalloprotein–ligand complexes.[36] Although this was an important pioneering step, the accuracy of the underlying methods, both for the vacuum part (AM1-D)[37] and for solvation (PB),[38] was not sufficient to yield quantitative results. More recently, Sulimov laboratory used the PM7 SQM method[39] in conjunction with the conductor-like screening model (COSMO) implicit solvent model[40] for identification of native ligand poses of 16 PL complexes in cognate docking.[41,42] They showed superior performance of their SQM/COSMO SFs over force-field-based scoring. We should note here that PM7 results for noncovalent interactions can be slightly improved by using the latest version of empirical corrections to the PM6 SQM method (PM6-D3H4X, see below).[43]

In our laboratory, we have been systematically developing empirical corrections to SQM methods to accurately treat an array of noncovalent interactions.[20] The latest version of empirical corrections for dispersion, hydrogen bonding, and halogen bonding yielded the PM6-D3H4X[37,44] method, which, coupled with the COSMO implicit solvent model,[40] forms the core of our SQM-based SF (eq ).[45,46]The score (an estimate of the PL binding free energy) is expressed as an unweighted sum of thermodynamic terms. It consists of the gas-phase PL interaction energy (ΔEint), the change in solvation/desolvation free energy upon complex formation (ΔΔGsolv), the change in the conformation “free” energies of the protein and ligand [ΔG′wconf(P,L)], and the interaction entropy change upon binding (TΔSint).[45,46] The PL complexes are optimized using the solution-phase SQM method before scoring. The ΔEint term is favorable for complex formation and usually is the largest in magnitude. It can reach a few hundreds of kcal/mol for charged or polar ligands. The ΔΔGsolv term opposes binding and can be nearly as large as the first term. These two dominant terms thus partially compensate for each other, and the final score is an order of magnitude smaller. Using this SF (eq ), we have rationalized the binding of series of ligands to a dozen of protein targets,[21,22,46−50] including covalent ligand binding.[27] It should be noted that this SF can also be extended to evaluate explicit solvent effects.[47,48,51]

Recently, we have accelerated our SQM-based SF by considering only the first two dominant terms and replacing the time-demanding SQM optimization with a quick MM relaxation of hydrogens.[23] We have shown in four difficult PL complexes that this SQM/COSMO SF at the PM6-D3H4X level outperforms eight widely used SFs in native ligand pose identification in cognate docking. The number of false-positive (FP) solutions (i.e., those poses that scored better than the native one) was up to 1 order of magnitude lower than that for the classical SFs.[23] In three PL cases, it was even 0. In the challenging system of the tumor necrosis factor-α converting enzyme (TACE) metalloprotein featuring Zn2+ in the active site, which is bound by the thiolate group of the ligand, 39 FPs were found.[23] A major improvement (FP = 0)[52] was observed when the ΔEint term was calculated with a more robust SQM method, the self-consistent-charge density-functional tight-binding method augmented with empirical dispersion (previously shown to be useful for the description of biomolecules)[53] and hydrogen-bonding corrections (SCC-DFTB3-D3H4, abbreviated DFTB3-D3H4).[54] The high-quality description of the other three PL systems was retained.[52] The price for the improvement was a higher but not unsurmountable computational cost. It should be noted that two recent studies used the uncorrected SCC–DFTB method in a QM/MM setup and reported success on the correct ligand binding geometries toward metalloproteins.[31,32] Their approach toward the computationally expensive task was to use a rather small QM part consisting of Zn2+, its coordinating protein side chains, and the ligand on a large number of PL systems.[31,32]

In this study, we aim to validate our SQM/COSMO SFs[23,52] for native pose identification in cognate docking on a data set consisting of 17 PL complexes from five diverse classes, selected using strict criteria for physics-based scoring. We apply two variants of the SQM/COSMO SF (ΔEint term at the PM6-D3H4X or DFTB3-D3H4X level)[23,52] and compare them with four standard SFs (Glide XP,[55] AutoDock4,[56] AutoDock Vina,[57] and UCSF Dock[58]). The performance criterion is the number of FPs[23,52] with an extended definition presented here. We show here that the unique behavior of the SQM/COSMO SFs observed in our recent studies[23,52] hold across 17 diverse PL complexes and gives promise of generality after comprehensive large-scale testing in the near future.

Results and Discussion

Data Set

In this work, we extend our previous pilot studies on four difficult PL systems[23,52] to 17 pharmaceutically relevant and diverse PL complexes from five classes, including three enzyme classes (transferase, hydrolase, and lyase), one chaperone, and two nuclear receptor classes from the PBDbind “core set”[59] (for details, see Methods section). We apply the strict criteria needed for physics-based scoring. Specifically, the crystal structures of the complexes have resolutions better than 2.5 Å, well-resolved electron densities for the ligands, and protein active sites. The ligands have variable chemistries, sizes (molecular weight of 305–666 Da), charges, and flexibilities (for details, see Methods section). Their binding constants toward their targets range from micro- to picomolar.

The crystal poses of the ligands were scored as reference. The ligand poses generated previously by docking[4] with seven docking programs (for details, see Methods section) totaled 4566 poses. RMSD-based clustering (see Methods section) of the poses was carried out to avoid pose redundancy. After this, the number of poses decreased to 3328, corresponding to approximately 250 poses per target. A comprehensive evaluation of the recognition of near-native poses requires a balanced distribution of RMSDs of the docked ligand poses with respect to the crystal (native) geometry from very similar to very dissimilar (up to 10 Å). In most of the PL systems, 20–60% of poses had RMSD <2 Å (Figure S1A). Furthermore, poses were evenly distributed in RMSD ranges of 2–5 and 5–10 Å (roughly 20–40% for each category). Figure S1B shows minimal RMSD (RMSDmin) for the poses studied. Near-native poses within the experimental accuracy of X-ray crystallography of 0.5 Å[60,61] were found in all but two cases (10GS: RMSDmin = 0.85 Å and 2VOT: RMSDmin = 0.70 Å). However, these cases only slightly exceeded the threshold.

Scoring

For each of the six SFs (two variants of SQM/COSMO SF and four standard SFs; for details, see Methods section), the scores of the docked ligand poses in their respective target proteins were calculated, transformed to relative scores with respect to the score of the crystal pose, and normalized (see Methods section).

The overall sampling power of all the SFs is shown as the enrichment plot (Figure S2), that is, the percent of PL cases (y axis) in which the best-scoring ligand of a given SF has defined RMSD (x axis) to the crystal pose. In the standard range of RMSD up to 2 Å, the SQM/COSMO SFs at the DFTB3-D3H4X level perform the best (88% of PL systems), followed by SQM/COSMO at the PM6-D3H4X level together with UCSF Dock (82% of PL systems). Slightly worse is the performance of Glide XP (76%), followed by AutoDock4 (71%), and AutoDock Vina (65%) (Figure S2). In recognition of near-native poses (RMSD < 0.5 Å), the two SQM/COSMO SFs together with AutoDock Vina perform the best (47%), followed by AutoDock4 and UCSF Dock (41%) and Glide XP, which recognize the poses only in 29% of cases.

The SQM/COSMO SFs also had the lowest number of PL complexes (two for DFTB3-D3H4X/COSMO, Table S1) for which the best-scoring pose exceeded the threshold for success of 2 Å. This was closely followed by SQM/COSMO at the PM6-D3H4X/COSMO level, UCSF Dock, and Glide XP (three cases). Five and six failures were found for Autodock4 and AutoDock Vina, respectively (Table S1). Averaging the RMSDs of the best-scoring poses across all 17 PL complexes (and counting all the failures >2 Å as 2.1 Å), DFTB3-D3H4X/COSMO was the winner (0.71 Å), closely followed by PM6-D3H4X/COSMO and UCSF Dock (0.77 and 0.79 Å, respectively; Table S1). Worse results (around 1 Å) were obtained for AutoDock Vina, AutoDock4, and Glide XP.

For detailed performance evaluations, we use the number of FP solutions criterion[23,52] with an extended definition presented here. Previously, FPs were defined as those poses that scored better than the native pose (defined by a 0.5 Å RMSD cutoff from the crystal pose due to inaccuracies of crystal structures).[23,52] Here, we allow room for larger uncertainties of native pose recognition by defining “hard FPs” (HFP) in which the cutoffs were increased to RMSD >2 Å and score better than −1 kcal/mol. The RMSD cutoff now also includes the effects of flexible parts of the ligands sticking out to the solvent, and the score cutoff corresponds roughly to 2–3 kcal/mol of unscaled energies, which are rough error bounds of the physics-based method. The HFPs for AutoDock Vina, AutoDock4, Glide XP, and UCSF Dock were high—211, 350, 425, and 635, respectively (Figure A). The SQM/COSMO SFs performed much better with the numbers of HFPs being up to 1 order of magnitude smaller—40 and 42 for the DFTB3-D3H4X and PM6-D3H4X levels, respectively (Figure A).

Figure 1

Number of HFP solutions for the six SFs used here across all the 17 PL systems studied. (A) Number of HFPs and (B) HFPs for individual PL complexes sorted by ligand charge: neutral (left) and charged (right).

The number of HFPs for individual PL complexes (Figure B and Table S2) differed markedly with respect to the ligand charge: in the case of the neutral ligands (Figure B, left), they were by 1 order of magnitude smaller than that for the charged ones (Figure B, right). For SQM/COSMO at the PM6-D3H4X and DFTB3-D3H4X levels, the numbers of HFPs for neutral ligands were single-digit values (1 and 2, respectively). The classical SFs performed worse, with the number of HFPs ranging from 18 to 85 for neutral ligands (Figure B, left and Table S2). The complex with the largest number of HFPs was the RNA-dependent RNA polymerase/ligand complex (3GNW) with 71, 28, and 14 HFPs calculated with AutoDock4, UCSF Dock, and AutoDock Vina, respectively. A large number of HFPs (40) was also observed for the cyclin-dependent kinase 2 (CDK2)/ligand complex (2FVD) for Glide XP (Table S2).

The results show that the classical SFs had larger troubles in identifying the native binding poses for charged ligands (for the classical SFs, more than 90% of HFPs were found for charged ligands). The largest number of HFPs (140) was found with Glide XP for the α-l-fucosidase (2ZX6) PL complex, which had a positively charged ligand. For UCSF Dock, four systems, 2P4Y, 4GID, 2VOT, and 3NOX, yielded in total 403 HFPs, which is 70% of HFPs for the charged ligands in that method (577; Table S2). In contrast, the number of HFPs for the charged ligands for the SQM/COSMO was in total 38 and 41 for DFTB3-D3H4X and PM6-D3H4X, respectively. This is considerably lower than the classical SFs (193–577 HFPs) (Table S2). For SQM/COSMO at the DFTB3-D3H4X level, the largest number of HFPs was 20 and 8 for 2P4Y and 3NOX, respectively. Also, PM6-D3H4X/COSMO had some troubles with these systems (5 and 10 HFPs, respectively). In both 2P4Y and 3NOX complexes, the HFP poses have the ligand cores placed at very similar positions as the crystal pose, whereas moieties sticking out to the solvent (the benzisoxazol and morpholino groups, respectively) had fewer noncovalent interactions with the protein. This can be one reason why poses with higher RMSD could score well. Other reasons can be some of the approximations embedded in our protocol for speed, such as the neglected terms in the SQM/COSMO SF (change of conformational energy, entropy) or explicit water molecules, which may need to be included in some PL systems for reliable description of the energetics.[47,51]

Conclusions

The sampling (docking) power, that is, the ability to recognize a ligand native pose in cognate PL docking, of two variants of quantum-mechanics-based SQM/COSMO SFs is tested here on 17 PL systems from five diverse protein families carefully selected for physics-based SFs. For comparison, four standard SFs—Glide XP, AutoDock4, AutoDock Vina, and UCSF Dock, are used. The SQM/COSMO SFs at the PM6-D3H4X and DFTB3-D3H4X levels markedly outperform the standard SFs as judged by the number of HFP poses. The time requirements for the SQM/COSMO SF (Table S3) are higher than those for classical SFs, but given the supercomputer power, thousands of docking poses can be evaluated in a reasonable time. The results of the freely available SQM/COSMO SFs give promise of generality, and after comprehensive large-scale testing in the near future, this method could serve as a useful tool in structure-based drug design and reference for SF development.

Methods

QM-based interaction energy calculations require sensible geometries and, therefore, we needed good-quality structures of PL complexes. The crystallographic structures should have fair resolution (<2.5 Å) with fully resolved electron density for the entire ligand and surrounding binding site residues. These criteria are fulfilled by the docking/scoring benchmark set PDBbind core set.[3,59,62] In our study, 17 PL complexes (Figure and Table ) were used with targets from diverse protein families: three enzyme classes (transferase, hydrolase, and lyase), chaperone, and nuclear receptor (Table ). The ligand structures are shown in Figure .

Figure 2

Two-dimensional structures of the ligands studied.

Table 1

Summary of the 17 PL Complexes Studied

PDB code	resolution (Å)	protein name	class	ligand charge	rotatable bonds in ligand
2FVD	1.8	CDK2	transferase (E.C.2)	0	6
10GS	2.2	glutathione S-transferase		–1	13
3PE2	1.9	casein kinase IIα		–1	4
3GCU	2.1	mitogen-activated protein kinase 14		0	6
2OBF	2.3	phenylethanolamine N-methyltransferase		+1	4
3JVS	1.9	checkpoint kinase 1		–1	5
3GNW	2.4	hepatitis C virus NS5B RNA-dependent RNA polymerase		0	5
2CET	1.9	β-glucosidase A	hydrolase (E.C.3)	+1	4
4GID	2.0	β-secretase I		+1	16
2ZX6	2.4	α-l-fucosidase		+1	4
3NOX	2.3	dipeptidyl peptidase 4		+1	3
2VOT	1.9	β-mannosidase		+1	4
2XB8	2.4	3-dehydroquinate dehydratase	lyase (E.C.4)	–1	4
2VW5	1.9	heat shock protein Hsp82	chaperone	0	3
2YKI	1.6	heat shock protein Hsp90-α		0	3
2P4Y	2.2	peroxisome proliferator-activated receptor γ	nuclear receptor	–1	9
3G0W	1.9	androgen receptor		0	2

Docking Poses

Ligand poses obtained by seven commonly used docking programs were collected from previously published work.[4] These programs were AutoDock (version 4.2.6),[56] AutoDock Vina (version 1.1.2),[57] LeDock (version 1.0),[63] UCSF Dock (version 6.7),[58] Glide SP (version 67011),[55] Glide XP (version 67011),[55] and Surflex Dock (version 2.706.13302).[64] For each target, the ligand poses were pooled, which amounted to approximately 350 poses per target. To reduce the redundancy, all poses per target were clustered using the “cluster_conformer” script in the Schrödinger suite[65] with an RMSD cutoff of 0.5 Å. The number of poses was thus reduced to approximately 250 poses per PL system. Each ligand pose, as well as X-ray reference geometry, was scored “in-place” using four classical SFs (AutoDock,[56] AutoDock Vina,[57] UCSF Dock,[58] and Glide XP[55]) and compared to that of two variants of SQM/COSMO SFs, see below.[23,52]

Protein and Ligand Preparation

Protein Structure Preparation

Following the standard virtual screening protocol,[4] all the crystal waters were removed from the PL complexes. As noted previously,[23] physics-based SFs require special care in preparing the PL structures. For all proteins, which were not deposited as monomers, chain A was used for protein preparation except 10GS and 2XB8 where the dimer interface makes important contributions to the binding. We used the LEaP program, which is part of the AMBER14 suite,[66] to protonate the proteins. The protonation state of histidine residues was assigned manually on the basis of hydrogen-bonding patterns. Cysteine disulfide bonds were assigned manually on the basis of the sulfur–sulfur distance. In the case of the 3PE2 complex, the B conformation of M163 was used because it forms interactions with the ligand. Hydrogen atom positions in PL complexes were relaxed by the simulated annealing protocol using short molecular dynamics (MD) (for details see Supporting Information).

Ligand Preparation

The protonation states of the ligands were carefully checked by pKa calculations at pH 7 using Schrödinger “Propka”.[65] The collected docking poses from seven different programs had different output file formats. Each ligand was made into one common MOL2 file format without any changes in X, Y, and Z coordinates. Partial charges were derived at the AM1-BCC level using RESP.[67−69]

RMSD Measurements

RMSD values of all the ligand poses were calculated with respect to the corresponding X-ray geometry of the ligand (without any further optimization) with the “heavy atom” option using the “rmsd.py” script by Schrödinger.

SQM/COSMO SFs

All of the docked PL complexes with close contacts (cutoff of 1.5 Å) between the protein and the ligand were relaxed by short AMBER/GB optimization as in previous studies.[52] Next, optimal hydrogen positions were localized in each complex using a short MDs run using AMBER/GB as in our previous studies.[52] The SQM/COSMO score is a sum of ΔEint and ΔΔGsolv terms. For speed-up and without compromising the reliability,[23] the former term was calculated on large parts of the protein (typically the ligand plus 10 Å protein surroundings) using two approaches: (i) the corrected PM6-D3H4X[44] and (ii) DFTB3-D3H4X method, a third-order DFTB[70,71] with the 3OB parameter set[72,73] and the latest version of the D3H4X corrections for noncovalent interactions.[74] The solvation free energy was calculated on the same truncated system as above using a COSMO implicit solvent model at the PM6 level.[40]

Glide XP Score

All scoring calculations were performed with Glide XP[55] and run in the extra precision (XP) workflow framework. Docking grids were generated by Glide using the cocrystallized ligand at the center of the grid box. The compounds were scored with the option “score in place only”.

AutoDock4 and AutoDock Vina

For both AutoDock4[56] and AutoDock Vina,[57] the centers of grid boxes were arranged according to the centers of the crystal ligand poses. The grid box sizes were adjusted to make scoring possible for all combinations of ligands and conformations. AM1-BCC RESP partial charges were used.

UCSF DOCK

The grid spacing was 0.3 Å. The cutoff for nonbonded interactions was not used. We used AMBER parameters. For ligands, we used AM1-BCC RESP partial atomic charges.

Score Scaling

The scores of all the poses of the 17 PL complexes obtained by the 6 SFs were transformed into relative numbers with respect to the score of the crystal pose and normalized as done previously.[23,52]