Exploring the Stability of Ligand Binding Modes to Proteins by

26-1, Muraoka-Higashi 2-chome, Fujisawa, Kanagawa 251-8555, Japan. J. Chem. Inf. Model. , 2017, 57 (10), pp 2514–2522. DOI: 10.1021/acs.jcim.7b0...
2 downloads 16 Views 3MB Size
Download PDF
Article Cite This: J. Chem. Inf. Model. 2017, 57, 2514-2522

pubs.acs.org/jcim

Exploring the Stability of Ligand Binding Modes to Proteins by Molecular Dynamics Simulations: A Cross-docking Study Kai Liu† and Hironori Kokubo*,‡ †

Drug Discovery Chemistry Laboratories, CNS Drug Discovery Unit, and ‡Partnership Research Center, Takeda Pharmaceutical Company Limited, 26-1, Muraoka-Higashi 2-chome, Fujisawa, Kanagawa 251-8555, Japan S Supporting Information *

ABSTRACT: Docking has become an indispensable approach in drug discovery research to predict the binding mode of a ligand. One great challenge in docking is to efficiently refine the correct pose from various putative docking poses through scoring functions. We recently examined the stability of self-docking poses under molecular dynamics (MD) simulations and showed that equilibrium MD simulations have some capability to discriminate between correct and decoy poses. Here, we have extended our previous work to cross-docking studies for practical applications. Three target proteins (thrombin, heat shock protein 90alpha, and cyclin-dependent kinase 2) of pharmaceutical interest were selected. Three comparable poses (one correct pose and two decoys) for each ligand were then selected from the docking poses. To obtain the docking poses for the three target proteins, we used three different protocols, namely: normal docking, induced fit docking (IFD), and IFD against the homology model. Finally, five parallel MD equilibrium runs were performed on each pose for the statistical analysis. The results showed that the correct poses were generally more stable than the decoy poses under MD. The discrimination capability of MD depends on the strategy. The safest way was to judge a pose as being stable if any one run among five parallel runs was stable under MD. In this case, 95% of the correct poses were retained under MD, and about 25−44% of the decoys could be excluded by the simulations for all cases. On the other hand, if we judge a pose as being stable when any two or three runs were stable, with the risk of incorrectly excluding some correct poses, approximately 31−53% or 39−56% of the two decoys could be excluded by MD, respectively. Our results suggest that simple equilibrium simulations can serve as an effective filter to exclude decoy poses that cannot be distinguished by docking scores from the computationally expensive free-energy calculations.

■

INTRODUCTION Predicting the binding pose of a ligand to a protein is as important as evaluating its binding affinity in drug design. There are many docking software tools1 available, the accuracy of which has been gradually improved by introducing new techniques, such as induced fit docking (IFD).2,3 As a result, current docking software can generate correct binding poses in many cases if a reasonable template structure of a protein is used. However, as the software typically generates many putative candidate poses, the challenge of how to efficiently select a correct pose among these candidates remains. Because a correct pose should have a stronger binding affinity to a protein than a decoy pose, it is considered that in principle, a correct pose can be distinguished by calculating its binding affinity. The molecular mechanics/Poisson−Boltzmann surface area (MM/PBSA) and molecular mechanics/generalized Born surface area (MM/GBSA) methods4−6 are often used to predict binding affinity by treating the solvent effect implicitly without atomistic detail. However, as suggested by Lu et al.7 in a study of 392 high-resolution crystal structures, about 3.5 water molecules on average are involved in binding through bridging the hydrogen bonds between a protein and its ligand. © 2017 American Chemical Society

Thus, it may be difficult to predict the affinity with high accuracy by MM/PBSA or MM/GBSA in these cases. Recently, accurate binding free-energy calculations based on perturbation theory and/or thermodynamics integration are becoming more practical, and many applications have been reported.8−13 However, although these methods may be applicable for refining the correct poses from putative poses, their application to a large number of poses for each ligand is in practice computationally too heavy. We considered that if the simulation force field parameters are accurate, then a correct pose should be stable enough to hold the initial pose more steadily than decoy poses would, even in equilibrium simulation runs. Therefore, we had systematically investigated the stability of ligand binding poses by molecular dynamics (MD) simulation in our previous study.14 We found that about 94% of the native poses maintained their stability during the simulations, which suggests that MD simulations are accurate enough to judge most experimental binding poses as being properly stable. InterestReceived: July 7, 2017 Published: September 13, 2017 2514

