Conformation Search Across Multiple-Level Potential-Energy Surfaces

Article Cite This: J. Chem. Theory Comput. 2019, 15, 4264−4279

pubs.acs.org/JCTC

Conformation Search Across Multiple-Level Potential-Energy Surfaces (CSAMP): A Strategy for Accurate Prediction of Protein− Ligand Binding Structures Lin Wei,†,‡ Bo Chi,†,‡ Yanliang Ren,†,‡ Li Rao,*,† Jue Wu,† Huan Shang,† Jiaqi Liu,† Yiting Xiao,† Minghui Ma,† Xin Xu,*,§ and Jian Wan*,†

Downloaded via GUILFORD COLG on July 18, 2019 at 12:53:34 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

International Cooperation Base of Pesticide and Green Synthesis (Hubei), Key Laboratory of Pesticide & Chemical Biology (CCNU), Ministry of Education, Department of Chemistry, Central China Normal University, Wuhan 430079, China § Shanghai Key Laboratory of Molecular Catalysis and Innovative Materials, Ministry of Education (MOE) Laboratory for Computational Physical Science, Department of Chemistry, Fudan University, Shanghai 200433, People’s Republic of China S Supporting Information *

ABSTRACT: Accurate protein binding structure determination presents a great challenge to both experiment and theory. Here, in this work, we propose a new DOX protocol which combines the ensemble molecular Docking as the coarse-level, structure Optimization with the semiempirical quantum mechanics methods as the medium level, and the eXtended ONIOM (XO) calculations as the fine level. The fundamental of the DOX protocol relies on the Conformation Search Across Multiplelevel Potential-energy surfaces (CSAMP) strategy, where the conformation spaces of a funnel-like structure are searched from the coarse level with hundreds of candidates to the medium level with around 10 top candidates to the fine level with the final top 1 or 2 binding modes. An in-depth test for the protocol set up against 28 crystallographic data consisting of HMGR-statins, SDase-inhibitors, 3HNRase-inhibitors, and NA-inhibitors yielded a satisfactory result with ∼0.5 Å root-mean-square deviations (RMSDs) on geometries and ∼0.8 kcal/mol absolute error of relative binding energies on average. A further larger scale validation on the Astex test set (including 85 diverse structures) revealed an impressive performance with a RMSD < 2 Å success rate of 99%, suggesting DOX is a promising computational route toward accurate prediction of the protein−ligand binding structures. cheapness.10,11 Recently,10,12 it has been noticed that the docking algorithm for conformation search is sufficiently good to generate a set of binding structures which encompasses the reference crystal structure; nevertheless, the employed empirical scoring functions often fail to acknowledge the “right poses” and credit a significantly worse binding pose as the optimal docking result. Hence, in our previous work, we introduced a protocol where some high-level quantum chemistry methods replaced the empirical scoring function for the ranking of dockinggenerated binding poses to predict the optimal binding structure. The resulting DOX method12 is a sophisticated combination of molecular Docking,13 ONIOM3,14 and eXtended ONIOM (XO)15,16 calculations. In the first step of the DOX calculation, Surflex-docking@SYBYL13 computation is used to generate a binding conformation space that consists of 100 protein−ligand binding structures rather than only one top score binding pose. In the next step of the DOX

1. INTRODUCTION Binding interactions (noncovalently and reversible covalently) are ubiquitous in chemical and biological systems involving molecular recognition and assembly, protein−protein bindings, protein/DNA/RNA-ligand bindings, etc., which are strongly responsible for the physical, chemical, biological, and physiological functions of these systems.1−6 In particular, the primary goal of drug research and development is to find a drug-like molecule that binds strongly to a given receptor cavity. Hence, it is very beneficial if the protein−drug binding structure could be accurately determined. Such information is a basic prerequisite for the structure-based drug design. Currently people solely rely on crystal structure for accurate protein−drug binding structure determination. Nonetheless, protein crystallization with X-ray structural analysis is extremely difficult and expensive or even impossible. Lack of information on crystallographic binding structures can result in failure on virtual screening or modification and optimization of the potential lead structures. A useful alternative is provided by molecular docking methods,7−9 although there is always a trade-off from accuracy to computational speediness and © 2019 American Chemical Society

Received: November 13, 2018 Published: May 29, 2019 4264

DOI: 10.1021/acs.jctc.8b01150 J. Chem. Theory Comput. 2019, 15, 4264−4279

Article

Journal of Chemical Theory and Computation calculation, the ONIOM314 (ωB97X-D17/6-31G*:PM6:18AmAmber19) single-point energy calculation is performed to rerank the entire binding conformation space from docking, where the top 10 binding poses credited by ONIOM3 come into the short list. In the last step of the DOX calculation, each of the top 10 structures is optimized at the XO2 (ωB97X-D/631G*:PM6) level of theory, and the final optimal binding structure is determined according to the order of the relative binding energy (RBE) calculated at the XO2 (ωB97X-D/6311+G**:PM6) level on top of the XO2-optimized geometry. A test against 15 crystal structures of HMGR (3-hydroxy-3methylglutaryl coenzyme-A reductase)-statin complexes showed that DOX can predict the right geometries for all HMGR-statin complexes without exception.12 The RMSD (root-mean-square deviation) between the DOX-predicted optimal binding pose and the atomic coordinates of the ligand molecule in the corresponding crystal structure of the binding complex is merely 0.54 Å on average, suggesting DOX as a promising computational route toward predicting the accurate protein−ligand binding structures, which can be used as an effective alternative to the experimental counterpart. Furthermore, it has been demonstrated that the DOX-predicted binding structures are also energetically accurate with an average energy difference of 1.13 kcal/mol as compared to those calculated at the crystal binding modes. On the other hand, the Surflex-docking@SYBYL-predicted top binding modes are on average ∼6.5 kcal/mol less stable. It has to be noted that a common practice in docking assessment studies is that the initial structure used for docking was taken from the bound ligand structure extracted from the corresponding crystal structure of the protein−drug complex, while such a choice of the initial structure was found to facilitate the generation of a better ensemble of binding conformations by the docking programs.20,21 Nevertheless, the bound ligand structure is not available for docking in a practical application when a theoretical prediction is needed for virtual screening or ligand optimization, because the corresponding crystal structure of the binding complex is absent. Our previous DOX protocol has also benefited from the bound ligand initial structure. Thus, in our continued study on the DOX method we employed the free ligand structure optimized by the molecular mechanics (MM) as the docking input. The performance of the previous DOX degraded, where the averaged RMSDs on geometries and energy differences increased to 0.89 Å and 2.66 kcal/mol for the same set of HMGR-statin complexes (see Table S8 in Supporting Information for more details). Hence, in this work, we carried out systematic validations and detailed analyses on the basic strategy, resulting in a significantly improved new DOX protocol which combines the ensemble molecular Docking, structure Optimization at the semiempirical quantum mechanics level and the XO calculations. The new protocol represents a promising accuracy at more realistic and strict conditions for structure-based drug design, where an extended test against crystallographic data of 112 protein−ligand binding complexes has been performed. Without relying on the bound ligand structures, the average RMSD of the present DOX prediction is just 0.68 Å and the RMSD < 2 Å success rate is as high as 99% for the whole set. Furthermore, the RMSD < 1 Å success rate is 79%, well beyond the reach of traditional molecular docking techniques.

The rest of the paper is organized as follows. In section 2 we present the details of the strategies and methods. The results and discussion are presented in section 3. The conclusion is summarized in section 4.

2. STRATEGIES AND METHODS 2.1. CSAMP Strategy. Conformation Search Across Multiple-level Potential-energy surfaces (illustrated in Figure 1), CSAMP in the following contents, is the fundamental of

Figure 1. Schematic for the CSAMP (conformation search across multiple-level potential-energy surfaces) strategy employed in the DOX protocol.