DOI: 10.1021/acs.jcim.7b00412 J. Chem. Inf. Model. 2017, 57, 2514−2522

Article

Journal of Chemical Information and Modeling

Figure 1. Physicochemical properties of all the ligands for the three targets: HSP90a, heat shock protein 90-alpha; CDK2, cyclin-dependent kinase 2.

change of a ligand was less than 2.0 Å in any one among five parallel MD runs with different initial velocities. In this study, we further examined the criteria of any one (= at least one), at least two, at least three, at least four, and all five MD runs.

ingly, incorrect decoy poses were much less maintained, and 38−44% of the decoys could be excluded just by performing equilibrium MD simulations, although 56−62% of the decoys were stable. These results suggested that MD simulations may be used as a filter of docking poses. However, these results were based on self-docking, and it was not clear whether they could hold true in practical situations of drug discovery research. Therefore, in this study, we investigated the stability of the ligand binding pose through MD simulation for both correct and decoy poses, generated by cross-docking under typical practical situations. Three different protein systems were used to take different situations into consideration for cross-docking, namely: simple cross-docking,15−17 IFD, and IFD with a homology model. We found that the results obtained by practical cross-docking were essentially similar to those obtained by self-docking in our previous study. We also examined how much the discrimination power changes depending on the criteria. In the previous study, we judged a pose as being stable if the root-mean-square deviation (RMSD)

■

MATERIALS AND METHODS

Target Selection. In our previous work,14 we created a data set of 1039 Protein Data Bank (PDB) files, which consisted of drug-like ligands with pharmaceutical properties. Based on this data set, three protein targets (viz., thrombin, heat shock protein 90-alpha (HSP90a), and cyclin-dependent kinase 2 (CDK2)) were selected for our current work. These three targets are of pharmaceutical interest and are already used in marketed drugs. In addition, the available numbers of PDB files were beneficial for the following systematic study. We selected 27, 41, and 50 complex structures for thrombin, HSP90a, and CDK2, respectively, for the cross-docking studies. Generation of Initial Poses. In order to obtain reasonable initial poses for the simulation, several docking protocols from 2515

DOI: 10.1021/acs.jcim.7b00412 J. Chem. Inf. Model. 2017, 57, 2514−2522

Article

Journal of Chemical Information and Modeling Schrödinger software18 were employed. Because of the high similarities among the PDB files of thrombin, two templates (PDB codes 3DA9 and 3QX5) were randomly selected for Glide19 standard precision (SP) mode docking. For HSP90a, some ligand-induced structural changes are known;20 therefore, aside from normal Glide SP mode docking, IFD with the extended sampling protocol was also performed for comparison. In the case of CDK2, although its crystal structures are available, homology modeling was conducted on the single template of CDK1 (PDB code 4Y72) by using Prime module18 with default settings. We built this test set as a practical example of when the experimental structure is not available. Because the side-chain conformations are expected to be inaccurate in the homology model, we performed IFD on the basis of this homology model structure. We think that these varieties of cross-docking protocols (simple docking, IFD, and homology model-based IFD) suitably reflect various practical docking situations. The initial conformations of the ligands were picked up from the corresponding PDB files and the protonated states were carefully examined. The complexes of the same target were aligned for the RMSD calculation.16 All residues of proteins were aligned by the super command in the PyMOL (version 1.4.1) software.20 We used the docking templates as the reference structures for alignment. Among the various docking poses generated by the different protocols, three poses (pose0, pose1, and pose2) were defined for each ligand, as done in our previous study (pose0 the closest pose to the crystal structure with RMSD < 2.0 Å; pose1 the smallest g-score pose with RMSD > 2.0 Å from both pose0 and the native pose; and pose2 the smallest g-score pose with RMSD > 2.0 Å from the native pose, pose0, and pose1). Pose1 and pose2 can be considered as two reasonable and challenging decoys. Here, the reference of the native pose for the RMSD calculation was obtained from the crystal structure after the alignment. Molecular Dynamics Simulation Protocol. We briefly describe our MD simulation protocol, which is the same as that used in our previous work.14 Two rounds of minimization were performed for the complex structure, each comprising 500 steps of steepest descent followed by 500 steps of conjugate gradient minimization. In the first round of minimization, a weak constraint (1.0 kcal/(mol·Å2)) was assigned to all heavy atoms. In the second round, no restraints were employed in order to allow minimization of the entire system. The system was then gradually heated from 50 to 300 K during 0.2 ns, while the heavy atoms of the protein were restrained with a force constant of 0.5 kcal/(mol·Å2). After that, the system was further equilibrated for 0.5 ns at a constant pressure (1.0 bar) with the same weak constraint force. In the final production step, a 10 ns simulation in the NPT (constant number− pressure−temperature) ensemble was performed without any constraints. For the CDK2 system, the simulation time was additionally extended to 100 ns because we considered that the longer equilibrium time would be necessary in the case of homology modeling. It is important to consider the dependency on initial velocities when we judge the stability of the simulations.21 We thus performed five independent simulations with different initial velocities. For the MD trajectory analysis, only the 500 snapshots extracted at even intervals from the last 1 ns trajectory of the MD simulation were used. During the RMSD calculation for a ligand, the protein backbone was initially aligned to the reference structure, which came from the