the DOX protocols, where three levels of theories, coarse, medium, and fine, are jointly used to locate the optimal binding pose. The CSAMP strategy is employed, aiming to approach the optimal binding pose accurately yet efficiently. We would like to emphasize that efficiency is not a synonym of speediness and cheapness. Instead, it counts on how many successful hits to achieve the desired goal at each level in an efficient way. The coarse level in the DOX protocols invokes docking at the MM level, which is very efficient so that sizable conformations across the MM potential energy surface (PES) can be searched with feasible computational cost, while the accuracy just needed is to make sure that the conformation space from the coarse level contains the real global minimum. The medium level narrows down the conformation space. The computational cost of the medium-level theory should be low enough to allow for efficient ranking in the conformation space from the coarse-level theory, while its accuracy needed is to make sure the right one is still included in a couple of top candidates at the medium level. Finally, the requirement for the fine-level theory mainly lies in its accuracy. Ideally, the finelevel theory ranks the right binding structure as the top 1 binding pose from the top candidates at the medium level. The above discussion reveals a fact that although the performance of the DOX protocol is a joint contribution from all three levels of theory, one should validate the protocol by checking if each level of theory reaches its own goal. Although combination of multiple PES search is a frequently used 4265


Article

Journal of Chemical Theory and Computation

Figure 2. Schematic illusitration of the computational models employed in DOX calculations as the medium level or the fine level: (a) ONIOM3/ SP, (b) PM7/OPT, and (c) XO2/SP//OPT.

level in the previous DOX protocol, i.e., 100 poses generated using the docking procedure with the bound ligand structure, denoted as Surflex (Bound1,100). Recall that the goal of docking here is to generate a binding conformation space that is as robust as possible so that the right one is most likely to be included. As Surflex(Confort1,100) does not lead to satisfactory results (Table S8), we expect Surflex(Confort6,300) from the ensemble docking procedure to be better (see below). Plewczynski et al. found that the ensemble docking procedure dramatically increased the possibility to include a “Good” pose (RMSD < 1 Å) in the generated binding pose conformation space in comparison with the single-docking procedure.10 An obvious explanation is that using more input structures helped in searching the conformation space more extensively. However, they also found that docking top score predictions did not benefit from such advantage. On the contrary, the CSAMP strategy allows DOX to take full advantage of the ensemble docking procedure. Further discussion concerning the ensemble docking procedure could be found in the Supporting Information (Table S1). 2.3. Medium Level: Details of ONIOM3/SP and PM7/ OPT. In our previous DOX protocol, the medium level is the single-point energy evaluation with ONIOM3(ωB97X-D/631G*:PM6:Amber) on top of the docking-generated binding structures from Surflex(Bound1,100). It is abbreviated as ONIOM3/SP in the following contents as schematically illustrated in Figure 2a. ONIOM3/SP aims to narrow the conformation space down to 10 top binding poses based on the calculated RBEs. In the continued study, we found that single-point energy evaluation on top of the Surflex(Confort1,100) results often does not reflect the nature of PES owing to the low quality of docking-predicted geometry parameters. However, considering the computational cost, geometry optimization at the ONIOM3 level is not allowed at this stage and at an extended conformation space from Surflex(Confort6,300). Hence, optimization at the PM7 level of theory (PM7/OPT for short) is introduced as a replacement of the previously used ONIOM3/SP method. The PM7

strategy in almost every topic of computational chemistry, detailed analysis for its own goal is seldom seen in the literature, while the above discussion could be useful in all similar situations. Here the CSAMP strategy is applied to address the challenge for predicting the accurate protein−ligand binding structure when the bound ligand structure is unavailable as the initial input in the real situation for the structure-based virtual screening and lead optimization. 2.2. Coarse Level: Details of Molecular Docking Calculations. In this work, Confort22 was used to generate the ligand initial structures for molecular docking without relying on its bound structure. Confort is a program implemented in SYBYL software for the conformation search of small organic compounds.21 Confort can be used to identify global minimum energy conformation or an ensemble of conformations within a specified energy range around the global minimum. An optimized and finely tuned version of the Broyden−Fletcher−Goldfarb−Shanno (BFGS) minimization algorithm23−25 is used for geometry optimization. The energies and gradients needed are obtained by the Tripos Force Field (TFF).26 All of the docking calculations were performed using the Surflex-docking@SYBYL module,13 the same as in our previous work.12 Before proceeding to the molecular docking procedure, all water molecules are removed and, unless they belong to a cofactor of the protein, all other nonprotein atoms are also removed. Since molecular docking results are sensitive to the input structures, we tested the quality of binding conformation space generated by using just one Confort ligand global minimum structure for docking input or using the first six ligand conformations. The latter is referred to as the ensemble docking procedure in the literature.10 The docking procedure employing only one Confort-generated ligand structure is referred to as Surflex(Confort1,100) with 100 poses generated in one docking calculation, while the docking procedure employing six ligand conformations is referred to as Surflex(Confort6,300) with 50 poses generated in each of the six docking calculations. This can be compared with the coarse 4266


Article


Figure 3. Comparison of the old (a) and new (b) DOX protocols.

residues within 3 Å from any atom of the ligand, which usually makes a cluster of more than 300 atoms. Further discussion concerning the range of the high-level region can be found in the Supporting Information (Table S5). All charged groups are neutralized as an empirical approximation to account for the solvation effect of water solution that strongly weakens the Coulomb interactions between charged groups; see Supporting Information Part 2 for a detailed illustration. It is a successful approach used extensively in QM/MM studies (e.g., according to Friesner et al.32) and our previous works employing the XO methods.12,30 For better efficiency, the high-level region is further divided into several fragments, each of which consists of the entire ligand molecule and a few protein residues surrounding one of the ligand function groups. Detailed graphic illustration for the XO2 calculations can be found in Figures S2−S5 in the Supporting Information. While each of the top 10 structures credited by the medium level is reoptimized at the fine level of XO2 with the 6-31G* basis set for DFT, the final ranking is determined according to the order of RBEs calculated with XO2 using the 6-311+G** basis set. An illustration of RBE could be found in the Supporting Information. Further discussion concerning the geometry optimization at the XO2 level can also be found in the Supporting Information (Table S6). Previously,12 RBE was calculated using the crystal binding structure as the reference, which is often unavailable in practice. In this work, we choose the binding pose with the lowest total energy (i.e., Top Pose) as the reference so that the RBE of a binding pose weighs how much it is higher than the predicted most stable binding pose in energy. Although the zero point itself corresponds to the optimal binding structure, we recommend that it is worthwhile to consider those binding poses with RBEs less than 2 kcal/ mol in the drug design. The original DOX protocol12 involves a sequence of molecular Docking, ONIOM3/SP, and XO calculations, which has been proved to work satisfactorily starting from the bound initial structure. The new DOX protocol proposed here invokes a sequence of ensemble docking, structure

semiempirical theory implemented in MOPAC2016 was claimed to be the most accurate method in the MOPAC family as a state-of-the-art semiempirical theory.27 In particular, PM7 exhibits a good performance on the description of nonbonded interactions with a reported MAD (mean absolute deviation) of 0.78 kcal/mol on the S66 intermolecular interaction energy set.18 Further discussion concerning the choice of semiempirical method for medium-level calculation can be found in the Supporting Information (Tables S2 and S3). Note that the optimized system only involved the ligand molecule and protein residues surrounding it (namely, the ligand binding pocket), while only ligand atoms were allowed to relax during optimization, as shown in Figure 2b. 2.4. Fine Level: Details of XO2, OPT, and SP. By introducing the divide and conquer strategy28 and inclusion− exclusion principle,29 the eXtended ONIOM (XO) method proposed by Xu et al.15,16 enables an accurate and efficient description of a very large quantum mechanics (QM) region. Previously, the XO method has been shown to exhibit promising performance on the protein−drug binding affinity estimation and the protein−drug binding structure predictions.15,16,30 In this work, we continue to use the XO method for the fine-level treatment of the binding structure as in the original DOX protocol. Currently, the XO method is implemented using the Gaussian external facilities. Hence, all calculations concerning XO are performed using the Gaussian09 package.31 As shown in Figure 2c, for each protein−ligand complex, the ligand molecule and protein residues within 7 Å from the ligand is considered in a two-level XO (XO2) calculation. The computation system which usually includes approximately 1000 atoms is first divided into a high-level region treated with a suitable density functional theory (DFT) method, e.g., ωB97X-D,17 and a low-level region treated with semiempirical PM618 (since PM7 is not available in the Gaussian09 package). Further discussion concerning choice of the low-level theory could be found in the Supporting Information, Table S4. The high-level region consists of the ligand molecule and all 4267


Article

Journal of Chemical Theory and Computation Table 1. Results of Different Docking Procedures: Redocking Test on 15 HMGR-statin Complexesa Surflex(Bound1,100)b Top Pose g