docking pose, and then the RMSD of heavy atoms of the ligand were calculated.

■

RESULTS AND DISCUSSION Basic Feature of the Ligands. Figure 1 shows the physicochemical properties of all the ligands. The general distribution of these properties of each target indicates the diversity of the ligands. Across the different targets, the distributions of ligand properties also varied. For example, the ligands of thrombin were relatively larger and more flexible than those of HSP90a. However, the estimated binding affinity, which was roughly calculated by −log10 K (K = Kd, Ki, or IC50)), indicated a weaker binding affinity to thrombin. The detailed properties of each ligand are listed in Table S1 of the Supporting Information (SI). Docking of Ligands from Thrombin. The binding pockets of 27 thrombin PDB files were initially aligned and compared with one another. The average RMSD value of core residues (residues within 6 Å of the ligand) was only 0.3 Å, which indicated the conservative movement of the pocket (shown in Figure S1 of the SI). However, the low average Tanimoto coefficient of the ligands (0.3) suggested their diversity (shown in Figure S2 of the SI). Here, we randomly selected two templates (PDB codes 3DA9 and 3QX5) for the Glide SP docking. Figure 2 shows the docking capability (or docking power)22 by the two different templates. Here, pose0 is the “best pose”

Figure 2. Root mean square deviation (RMSD) values of pose0 generated by docking with two different templates.

according to the definition in Materials and Methods. If the RMSD value between the best docking pose and the crystallographic one was larger than 2.0 Å, then we regarded the docking of that ligand as a failed case. As a result, only three and five ligands were failed for the docking based on templates 3DA9 and 3QX5, respectively. Therefore, the docking capability of Glide SP mode was quite good, and the dependence of the docking on the template was not obvious in the current test. At the same time, the g-score function also presented the high accuracy. Among the 46 cases successfully docked by the two templates, about 98% (45 ligands) were properly ranked in the top 10, and 43% (20 ligands) were correctly ranked at the top. On the other hand, the ranks of about 90% of the decoys were within the top 10, and seven ligands were incorrectly ranked at the top. These results indicate that the decoys are challenging poses for the current docking score function. Therefore, these poses would be very suitable for investigation of the discrimination power by MD simulation. The detailed RMSD values and rank of the best docking poses are listed in Table S2 of the SI. 2516

DOI: 10.1021/acs.jcim.7b00412 J. Chem. Inf. Model. 2017, 57, 2514−2522

Article