e

PDB ID

score

RMSD

1HW8 1HW9 1HWI 1HWJ 1HWK 1HWL 2Q6B 2Q6C 2R4F 3BGL 3CCW 3CCZ 3CD0 3CD5 3CD7 mean

13.83 14.87 13.33 16.96 17.89 14.84 16.54 15.71 15.95 15.24 15.53 14.16 14.81 16.46 16.95 15.54

1.75 0.82 0.64 1.06 0.47 0.73 0.69 1.06 0.50 0.79 0.97 0.75 0.83 1.70 0.93 0.91

Surflex(Confort1,100)c

f

Best Pose h

Top Pose I

Surflex(Confort6,300)d

Best Pose

Top Pose

Best Pose

score

RMSD

Good

score

RMSD

score

RMSD

Good

score

RMSD

score

RMSD

Good

12.56 12.22 13.33 16.21 17.43 14.39 16.04 13.26 14.90 15.01 14.94 13.81 14.81 15.50 14.82 14.62

0.89 0.28 0.64 0.48 0.40 0.45 0.56 0.54 0.46 0.77 0.84 0.73 0.83 0.49 0.54 0.59

41 33 19 79 44 58 37 17 42 13 5 34 2 62 25 34.1

11.00 13.79 11.90 15.15 14.20 13.40 14.65 13.79 13.18 14.20 14.46 13.77 13.68 15.54 14.15 13.79

1.92 1.38 1.17 0.91 0.59 0.95 1.23 1.12 0.87 0.99 0.55 0.64 0.43 0.56 1.35 0.98

9.73 13.34 10.64 12.78 14.20 13.40 14.16 11.11 11.56 11.72 11.91 13.31 13.68 14.92 12.54 12.6

0.87 0.85 0.58 0.58 0.59 0.95 1.01 0.60 0.76 0.74 0.55 0.56 0.43 0.53 0.53 0.68

2 2 19 23 36 1 0 8 7 20 11 47 4 82 45 20.5

12.87 13.14 11.84 15.05 15.07 13.00 15.52 14.28 15.08 14.38 15.06 14.02 14.14 15.60 15.51 14.30

1.93 1.76 0.89 0.84 0.55 0.74 1.00 1.05 0.43 0.93 1.31 0.59 0.60 1.81 1.15 1.04

11.23 11.84 8.54 13.31 14.49 11.96 14.98 11.03 14.89 13.32 12.73 12.85 13.43 14.22 14.50 12.89

0.96 0.80 0.54 0.55 0.48 0.66 0.53 0.48 0.43 0.38 0.38 0.53 0.45 0.47 0.43 0.54

2 11 66 74 115 28 136 65 38 156 64 66 65 168 99 76.9

a

RMSD Unit: Angstroms. bOne hundred binding poses are genrated in one docking calcualtion with initial ligand structure taken as the bound ligand structure extracted from the corresponding complex crystal structure. cOne hundred binding poses are genrated in one docking calcualtion with initial ligand structure taken as the ligand global minimum structure obtained by the Confort@SYBYL program. dThree hundred binding poses are generated in six docking calculations with six initial ligand structures obtianed by the Confort@SYBYL program. eMolecular dockingpredicted optimal top score binding pose. For Confort(6,300), the top pose is taken as the top score binding pose in all of the 300 poses generated in six docking calculations. fBest pose is the binding pose most similar to the crystal binding structure, credited by RMSD, in all of the 100 or 300 binding poses. gBinding affinity of the corresponding binding pose scored by empirical scoring function. hRMSD between the predicted binding pose and the bound ligand geometry coordinates in complex crystal structure. INumber of binding poses with RMSD less than 1.0 Å generated in corresponding docking calculations.

interactions and thus the complex structure. We refer it as the “bound ligand redocking”. In order to suggest improvements to ligands or screen ligand databases for candidate binders, a theoretical prediction of the protein−ligand binding structure has to start from a free ligand. We refer the procedure where a free ligand is docked back to the protein structure extracted from the crystal structure of the corresponding protein−ligand complex as the “free ligand redocking”. In this work, free ligand redocking tests have been carried out on all 112 crystal structures, where the initial ligand structures are the low-energy structures obtained by Confort@SYBYL, which are usually quite different from the bound state. Without relying on the “correct” bound state conformation of the ligand, molecular docking calculations are expected to give worse results,10 which pose a challenge to the subsequent DOX calculations. In real condition of drug design where the crystallographic data are often not available for the binding complex of interest, the desired ligand has to be docked to a protein structure extracted from the other “same protein−different ligand” complexes, or an apo-state protein structure, or even protein structures generated by homology modeling. In this work, we also performed the cross-docking tests against the HMGRstatin complexes for validation of the new DOX protocol. In such cross-docking tests, the protein structure extracted from the 3CCZ crystal structure is used for prediction of the rest of the 14 HMGR-statin complexes, while the protein structure extracted from 2R4F is used for the 3CCZ ligand. As the protein structure of the important amino acid residues surrounding the ligand is optimized in the last stage of the DOX protocol via XO2/OPT, the flexibilities of the protein− ligand binding poses can be better taken into account (see

optimization at PM7 and XO calculations, starting from an ensemble of free ligand conformations. Both the old and the new DOX protocols (shown in Figure 3) represent the conformation space of a funnel structure from coarse to medium to fine levels based on the CSAMP strategy. 2.5. Testing Set and Validation Procedure with Redocking and Cross-Docking. The primary testing set used in this work includes a total of 28 protein−ligand complexes, consisting of 15 HMGR-statins (PDB ID 1HW8, 1HW9, 1HWI, 1HWJ, 1HWK, 1HWL, 2Q6B, 2Q6C, 2R4F, 3BGL, 3CCW, 3CCZ, 3CD0, 3CD5, 3CD7), 5 SDase (scytalone dehydratase) inhibitors (PDB ID 3STD, 4STD, 5STD, 6STD, 7STD), 4 3HNRase (trihydroxynaphthathalene reductase) inhibitors (PDB ID 1YBV, 1G0O, 1DOH, 1G0N), and 4 NA (Neuraminidase) inhibitors (PDB ID 2HTU, 2HTQ, 1F8B, 4KS1) whose crystal structures have been well examined in the related published and ongoing studies of our group.33,34 For a full-scale validation of its accuracy, efficiency, universality, and practicality, the new DOX protocol has also been tested on a well-established data set, the Astex diverse set.21 The Astex set contains crystallographic data of 85 protein−ligand complexes. All proteins are of direct interest to the pharmaceutical or the agrochemical industry, while all ligands are drug-like molecules. Note that the Astex set also includes entry 1HWI mentioned above, therefore the two sets contain 112 protein−ligand complexes together. “Redocking” often refers to the procedure where the ligand is docked back to the protein structure, both of which are extracted from the crystal structure of the corresponding protein−ligand complex. This procedure verifies the suitability of the docking method to recover the known binding 4268


Article


Figure 4. Comparison of docking and DOX performance on “free ligand redocking” test of 28 HMGR, SD, NAm and 3HNR systems. (a) Superposition of Surflex(Confort6,300)-predicted Best Pose (red) versus the bound ligand structure in the corresponding complex crystal (gray). (b) Superposition of the DOX-predicted ligand binding structure(cyan) versus the bound ligand structure in the corresponding complex crystal (gray). Data in parentheses are the corresponding RMSDs for the 28 systems.

As shown in Table 1, Surflex(Bound1,100), (Confort1,100), and (Confort6,300) lead to Top Pose with an average RMSD of 0.91, 0.98, and 1.04 Å, respectively, for the redocking tests on 15 HMGR-statin complexes. These may all seem reasonable as compared to our standard with RMSD < 1.0 Å.12 However, according to the CSAMP strategy, the goal of molecular dockings is to provide a wide funnel top at the coarse level, whose accuracy just needed is to encompass the right binding poses in the sampling space. For such a purpose, the Best Pose provides the suitable information. As shown in Table 1, Surflex(Bound1,100), (Confort1,100) and (Confort6,300) lead to Best Pose with an average RMSD of 0.59, 0.68, and 0.54 Å, respectively. The superpositions of the Best Poses from Surflex(Confort6,300) versus the bound ligand structures from the corresponding complex crystals, shown in Figure 4a, are particularly encouraging. The superpositions of the Top Poses from Surflex(Confort1,100) and (Confort6,300) versus the bound ligand structures from the corresponding complex crystals are provided as Figure S6 and S7, respectively, where RMSDs are significantly large. Even though that the average RMSD obtained by Surflex(Confort1,100) for the Best Poses is not significantly higher than that by Surflex(Bound1,100), i.e., 0.68 vs 0.59 Å, some detailed analyses of the respective conformation space show the difference. As shown in Figure 5a, the proportion of “Good Poses” with RMSD < 1.0 Å can be quite low from

Results and Discussion below). Hence, the DOX protocol goes beyond the rigid receptor approximation commonly used in structure-based drug design where the protein structure is held fixed. To allow tryout of the DOX protocol, a web server is provided: http://202.114.32.71:10280/wandox/DOXserver/ home.html.