Journal of Chemical Information and Modeling Only the successfully docked ligands (pose0) and the corresponding decoy poses (pose1 and pose2) were considered for the following MD evaluations. In other words, 24 and 22 ligands based on the templates 3DA9 and 3QX5, respectively, qualified and were selected for the following simulations. Docking of Ligands from HSP90a. Similar to the docking of thrombin, the Glide SP mode based on two randomly selected templates was initially performed for HSP90a. Considering the reported crashes among HSP90a structures,16 an additional IFD of the extended sampling protocol was also applied. We found that the IFD protocol showed better posegenerating capability than did the normal SP mode. The detailed results are shown in Table S3 of the SI. It should be noted that although we had advanced knowledge of the crashes of residues around the binding pocket, the default settings of SP mode docking were still employed. We did not define flexible residue(s) to get a better performance in the docking protocol of SP mode because it is not feasible in practice to reasonably and efficiently set the flexible residues without preknowledge about the structural information across proteins. In addition, the purpose of docking here was to generate the correct poses and comparable decoys for simulation. Therefore, we simply selected the poses according to the results from IFD. Figure 3 shows the distribution of minimum RMSD values of all ligands from the docking. Although IFD can generate a good

The template for the homology modeling was based on the structure of CDK1 (PDB code 4Y72, chain A of the highest resolution) because of its highest similarity with CDK2. The sequence alignment with CDK1 is shown in Figure S3 of the SI. Based on the single template, the Prime module of the default settings was used for the geometry building. After that, two loops far away from the pocket were refined by the generalized Born model with the OPLS3 force field.25 The final geometry of CDK2 through homology modeling is shown in Figure S4 of the SI. Surprisingly, the homology model in the current test was quite good. The RMSD values of the binding pocket (the residues within 6 Å around the ligand) and entire protein were about 0.8 ± 0.2 and 0.9 ± 0.4 Å across all 50 PDB files, respectively. In a real case, the accuracy of the structure may be improved by comparing with the structures generated by other popular homology modeling servers.26−29 Here, we did not use other such servers for the structure building as it may have become contaminated by the known structures of CDK2 in these servers. Docking of Ligands from CDK2. We applied IFD with the extended sampling protocol to generate reasonable poses by default settings in the same way as was done for HSP90a. Among the 50 ligands, 46 fulfilled our criterion, where the RMSD value between the best pose from the docking and crystallographic pose was less than 2.0 Å. Figure 4 presents the

Figure 3. Root mean square deviation (RMSD) values of pose0 generated by induced fit docking of heat shock protein 90-alpha.

Figure 4. Root mean square deviation (RMSD) values of pose0 generated by induced fit docking of cyclin-dependent kinase 2.

pose0, the rank from the docking g-score was poor and only one ligand was correctly ranked within the top 10 by g-score. IFD enhances the conformational sampling by allowing flexibility of both ligands and proteins. In such a case, the general score function with default settings might not be suitable, and other score functions like MM/GBSA may possibly provide better predictions.23 This also highlights the importance of our current work, in which the pose discrimination does not solely rely on the docking score function. Finally, based on the IFD results, 36 out of 41 ligands were selected for the following MD simulation. Homology Modeling of CDK2. Recently, Clark et al.24 reported that the binding free energy predicted by using homology modeling could provide comparable accuracy with the prediction based on high-resolution crystal structures. Thus, we examined the discrimination power of the poses by MD simulation based on the structure from the homology modeling as a challenging case. It is not an unusual situation for the structure of a protein to be unknown when a docking is performed. We hypothetically assumed that the structure of CDK2 is unsolved but other homologous proteins have already been reported.

distribution of the minimum RMSD value of the best poses for all ligands. Although the docking was based on the homology model structure, the docking software could also generate very good poses, which were close to those of the corresponding Xray structures. However, the rank by the docking g-score function was somewhat disappointing. About 48% (22 out of 46 ligands) of the ligands in pose0 were correctly ranked within the top 10, and only two ligands in pose0 were accurately ranked at the top. On the contrary, about 63% of the decoys (58 out of 92 ligands) were within the top 10 ranked by g-score function. Therefore, these decoys can be considered to be challenging poses for the score function, and they are appropriate for the examination of the discrimination power by MD simulation. Detailed RMSD values of docking pose0 are listed in Table S4 in SI. Definitions of Stability by MD Simulations. Because the simulation trajectory depends on the initial velocities, it is difficult and risky to judge the stability of the ligand binding mode only from one simulation in a limited simulation time. Therefore, as described in Materials and Methods, five parallel MD runs with different initial velocities were performed for each pose (pose0, pose1, and pose2) to measure the statistical 2517

DOI: 10.1021/acs.jcim.7b00412 J. Chem. Inf. Model. 2017, 57, 2514−2522

Article

Journal of Chemical Information and Modeling

Figure 5. Stability of the docking poses based on two different templates under molecular dynamics (MD) simulation for thrombin. A simulation is classified as “a qualified MD” if the root-mean-square deviation change of the ligand is less than 2.0 Å from the initial structure during the simulation.

stability of the ligand poses. When each MD run was finished, the 500 snapshots from the last 1 ns of the trajectory were extracted for analysis. First, the protein backbone was aligned to the reference structure, which was generated from docking. This is important in practical usage where the crystal structure of the complex is unknown. Then, the average RMSD value of the ligand (heavy atoms only) under MD simulation was calculated. Finally, this average RMSD value was used to evaluate the stability of the ligand. If the RMSD value of a ligand in the MD run was less than 2.0 Å, the ligand was judged to be stable by our definition, and we regarded this MD run as “a qualified run.” In our previous study based on self-docking, if any MD trajectory in the five parallel runs (at least one) was qualified, that ligand pose was then treated as being stable, which may result in overestimation of the ligand stability. This is because the ligand can be accidently stable depending on the initial velocities. We allowed this potential overestimation of the stability to avoid any accidental exclusion of the correct poses by inappropriate initial velocities. Despite that the stability of the decoys was also potentially overestimated by such a criterion, only 56−62% of the decoys were incorrectly considered as stable, which was much less than the values of correct pose (94% pose0 was correctly considered as stable). As a result, about a 32−38% discrimination capability was observed between the correct poses and decoys in our previous self-docking study. On the other hand, in the current crossdocking study, we further evaluated the relationship between the discrimination capability of MD and the number of qualified MD runs. Stability of Ligands from Thrombin under MD. The highest fractions of stable correct poses (pose0) of thrombin were observed for both templates when one qualified MD run was considered (92% and 95% for templates 3DA9 and 3QX5, respectively), as shown in Figure 5. This is reasonable because pose0 was very close to the crystallographic pose. In particular, by the same criterion of stability analysis, the percentage of stable ligands was similar to that found in our previous study based on self-docking (94% and 88% for native and pose0, respectively14). This suggests that our observation in the selfdocking study also holds true in this cross-docking set. Compared with pose0, the decoys (pose1 and pose2) were less stable under MD, no matter which template was used for docking. About a 29−52% discrimination capability of MD could be obtained on average in the current test, when one qualified MD run was considered. This highlights the general discrimination capability of MD simulation.

Figure 5 also shows the fractions of stable poses under MD by the criterion of using different qualified numbers, showing that as the number of qualified MD runs is increased, the fraction of stable poses is decreased for all three poses, as it should be. In this test set, if two or three qualified MD runs were used as the criterion of stability, the discrimination capability of MD seemed to be improved and few correct poses were excluded. The safest choice is obviously to evaluate the ligand stability on the basis of a single qualified MD run, so that pose0 would survive maximally. Maintaining a higher number of correct poses may be more important than excluding a higher number of decoys at the current step, if we use this kind of MD as a filter. In the next step, the accurate but computationally heavy free-energy calculation based on perturbation theory8 may be useful for the further refinement of the poses. However, if a better discrimination capability by MD is really required, then a greater number of qualified MD runs may be used as the criterion. For a successfully docked pose0, the fractions of stable pose0 generated by MD simulation are not perfect, even when one qualified MD is considered. The stability analyses of pose0 for ligand 1D4P (1.3 Å from docking) and ligand 1WAY (1.8 Å from docking) were failed by MD simulation based on templates 3DA9 and 3QX5, respectively. The minimum RMSD value of the 1D4P ligand was 2.4 Å from the MD simulation, which was slightly above the threshold (2.0 Å) may be more appropriate, considering that the protein structure is less accurate in homology modeling and some structural relaxation should be allowed. Another reason might be the poor quality of pose0. Although it is the closest pose to the crystallographic pose in our definition (