3. RESULTS AND DISCUSSION 3.1. Coarse Level: Surflex(Confort6,300) in Replacement of (Bound1,100). Abbreviated as “D” in the DOX protocol, binding conformation space search employing molecular docking is the first stage of the DOX calculations. Table 1 summarizes the performance of different docking protocols in the redocking tests on the 15 HMGR-statin complexes. Here the “Top Pose” refers to the dockingpredicted optimal binding pose with top score, while RMSD refers to the corresponding difference between the dockingpredicted binding pose and the bound ligand geometry in the crystal structure of the binding complex. The “Best Pose” refers to the docking-predicted binding pose which is the most similar to the corresponding ligand structure in the crystal. It is immediately clear from Table 1 that the Top Pose is usually not the Best Pose or the Best Pose usually does not have the highest score. Such an observation indicates that the scoring function used in our docking protocol may still have room for further improvement. 4269


Article


Figure 5. Proportion of Good Poses (RMSD < 1.0 Å) in (a) docking generated binding pose ensemble and (b) medium-level theory credited top 10 pose ensemble.

Surflex(Confort1,100) calculations for some systems. This is particularly true for 1HW8, 1HW9, 1HWL, 2Q6C, and 2R4F, while even zero occurrence of the Good Pose happens at 2Q6B. Undoubtedly, Surflex(Bound1,100) can always generate some Good Poses because of its satisfactory initial input, highlighting the importance of the quality of the initial input at the coarse level. Figure 5a also summarizes the numbers of the Good Poses from Surflex(Confort6,300). Encouragingly, Surflex(Confort6,300) are always able to capture more Good Poses without relying on the bound ligand structure. Six different initials can naturally have more chance, predicting better and more Good Poses in a larger pool (i.e., 300 vs 100 conformations). This provides a guarantee that the goal of the coarse level can be better fulfilled, which, in turn, ensures a better chance for the success of the subsequent procedures in the DOX protocol. 3.2. Medium Level: PM7/OPT in Replacement of ONIOM3/SP. The goal of the medium level is to taper the hundreds of docking poses to a couple of top candidates for final evaluation at the fine level. Hence, the CSAMP strategy demands the computational cost of the medium-level theory to be low enough to go through the conformation space from the coarse-level theory efficiently, while its accuracy is high enough to make sure the right one is still included in the 10 top candidates at the medium level.

The testing results up to the medium level on the 15 HMGR-statin complexes are summarized in Table 2 and Figure 5b. We present the RMSDs for the Top Poses and the Best Poses as well as the numbers of the Good Poses. Clearly Top Poses are still of no concern, as it is not the purpose for the medium level of theories to propose the Top Poses. In fact, as shown in Table 2, all mean RMSDs for the Top Poses overshoot the RMSD < 1.0 Å standard, while the individual RMSD for a specific ligand can be more than 2.4 Å even if it starts from Surflex(Bound1,100). On the other hand, RMSDs for the Best Poses are more relevant. In particular, the numbers of the Good Poses are of the major concern. In the previous DOX protocol, the medium level corresponds to ONIOM3/SP. As compared to the results on top of Surflex(Bound1,100), those based on Surflex(Confort1,100) are clearly less satisfactory (see Table 2 and Figure 5b). The latter leads to an increased mean RMSD to 0.72 Å for the Best Poses, while the corresponding value is 0.62 Å for the former. In particular, in 3 out of 15 cases (i.e., 1HW8, 1HW9, and 2Q6B), even the Best Poses exceed our standard of RMSD < 1.0 Å if the starting point is from Surflex(Confort1,100). Figure 5b shows a detailed analysis of the top 10 binding poses credited by the medium level of theory, which further reveals the barely acceptable performance of the previous version of the DOX protocol if the bound ligand 4270


Article

Journal of Chemical Theory and Computation Table 2. RMSD of the Medium Theories Credited Binding Poses: Redocking Test on 15 HMGR-Statin Complexes Bound+ONIOM3/SPa

Bound+PM7/OPTb

Confort1+ONIOM3/SPc

Confort1+PM7/OPTd

Confort6+PM7/OPTe

PDB ID

Top Posef

Best Poseg

Goodh

Top Pose

Best Pose

Good

Top Pose

Best Pose

Good

Top Pose

Best Pose

Good

Top Pose

Best Pose

Good


1.23 1.66 0.64 0.48 1.59 0.45 0.75 0.58 0.56 1.55 0.97 0.75 2.42 1.70 0.57 1.06

0.90 0.35 0.64 0.48 0.40 0.45 0.60 0.58 0.47 0.79 0.84 0.73 0.83 0.63 0.57 0.62

2 6 5 5 6 8 6 2 7 1 2 5 1 6 3 4.3

0.80 0.64 4.66 4.59 1.01 0.34 0.58 0.74 0.36 0.87 0.97 2.44 2.37 0.77 0.94 1.47

0.80 0.44 0.41 0.35 0.38 0.34 0.31 0.48 0.36 0.66 0.49 0.54 0.56 0.37 0.48 0.46

5 7 3 9 4 6 8 7 3 5 5 2 1 7 5 5.1

2.67 1.38 0.80 0.91 1.19 0.95 1.23 1.24 0.76 2.31 0.55 0.75 0.43 1.71 2.19 1.27

1.01 1.14 0.58 0.58 0.59 0.95 1.01 0.63 0.76 0.75 0.55 0.64 0.43 0.56 0.58 0.72

0 0 4 5 3 1 0 4 2 2 1 5 3 8 6 2.9

0.94 0.60 4.65 4.62 0.52 0.98 5.79 0.47 1.12 0.72 0.94 1.53 0.89 0.93 0.97 1.71

0.90 0.60 0.43 0.46 0.52 0.98 0.95 0.47 0.53 0.46 0.64 0.49 0.49 0.35 0.54 0.59

4 2 3 5 4 2 1 4 2 8 4 3 3 7 5 3.8

0.71 0.61 4.64 4.61 0.88 0.34 0.67 0.46 1.05 0.85 0.92 1.91 1.02 1.17 0.97 1.39

0.71 0.54 0.41 0.51 0.42 0.34 0.51 0.46 0.37 0.65 0.43 0.52 0.99 0.62 0.55 0.54

3 2 5 4 6 6 9 6 2 8 5 1 1 7 8 4.9

a

Suflex(Bound1,100) generated binding poses screened by ONIOM3 single-point calculation. bSuflex(Bound1,100) generated binding poses screened by PM7 optimization calculation. cSuflex(Confort1,100) generated binding poses screened by ONIOM3 single-point calculation. d Suflex(Confort1,100) generated binding poses screened by PM7 optimization calculation. eSuflex(Confort6,300) generated binding poses screened by PM7 optimization calculation. fThe binding pose with the lowest medium theory evaluated total energy. gThe binding pose most similar to the crystal binding structure among the top 10 poses with the lowest medium theory evaluated total energy. hNumber of binding poses with RMSD less than 1.0 Å in the top 10 poses.

structure is not allowed to be used. The Surflex(Confort1,100) + ONIOM3/SP combination provides no Good Pose at all in cases of 1HW8, 1HW9, and 2Q6B (see Table 2). In the new DOX protocol, PM7/OPT replaces ONIOM3/ SP as the medium level. The merits or demerits are clearly shown in Table 2. Again, Top Poses are of little concern. Regardless of the coarse-level results, the PM7/OPT-credited Top Poses can have a RMSD which is more than 4.5 Å. This is a reflection that, in terms of the energy accuracy per se, PM7 is inferior to ONIOM3(ωB97X-D/6-31G*:PM6:Amber). On the other hand, structure optimization is very beneficial to locate the Best Poses. When combined with Surflex(Bound1,100), PM7/OPT decreases the average RMSD for the Best Poses to only 0.46 Å, while when combined with Surflex(Confort1,100) and Surflex(Confort6,300), PM7/OPT yields the average RMSDs of only 0.59 and 0.54 Å for the Best Poses. In particular, none of the Best Poses credited by PM7/ OPT exceeds our standard of RMSD < 1.0 Å. Figure 5b presents a detailed analysis of the top 10 binding poses credited by the medium level of PM7/OPT, which further demonstrates the importance to invoke structure optimizations. As compared to the results from the corresponding ONIOM3/SP combinations, substitution of PM7/OPT undoubtedly increases the numbers of Good Poses. For example, zero Good Pose is credited for 2Q6B with the Surflex(Confort1,100)+ONIOM3/SP combination, whereas the PM7/OPT combination succeeds in improving one pose to a “good” one and put it in the top 10 (See Figure S8 in the Supporting Information for some details). Also zero Good Pose is credited for 1HW8 and 1HW9 with the Surflex(Confort1,100)+ONIOM3/SP combination, although the previous stage Surflex(Confort1,100) does capture two Good Poses (Figure 5a and Table 1). With the help of PM7/OPT, four Good Poses are credited for 1HW8 and two Good Poses are credited for 1HW9 (Figure 5b and Table 2).

The average performance of Surflex(Confort6,300)+PM7/ OPT is better than Surflex(Confort1,100)+PM7/OPT. The combination of Surflex(Confort6,300)+ONIOM3/SP is not considered in the present work because PM7/OPT clearly outperforms ONIOM3/SP and also we do not intend to spend too much computation resource at the medium level. Figure 6 shows the computational cost of each “coarse + medium” combination. In fact, the Surflex(Bound1,100) or Surflex(Confort1,100)+PM7/OPT combinations are even less costly than the corresponding ONIOM3/SP combination. Certainly,

Figure 6. Average computational cost of coarse + medium combination for the process of one HMGR-statin complex. Calculations were performed on E5-2683 dual-processor machines, and computational cost is expressed in core × hours (32 cores × real time). 4271


Article

Journal of Chemical Theory and Computation ONIOM3/OPT is superior to PM7/OPT, but is too expensive to be the medium-level theory. It is true that the computation cost is linearly increased when the Surflex(Confort6,300)+PM7/OPT combination is adopted, as the sample size is then three times larger. However, these calculations only account for ∼2% of the total amount of computer time as compared to the fine level of XO2 calculations. In fact, the combination of Surflex(Confort6,300)+PM7/OPT is our recommendation in terms of accuracy and efficiency to better fulfill the goals of the coarse and medium levels in the CSAMP strategy. 3.3. Performance of the New DOX Protocol: Redocking Tests on 28 HMGR, SD, NA, and 3HNR Complexes. In this work, a new version of the DOX protocol is introduced, with previous coarse and medium levels of theories being replaced by Surflex(Confort6,300) and PM7/OPT. As before, the fine level of the DOX protocol invokes the XO2 calculations, where each of the top 10 structures credited by medium-level theory is optimized using XO2 (ωB97X-D/631G*:PM6), and the final Top Pose is determined according to the order of the RBEs calculated using XO2(ωB97X-D/6311+G**:PM6). The new DOX protocol was first tested on the same set of 15 HMGR systems employed in our previous work and the other 13 complexes from SD, NA, and 3HNR systems. The final DOX-predicted binding structures for these 28 complexes are shown in Figure 4b, which are superimposed on the bound ligand structures in the corresponding crystal structures. With the help of the fine level of XO2 calculations, the overall performance of the DOX prediction is satisfactory. The average RMSD for 15 HMGR-statin complexes is only 0.50 Å, which is just one-half of the average RMSD from molecular docking top score predictions (i.e., 1.04 Å obtained by Surflex(Confort6,300)), one-half of our Good Pose standard (i.e., RMSD < 1.0 Å),12 and one-quarter of the standard (i.e., RMSD < 2.0 Å) usually employed in the literature.10 It is also slightly better than the average RMSD of 0.54 Å obtained in our previous work.12 The present work starts from the free ligand redocking, while the previous work starts from the bound ligand redocking. Hence, the new DOX protocol is more robust. Just as that in the HMGR-statin cases, the average RMSD for the DOX-predicted SD, NA, and 3HNR structures is only one-half of the 1.0 Å standard, which is a significant improvement as compared to the docking credited geometries. Figure 4b also displays the RMSD for each of the 28 complexes. In none of the cases, there is a RMSD that exceeds the 1.0 Å standard. The largest RMSD is 0.89 Å, which occurs in 3CD0, where there is a rotatable aromatic ring which is considerably flexible. All in all, the DOX-predicted binding mode is essentially the same as that appearing in the corresponding crystal structure of the binding complex. Note that the fine level of the DOX protocols corresponds to the XO2 SP//OPT calculations, where the basic requirement mainly relies on its accuracy. Hence, we further verify the accuracy of the XO2-calculated RBEs, comparing those on top of the crystal binding structures, as well as those on top of the Top Poses from Surflex(Confort6,300) and PM7/OPT. As the binding pose credited by DOX is taken as the reference, the corresponding RBE for each binding complex is 0.00 kcal/mol by definition as shown in Table 3. While a positive RBE with respect to the reference indicates that the referred binding pose is less stable than the DOX prediction, a

Table 3. Results of New DOX Version Free Ligand Redocking Test on HMGR, SD, NA, and 3HNR Complexes: Relative Binding Energy (RBE)a PDB ID HMGR 1HW8 1HW9 1HWI 1HWJ 1HWK 1HWL 2Q6B 2Q6C 2R4F 3BGL 3CCW 3CCZ 3CD0 3CD5 3CD7 mean absf SD 3STD 4STD 5STD 6STD 7STD NA 2HTU 2HTQ 1F8B 4KS1 3HNR 1YBV 1DOH 1G0N 1G0O mean absg

DOX RBEa

DOX RMSDb

crystal RBEc

docking RBEd

PM7/OPT RBEe

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.72 0.67 0.54 0.58 0.41 0.32 0.25 0.39 0.35 0.50 0.36 0.51 0.89 0.49 0.47 0.50

−0.53 −2.45 1.26 −1.80 1.40 0.66 −0.18 0.66 −0.67 0.15 0.58 0.84 0.02 −2.09 −0.55 0.92

25.89 19.69 7.99 4.27 2.23 6.08 21.59 12.52 0.72 12.78 23.47 0.87 1.10 6.80 8.06 10.27

0.00 0.00 10.51 12.05 10.85 0.00 5.80 0.00 1.92 1.89 4.40 15.83 6.74 12.32 2.77 5.67

0.00 0.00 0.00 0.00 0.00

0.56 0.47 0.75 0.29 0.58

0.67 −1.42 0.23 −0.89 −0.08

5.50 4.07 13.28 −2.57 0.79

0.00 2.83 1.25 0.00 0.00

0.00 0.00 0.00 0.00

0.34 0.71 0.79 0.57

0.90 −1.91 −1.55 0.49

4.79 9.86 5.30 11.62

2.02 20.53 0.11 13.95

0.00 0.00 0.00 0.00

0.28 0.24 0.68 0.21 0.50

0.13 0.21 −0.59 0.43 0.83

−0.23 5.05 0.54 1.15 7.81

0.00 7.49 1.02 0.00 4.80

a

Slightly different from our previous work,12 the reference for RBE calculation is always taken as the DOX-predicted top-ranking binding pose in this work. Hence, the RBE of all DOX-predicted top-ranking binding poses are naturally 0. We still put these zeroes here for an emphasis of RBE definition. bRMSD of DOX-predicted top-ranking binding pose. cCrystal binding structure, optimized with the XO2 method for the convenience of energy comparison. dSurflex(Confort6,300) docking procedure-predicted top score binding pose, optimized with the XO2 method for the convenience of energy comparison. ePM7/OPT-predicted optimal binding pose, optimized with the XO2 method for the convenience of energy comparison. f Mean absolute value of 15 HMGR systems for comparison with our previous work. gMean absolute value of all 28 systems.

negative RBE indicates that the referred binding pose is more stable than the DOX prediction. The RBE results credited by the XO2 are summarized in Table 3. The mean absolute RBEs for the crystal structures of these 28 systems is 0.83 kcal/mol, which suggests that the accuracy of XO2 has achieved a chemical accuracy of 1.0 kcal/mol. In most of the cases (26 out of 28), the calculated RBEs of crystal binding structures are within the range of ±2.00 kcal/mol, indicating that the DOXpredicted binding poses are usually energetically closed to the 4272


Article


Table 4. RMSD of Docking, PM7/OPT, and DOX-Predicted Top Pose in Free Ligand Redocking Tests on Astex Diverse Set (Angstroms) PDB ID

dockinga

PM7/OPTb

DOXc

PDB ID

dockinga

PM7/OPTb

DOXc

1G9V 1GKC 1GM8 1GPK 1HNN 1HP0 1HQ2 1HVY 1HWId 1HWW 1IA1 1IG3 1J3J 1JD0 1JJE 1JLA 1K3U 1KE5 1KZK 1L2S 1L7F 1LPZ 1LRH 1M2Z 1MEH 1MMV 1MZC 1N1M 1N2J 1N2V 1N46 1NAV 1FO1 1OF6 1OPK 1OQ5 1OWE 1OYT 1P2Y 1P62 1PMN 1Q1G 1Q41 1Q4G 1R1H 1R55

0.97 1.39 2.21 4.32 5.49 2.08 0.17 2.79 0.89 1.94 0.82 1.48 0.53 5.87 8.17 1.02 0.60 1.15 0.97 1.54 0.70 1.84 2.06 0.32 2.55 2.19 1.54 1.04 2.48 0.79 0.60 1.00 0.72 1.30 0.90 1.23 2.21 1.02 0.72 0.61 1.45 0.83 0.53 0.88 0.79 0.81

1.45 1.13 1.63 0.43 2.13 0.47 0.51 3.72 4.64 0.24 1.09 3.75 0.53 2.74 0.37 1.02 0.90 0.50 0.97 2.52 1.13 1.32 1.11 0.38 2.05 1.14 0.59 4.66 0.39 1.62 0.63 0.33 0.69 0.37 0.59 4.27 2.26 1.62 0.55 0.55 1.44 0.78 0.54 0.30 1.20 0.39

1.50 0.86 1.70 0.21 0.93 0.34 0.57 1.09 0.65 0.27 1.13 0.66 0.46 2.25 0.42 0.64 0.86 0.48 0.44 1.58 0.37 1.02 0.25 0.52 0.49 0.32 0.51 1.19 0.35 0.55 0.43 0.75 0.60 0.39 0.35 1.01 1.83 0.43 0.38 0.42 0.56 0.93 0.91 0.26 1.09 0.67

1R58 1R9O 1S19 1S3V 1SG0 1SJ0 1SQ5 1SQN 1T40 1T46 1T9B 1TOW 1TT1 1TZ8 1U1C 1U4D 1UML 1UNL 1UOU 1V0P 1V48 1V4S 1VCJ 1W1P 1W2G 1X8X 1XM6 1XOQ 1XOZ 1Y6B 1YGC 1YQY 1YV3 1YVF 1YWR 1Z95 2BM2 2BR1 2BSM mean

2.84 0.53 0.69 0.81 0.99 2.01 6.26 0.52 1.51 1.01 2.10 0.86 3.47 1.13 1.59 0.73 7.06 0.57 0.38 0.58 0.55 0.75 0.74 2.94 1.15 1.04 2.07 2.61 0.63 7.60 2.57 0.66 0.64 6.31 5.84 0.62 0.68 6.68 0.65 1.83

3.31 0.72 0.46 6.80 7.43 0.96 2.32 0.40 3.52 0.55 5.01 0.78 4.19 0.66 1.36 5.08 7.06 0.65 0.51 1.94 0.98 0.71 0.75 2.94 0.94 0.58 1.32 0.52 0.19 0.41 2.31 0.51 0.34 0.89 5.84 0.47 0.73 0.46 1.12 1.61

0.42 0.50 0.54 0.65 0.73 1.44 1.36 0.65 0.26 0.70 1.20 1.19 0.39 0.71 1.70 0.82 1.16 0.67 0.58 0.64 1.05 0.63 1.11 0.64 1.00 0.67 1.23 0.48 0.38 0.39 1.00 0.83 0.58 0.51 0.80 0.28 1.07 0.80 1.34 0.76

a Surflex(Confort6,300) docking protocol-predicted optimal top score binding pose. bSuflex(Confort6,300)-generated binding poses screened by PM7 optimization calculation; the pose with the lowest PM7 energy. cThe Top Pose credited by XO2 opt/energy. d1HWI is also included in the 15 HMGR systems discussed above.

flexible hydrophobic functional group. There are four cases (1HWJ, 4STD, 2HTQ, and 1F8B) where RBEs of crystal binding structures are 1.4−1.9 kcal/mol more stable, indicating that there is still room to increase the number of Good Poses by improving the coarse and medium levels. On the other hand, the average RBE for Top Poses credited by molecular docking is 7.81 kcal/mol, indicating that the docking-predicted “optimal” binding structures are on average ∼8 kcal/mol less stable than the DOX-predicted structures for these 28 complexes. It has to be noted that in some worst cases (8 out of 28) the docking credited Top Poses could be 10−20 kcal/mol higher in energy than the DOX-predicted structures

corresponding crystal structures. This observation echoes that they are also structurally similar (see Figure 4b and also Table S9 in the Supporting Information for some more details). The largest RBE is found to be −2.45 kcal/mol in the case of 1HW9, indicating that the crystal structure is more stable. Hence, if a similar pose to the crystal structure has been encompassed in the sampling space, XO2 calculation would have credited it as the top 1 pose. Most likely the large and flexible hydrophobic functional group in the 1HW9 ligand has brought about challenge to the molecular docking used for the initial conformation space generation. The second largest RBE is −2.09 kcal/mol for 3CD5, which also possesses a large and 4273


Article


rate is 72%, suggesting that in most of the cases the DOXpredicted top binding pose is comparable with the crystal structure. Note that the popular docking tools like SurflexDock@SYBYL can only achieve a RMSD < 2 Å success rate of 70% and a RMSD < 1 Å success rate of 46% (by Top Poses). However, Surflex docking Best Pose can achieve a RMSD < 2 Å success rate of 100% and a RMSD < 1 Å success rate of 90%, which shows that the sampling space of binding conformations generated by Surflex is of high quality that makes a solid foundation for the good performance of DOX. As compared to other popular docking tools, GOLD can achieve a RMSD < 2 Å success rate of 75% (by Top Poses) on the same Astex set,21 employing a similar free ligand redocking protocol. Glide36−38 can achieve a RMSD < 2 Å success rate of 82% and a RMSD < 1 Å success rate of 66% also on the Astex set39 (by Top Poses). We also carried out a test on the Astex set using Autodock Vina,40 with which the RMSD < 2 Å success rate is 72%, while the RMSD < 1 Å success rate is 42% (by Top Poses, see Supporting Information Table S54). All of them are inferior to the DOX protocol. The only case in which the DOX-predicted Top Pose does not meet the RMSD < 2 Å standard is 1JD0. We found that nearly one-half of the small ligands are surrounded by crystallized water molecules in 1JD0. Thus, its binding structure is significantly affected by the water molecules bridging the ligand and the protein which are, nonetheless, deleted in the DOX calculations. Such a water problem, which also exists in the SD systems mentioned in the previous section, shall be overcome in our future study. Another case worthy of mention is 1YVF, which shows the importance of proper treatment of ligand intramolecular interactions. The ligand in 1YVF possess a π-conjugated rigid backbone, including the carboxylic acid. Therefore, the ligand intramolecular interactions are critical to the ligand structure, which are properly addressed by introducing QM optimizations in the middle and fine levels of the DOX protocol. It explains why DOX has exhibited good performance in this case (RMSD = 0.51 Å). On the contrary, with poor consideration of ligand intramolecular interactions at the MM level, the RMSD of the Top Pose as credited by Surflex scoring function is as high as 6.31 Å. The overlay of the docking prediction with the crystal structure in Figure S9 shows that the π-conjugated interaction is broken. Autodock vina also exhibits a poor performance in this case (RMSD = 2.1 Å, see Table S54). Note that the solvation effect is not considered in the current DOX protocol mainly owing to the difficulty to implement an accurate and affordable solvation model. On the other hand, luckily, we found that the solvation effect is less important for binding pose ranking in our present cases, see Table S7 in the Supporting Information. This partly supports the generally good performance of the current DOX protocol as well as the reasonable performances of the docking tools like Glide, Autodock vina, and Surflx against the Astex set. We also notice that in about 20% of the cases the DOX protocol failed to pick out the docking Best Poses, likely owing to the lack of solvation. We will try to address this problem in our future work. 3.5. Performance of the New DOX Protocol: CrossDocking Test on 15 HMGR-Statin Complexes. In real condition of drug design, the crystallographic data are often not available for the binding complex of interest. Hence, the desired ligand would have to be docked to a protein structure

or the crystal structures. This observation downplays the role of using Top Poses credited by molecular docking in the structure-based drug design. Table 3 shows that the PM7/ OPT credited Top Poses are significantly better with an average RBE of 4.80 kcal/mol. There are some favorable cases (9 out 28) that the PM7/OPT-credited Top Poses are identical to the DOX predictions in energy. Still, there are some unfavorable cases that the PM7/OPT-credited Top Poses are far from the predictions of the DOX or the crystal structures. For example, RBEs of Top Poses credited by PM7/ OPT overshoot by 12.05 kcal/mol for 1HWJ, 15.83 kcal/mol for 3CCZ, 12.32 kcal/mol for 3CD5, 13.95 kcal/mol for 4KS1, and 20.53 kcal/mol for 2HTQ. In addition, we point out that, in all of the SD-ligand complexes, one water molecule is involved which forms a H bond with both the ligand and the protein, acting as a bridge between the two (see Figure S9 in the Supporting Information for a graphic illustration). However, as a standard step suggested in the literature,10 all water molecules are removed at the beginning of the molecular docking procedure. Hence, the bridge water molecule in the SD systems is absent in all steps of the DOX calculations which more or less leads to a worsened prediction. Even though such a problem has not affected the results severely, as indicated by the corresponding RMSDs and RBEs, it is a problem which shall be addressed in a future study. 3.4. Large-Scale Test of the New DOX Protocol: Redocking Tests on the Astex Set. To make further assessment and comparison with traditional molecular docking methods, we performed free ligand redocking tests on the Astex set,21 which is a diverse set containing 85 protein−ligand complexes that cover various pharmaceutical or agrochemical targets and ligand molecules ranging from 10 to 60 heavy atoms. The results are summarized in Table 4, which includes the RMSD of Top Poses as credited by the coarse level, medium level, and fine level of theories. The coarse level of Surflex-Dock predicted the Top Poses to have an average RMSD of 1.83 Å. Seventy percent of the Top Poses meet the loose yet commonly referred to standard for RMSD < 2 Å. We noticed that Jain reported the same 70% success rate for Surflex-Dock on the same Astex set,35 employing a similar free ligand redocking protocol. If the success standard is tightened to RMSD < 1 Å, the success rate declines to 46% (50% in Jain’s work34). There are some worst cases (9 of 85), where the RMSDs for the docking credited Top Poses exceed 5 Å, being of no reference value at all. On the other hand, the docking-credited Best Poses (Table S53) represent an average RMSD of 0.60 Å, with no individual RMSD that is larger than 2 Å. This observation once again indicates that docking is able to generate good poses, although it often fails to credit the Top Poses correctly. The medium level of PM7/OPT generally leads to a better, albeit not by much, prediction of the Top Poses than the coarse level. The average RMSD for the Top Poses is 1.61 Å, and the RMSD < 2 Å success rate is 75%. On the other hand, the Top Ten Best Poses (Table S53) have an average RMSD of 0.67 Å, with no individual RMSD that is larger than 2 Å. This observation once again suggests that the medium level has also achieved its goal in the DOX protocol. The most encouraging result has been obtained at the final step of the DOX protocol. The Top Poses credited by the fine level have an average RMSD of 0.76 Å. The RMSD < 2 Å success rate is as high as 99%, while the RMSD < 1 Å success 4274


Article

Journal of Chemical Theory and Computation extracted from the binding complex with some other ligand. In this work we explore the performance of the DOX protocol in such situations by using cross-docking tests in which 14 statin molecules (except the ligand in 3CCZ) are docked to the protein structure extracted from 3CCZ, while the 3CCZ ligand is docked to the 2R4F protein. The Surflex(Confort6,300) results are summarized in Table 5. As one could expected, the docking performance degrades as Table 5. Results obtained in HMGR-statin complexes crossdocking testsa: Surflex molecular docking. Unit: Å Top Poseb d

Best Posec e

PDB ID

score

RMSD (Å)

score

RMSD

Goodf


12.98 12.75 12.65 12.82 15.92 12.30 14.69 13.77 15.06 14.45 13.94 14.49 14.66 15.08 16.06 14.11

2.48 3.12 1.63 0.84 0.61 0.65 0.90 0.85 0.51 1.35 1.37 0.72 0.81 0.70 0.65 1.15

11.68 11.53 11.03 10.94 15.48 10.21 13.15 11.89 15.06 12.09 13.70 12.80 13.27 13.67 14.54 12.73

1.00 0.93 0.63 0.44 0.47 0.63 0.54 0.69 0.51 0.63 0.49 0.54 0.52 0.42 0.50 0.60

0 8 40 47 97 21 86 32 28 82 18 61 34 150 85 52.6

a

All of the 14 statin molecules in the listed crystal structure were docked to the protein extracted from the 3CCZ crystal structure employing the Confort(6,300) procedure, while the ligand of 3CCZ was docked to the protein extracted from the 2R4F crystal structure. Ligand RMSD’s were calculated after protein−protein superimposition according to coordinates of all backbone heavy atoms. b Surflex(Confort6,300) molecular docking-predicted optimal top score binding pose. cThe binding pose most similar to the crystal binding structure. dBinding affinity of the corresponding binding pose scored by empirical scoring function. eRMSD between the predicted binding pose and the bound ligand geometry coordinates in complex crystal structure. fNumber of binding poses with RMSD less than 1.0 Å generated in corresponding docking calculations.

Figure 7. Superposition of cross-docking-predicted ligand binding structure and native ligand structure in the corresponding complex crystal, obtained by protein−protein superimposition according to coordinates of all backbone heavy atoms. (a) Surflex(Confort6,300)predicted top score binding structure (yellow) versus crystal (gray). (b) DOX prediction (cyan) versus crystal (gray).

opening the possibility of further improvement in the later DOX procedures. Table 6 and Figure 7b demonstrates the promising performance of the DOX protocol in the cross-docking tests. The Top Poses credited by DOX processes averaged a RMSD of 0.63 Å with none of the cases where the individual RMSD exceeds the 1.00 Å standard as a Good Pose. The DOX protocol goes beyond the rigid receptor approximation, while the protein structure of the important amino acid residues surrounding the ligand is optimized at the last stage of the DOX protocol via XO2/OPT. We believe that it is important that the flexibilities of the protein−ligand binding poses can be properly taken into account in the structure-based drug design where the desired crystallographic data for the binding structure is often unavailable. Table 6 also summarizes the RBEs calculated by the fine level of the XO2 method based on the Top Poses from the coarse and medium levels, where the corresponding RBE from DOX is taken as the reference. In fact, all of the Top Poses from the low levels are inferior to the fine level. The mean RBEs are 7.77 and 9.42 kcal/mol for the coarse and medium

a result of the rigid receptor approximation. The average RMSD for the Top Poses increases from the original 1.04 Å obtained in the redocking tests (see Table 1 and Figure 7a) to 1.15 Å in the cross-docking tests. Comparing with each individual case in the HMGR-statin set, the most significantly worsened predictions in the cross-docking tests occur at 1HW8 (RMSD 2.48 Å), 1HW9 (3.12 Å), and 1HWI (1.63 Å), most likely because these ligand molecules are quite different from the original ligand molecule of 3CCZ. Under the rigid receptor approximation, the protein residues in 3CCZ are not allowed to relax in the docking calculations; hence, they are not perfectly positioned to accommodate the other ligand molecules. Of course, the main concern, according to the CSAMP strategy, is the Best Poses and the number of the Good Poses rather than the Top Poses. The average RMSD of the Best Poses from the cross-docking tests is 0.60 Å (see Table 4), which can be compared to the average RMSD of 0.54 Å from the redocking tests (see Table 1). Still the number of Good Poses for each system is sufficiently good (see Table 5), 4275


Article

Journal of Chemical Theory and Computation Table 6. Results Obtained in HMGR-Statin Complexes Cross-Docking Testsa: DOX PDB ID

DOX RMSDb

docking RBEd

PM7/OPT RBE

1HW8 1HW9 1HWI 1HWJ 1HWK 1HWL 2Q6B 2Q6C 2R4F 3BGL 3CCW 3CCZ 3CD0 3CD5 3CD7 mean abs

0.77 0.57 0.65 0.62 0.43 0.39 0.81 0.73 0.43 0.79 0.70 0.73 0.94 0.45 0.37 0.63

27.66 23.04 2.60 4.14 1.80 3.91 1.81 13.37 1.27 15.40 10.89 3.52 5.13 2.04 0.04 7.77

11.28 0.91 0.87 4.05 12.60 27.81 0.00 15.27 18.09 5.54 8.54 14.09 12.24 3.08 6.87 9.42

Table 7. Ranking of Best Pose According to Docking Score of SYBYL, PM7 Energy, and XO2 Energy: Redocking Test on HMGR, SD, NA, and 3HNR Complexes PDB ID HMGR 1HW8 1HW9 1HWI 1HWJ 1HWK 1HWL 2Q6B 2Q6C 2R4F 3BGL 3CCW 3CCZ 3CD0 3CD5 3CD7 SD 3STD 4STD 5STD 6STD 7STD NA 2HTU 2HTQ 1F8B 4KS1 3HNR 1YBV 1DOH 1G0N 1G0O

a

All of the 14 statin molecules in the listed crystal structure were docked to the protein extracted from the 3CCZ crystal structure employing the Confort(6,300) procedure, while the ligand of 3CCZ was docked to the protein extracted from the 2R4F crystal structure, at the first stage of DOX calculation. bRMSD between the DOXpredicted top binding pose and the bound ligand geometry coordinates in complex crystal structure. Ligand RMSD were calculated after protein−protein superimposition according to coordinates of all backbone heavy atoms. cSurflex(Confort6,300) docking procedure-predicted top score binding pose, optimized with the XO2 method for the convenience of energy comparison. dPM7/ OPT-predicted optimal binding pose, optimized with the XO2 method for the convenience of energy comparison.

levels, respectively. It once again suggests that the docking scoring function is unreliable to predict correctly the optimal binding pose and that high-level QM theory is necessary to address the challenge. 3.6. Discussion: When to use DOX? Molecular docking has become a powerful tool for structure-based drug design. A reported benchmark study on molecular docking10 is based on a large sample size (∼1000) for statistical significance. Indeed, the DOX protocol is much more computationally demanding. Hence, DOX is not designed to compete with molecular docking for large-scale applications, e.g., high-throughput virtual screening. On the other hand, DOX is much less empirical. The medium- and high-level theories implemented in DOX rely on QM with a limited number of parameters for general purpose. The molecular docking employed in DOX intends to encompass the Best Poses rather than to score the Top Pose. Data in Table 1, Figure 4, and Table 4 reveal that the Top Poses as credited by molecular docking can be very good or very bad on different systems, while the Best Poses are always predicted satisfactorily. The good performance of the DOX protocol is a joint contribution from all three levels of theory, each of which achieves its own goal based on the CSAMP strategy. Table 7 shows how the Best Poses are tracked down steadily from the coarse level to the medium level and finally are located as the Top Poses by the fine level of theory. In about one-half of the cases (14 out of 28), the Best Poses are ranked after 100 according to the docking scores, while only in 2 cases (i.e., 2R4F and 1G0O) the Best Poses are ranked as top 10. This suggests that for the docking to more safely encompass the

dockinga

PM7b

XO2c

167 137 225 24 18 20 51 198 2 63 130 29 15 179 19

1 7 5 2 3 1 6 1 4 6 5 7 5 6 5

2 2 2 1 1 1 1 1 1 1 1 1 1 1 1

149 221 145 169 204

2 3 8 1 2

2 1 3 1 2

74 245 69 116

1 3 8 2

1 1 1 2

16 36 103 5

1 7 5 1

1 2 1 1

a

The ranking of the binding pose most similar to the crystal binding structure among all of the docking generated poses, according to the docking score of SYBYL. bThe ranking of the binding pose (after optimization) most similar to the crystal binding structure among all of the docking generated poses (after optimization), according to the PM7/OPT-estimated total energy. cThe ranking of the binding pose most similar to the crystal binding structure among all of the XO2evaluated binding poses, according to the XO2-estimated total energy.

Best Poses, the conformation space is better to include the first 300 poses. When coming to the medium level of PM7/OPT, the funnel is significantly narrowed down. In all of the cases, the Best Poses are ranked within 10, although only in 7 out of the 28 cases the Best Poses are actually the Top Poses. The advantage of the fine level of the XO2 theory is consolidated by the fact that in 20 out of the 28 cases the Best Poses are credited as top 1, while in the rest of the 7 cases the Best Poses are credited as top 2. Only in one case (i.e., 5STD) the Best Pose is credited as top 3. For free ligand redocking tests, more details on the top 10 binding poses as credited by the DOX protocol can be found in Tables S10−S37 is the Supporting Information, where the binding score and the corresponding RMSD (Å) given by SYBYL, PM7, and XO2 are listed for comparison. It is interesting to see that none of the Best Poses as credited by SYBYL ends up to be the top 10 poses as credited by DOX. A technical detail one has to note is that, for the sake of saving 4276


Article


conformation space of a funnel-like structure is searched from the molecular docking with hundreds of candidates to PM7/OPT with a couple of candidates to XO2 with the final optimal binding mode. The CSAMP strategy is employed as in the DOX protocol, aiming to approach the optimal binding pose accurately yet efficiently. Benchmark tests have been carried out against crystallographic data of 112 protein−ligand complexes, covering various pharmaceutical or agrochemical targets and drug-like ligand molecules. DOX predictions starting from either redocking or cross-docking have been performed without relying on the bound ligand structure and/or the corresponding protein structure extracted from the specific protein−ligand complex. Satisfactory results have been obtained with an impressive RMSD < 2 Å success rate of 99% and RMSD < 1 Å success rate of 79%. The present study suggests that the DOX protocol opens the possibilities toward accurate prediction of the protein−ligand binding structures, which shall prove useful for lead optimization in structure-based drug design. In addition, it is anticipated that the CSAMP strategy and the DOX protocol with possible variations hold the promise in accurate conformation searches for other chemical and biological systems, involving protein/DNA/RNA-ligand bindings, protein−protein bindings, as well as molecular recognition and assembly, etc.

the computational costs, there is a screening process before the medium level (See Figure 3). The survival poses are actually the good poses having a pose−pose RMSD smaller than 0.3 Å as compared to the Best Pose but with higher binding scores. The only exception in this set of free ligand redocking test is 1HW8, where the docking-generated Best Pose (see Figure 4a) is significantly worse than the DOX-credited Top Pose which has to be optimized using a high-level theory. Similar details to those in Tables S10−S37 are provided as Tables S38−S52 in the Supporting Information for crossdocking tests. One case worthy of mention is 3CD0 in the cross-docking test (Table S50. The top 1 pose predicted by DOX has a RMSD of 0.94 Å, whereas the Best Pose which possesses a RMSD of 0.51 Å is actually ranked as the second. The energy difference between them is very small (0.73 kcal/ mol). We found that in all cases of cross-dockings when the Best Poses are not ranked as top 1 they are actually energetically very close to the corresponding top 1 pose. It is possible that the current DOX protocol is not yet accurate enough to discriminate such situations or both binding structures are all possible candidates. Hence, it is understandable that the fine level of the XO2 method may sometimes miss to rank the Best Pose as top 1 pose as shown in Table 6. The users of the DOX method are encouraged to look into not only the XO2-credited top 1 pose but also the other poses within a small range of RBEs (

Conformation Search Across Multiple-Level Potential-Energy Surfaces

Recommend Documents