|Home | About | Journals | Submit | Contact Us | Français|
Hydrogen atoms are not typically observable in xray crystal structures but inferring their locations is often important in structure-based drug design. In addition, protonation states of the protein can change in response to ligand binding, as can the orientations of OH groups, a subtle form of “induced fit”. We implement and evaluate an automated procedure for optimizing polar hydrogens in protein binding sites in complex with ligands. Specifically, we apply the previously described ICDA algorithm (Proteins 66: 824–837), which assigns the ionization states of titratable residues, the amide orientations of Asn/Gln side chains, the imidazole ring orientation in His, and the orientations of OH/SH groups, in a unified algorithm. We test the utility of this method for identifying native-like ligand poses using 247 protein-ligand complexes from an established database of docked decoys. Pose selection is performed with a physics-based scoring function based on a molecular mechanics energy function and a Generalized Born implicit solvent model. The use of the ICDA receptor preparation protocol, implemented with no knowledge of the native ligand pose, increases the accuracy of pose selection significantly, with the average RMSD over all complexes decreasing from 2.7 to 1.5 Å when applying ICDA. Large improvements are seen for specific classes of binding sites with titratable groups, such as aspartyl proteases.
The two major challenges in virtual ligand screening for drug design are sampling—enumerating possible conformations of ligands in the receptor—and scoring, i.e., estimating the relative binding affinities among different poses of the same ligand and among different potential ligands [1–18]. In general, current docking programs are effective at generating native like binding modes but are less effective at identifying the correct binding mode among alternatives and rank-ordering ligands according to predicted binding affinity [3, 6, 8, 10, 18]. A great deal of effort has been devoted to developing and improving scoring functions for these purposes, using a variety of strategies, often categorized as knowledge based [11, 19, 20], empirical [4, 7, 15, 21, 22], consensus [23–26], or force field based [12, 14, 16, 27–34].
Thus far, knowledge-based and empirical scoring functions have been much more extensively used than force-field based methods. This has been due in part to computational expense associated with computing force field based energies, especially when treating the effects of water using implicit solvent models. However, recently there have been a number of studies evaluating the accuracy of these types of scoring functions using different metrics of success [12, 14, 16, 27–34]. In our own studies, we have noted that the results of using MM-GBSA scoring functions for protein-ligand complexes can be exquisitely sensitive to the treatment of hydrogen atoms in the binding site, a conclusion also reached by other authors . Receptor preparation, including the treatment of hydrogen atoms, is frequently accomplished by manual inspection, but here we investigate whether it is possible to perform this task in an automated way, increasing the robustness and reproducibility of physics-based scoring of protein-ligand complexes.
The problem of assigning hydrogen atoms to crystal structures is an old one and we build upon much prior work. For most experimental structures, the positions of hydrogen atoms cannot be resolved and thus protonation states of titratable residues, the tautomeric states of His, and the orientations of hydroxol or thiol groups are undefined, and may in fact vary depending on what ligand is bound. The orientations of the side chains of Asn, Gln and His may be similarly undefined because of the difficulty in differentiating between N, C, and O atoms from electron densities obtained from crystallography. A number of algorithms have been described for generating all-atom protein models from heavy-atom-only crystal structures, which deal with some or all of these ambiguities, e.g., REDUCE , WHATIF [36,37], and various methods for predicting protonation states/pKa’s based on molecular mechanics and implicit solvent [38–42]. Recently, a new method called the Independent Cluster Decomposition Algorithm (ICDA)  has been described that assigns the ionization states of titratable residues, the amide orientations in Asn/Gln, the imidazole ring orientation in His, and the orientations of OH/SH groups in a unified algorithm that utilizes an all atom protein force field (OPLS) and a Generalized Born continuum solvation model for scoring different possibilities.
Here we use the ICDA algorithm in the context of docking by applying it to protein-ligand binding sites. The atomic level details predicted by ICDA can significantly impact scoring, e.g., errors can result if hydrogen bond donors are incorrectly treated as acceptors, charged side chains are treated as neutral, etc. A common practice in docking and scoring protocols has been to assign protonation and tautomeric states by visual inspection [10, 14, 16], e.g., protonating one of two neighboring aspartic acid side chains in the active site of HIV-protease . Others have modified titratable residues based on their dominant states at neutral pH [1, 5, 45] or based on the pH of the experimental crystal structure or binding assay [44–46]. The drawback of these approaches is that the pKa of a given group depends on its local environment (electrostatic potential, solvent exposure) which may change in response to ligand binding [47–51]. Similarly, the orientations of OH groups in Ser, Thr, and Tyr, and amide groups in Asn and Gln, can also change in response to ligand binding (a subtle form of “induced fit”) .
In this study we examine whether ICDA can be used to improve the identification of native-like binding poses. We use ICDA to determine positions of all atoms undefined by the crystal structure in an automated protocol which does not involve user intervention. Pose selection is performed with a physics-based scoring function based on a molecular mechanics energy function and a Generalized Born implicit solvent model [53–55], the same energy function used in the ICDA procedure. This type of scoring function has been shown to be effective in identifying binders in virtual screening [12, 16, 28] and rescoring conformational ensembles to predict the native ligand pose . Similar Molecular Mechanics/Poisson Boltzmann approaches have been used for the calculation of binding affinities [56–59].
It is important to note that there are several alternative approaches to “physics-based” scoring of protein-ligand complexes, including the problem of assigning protonation states, other than the approaches pursued here. Notably, the PDLD/S-LRA method, which uses the Protein Dipoles Langevin Dipoles approach to treat solvated electrostatics rather than the more approximate Generalized Born implicit solvent models used here, has been shown to accurately reproduce binding affinities in a variety of protein-ligand systems [32–34]. The PDLD approach has also been used to accurately predict pKa’s of groups in proteins [41,42] suggesting that it could provide a self-consistent, physics-based approach to both assigning protonation states and estimating binding affinities.
Our focus is more modest, and we judge success of the ICDA method for preparing the binding site by the ability to subsequently identify near-native conformations of ligands. To focus on the problem of scoring, as distinct from ligand sampling, we test our scoring protocol using an established database of docked decoys [14, 60]. The ability to select native-like poses is tested for a diverse set of 247 protein-ligand complexes taken from the Ligand Protein Database (LPDB) of Brooks et al. . We perform two types of tests. First, we apply the ICDA algorithm to the crystal structure of the protein-ligand complex, and then use that prepared receptor to score the decoys. The native complex would not be known in any realistic application, but this test provides a benchmark for how much improvement may be possible by applying ICDA. We refer to this test below as “biased” (by knowledge of the native structure). Second, we apply ICDA in an “unbiased” way by applying it to every decoy for all 247 cases in the test set (~12,000 protein-ligand structures in total). We compare the results for both of these tests to a “default” receptor preparation procedure that does not optimize protonation states or other optimizations performed by ICDA. We demonstrate that the use of the receptor preparation protocol, in either the “biased” or “unbiased” manner, increases the accuracy of pose selection significantly, with large improvements seen for specific classes of binding sites with titratable groups, such as aspartyl proteases. Overall, our results highlight the importance of optimizing protonation states and the positions of polar hydrogens in binding sites, and present a promising automated method for addressing this challenge as part of a virtual screening protocol.
Our test set is taken from the Ligand Protein Database (LPDB) of Brooks et al. . That database contains 262 receptor ligand complexes, drawn from the Protein Databank (PDB) , and includes 76 unique proteins and 240 unique ligands. For each complex, the LPDB contains 10–120 (50 on average) ligand poses in the receptor binding site, referred to as “active site” decoys. The following complexes are excluded from our study: 1c2t, 1cbs, 1dy9, 1jld, 1mcb, and 1hsh because the ligand file (1mcb), the decoy files (1c2t, 1cbs, 1hsh) or all files (1dy9, 1jld) are missing from the LPDB; 1ac0 and 2acs because we could not create parameter files for the ligand glc (1ac0) and the cofactor nap (2acs); 1bid because the optimization of the native receptor failed; and 1cvu, 1hge, 1hgh, 1hgj, 4hmg, and 5cna because application of the ICDA algorithm failed due to the large size of the receptor (> 900 residues). In total, 247 receptor ligand complexes are included in our study. These complexes include a range of receptor types, including 85 aspartyl proteases, 38 oxioreductases, 31 serine peptidases, 10 immunoglobulins, 14 metallo peptidases, 8 transferases, 6 isomerases, and several others.
All calculations are run in PLOP (free academic version 8.1; commercial version marketed as Prime by Schrodinger Inc.) [63–65] using the all-atom OPLS force field OPLS-AA [53,54] and the Surface Generalized Born implicit solvent model  described in earlier works. Three sets of calculations were performed:
The ICDA algorithm is described in an earlier work  and determines positions of all atoms undefined by the crystal structure, including protonation states for all ionizable residues (His, Glu, Asp, Lys, Arg, Cys, Tyr), ambiguous side chain orientations (His, Asn, Gln) and ambiguous OH/SH group orientations (Tyr, Cys, Ser, Thr). The program assigns polar states in clusters based on neighboring polar residues. The assignments of each cluster are combined and all residues are subject to energy minimization . A pH value of 7.0 is assumed in all calculations.
The unbiased protocol is significantly more computationally intensive than the default or biased protocols, as it requires application of the ICDA algorithm to every ligand pose. For example, for the complex 1a07 with a set of 60 decoys, the total computational expense for the default, “biased”, and “unbiased” protocols were 29 cpu-minutes, 48 cpu-minutes, and approximately 17.5 cpu-hours, respectively. In practical application, we expect that ICDA would not be applied to as many poses per compound as we have done here. Our primary goal here is to provide a proof-of-concept, and we did not focus on algorithmic speed. We expect that it will be possible to improve the speed of the ICDA algorithm as well.
Energies and RMSDs are calculated with the lowest energy decoy being ranked as the top scoring pose. The native pose is not included in the ranked structures.
In the case of the “unbiased” protocol, the protonation states for different decoy poses can vary. To account for this when comparing energies for different poses, we used a standard thermodynamic cycle where the reference state is amino acid dipeptides in solution. For example, a correction factor of −60.65 kcal, the energy difference of charged and neutral glutamate dipeptide, is added for every protonated Glu residue. The pH is accounted for using the term 2.3RT*(pH-pKa). Similar correction factors are applied for all changes in protonation state. The relevant parameters are described in the ICDA paper .
RMSD values are calculated in PLOP and are reported for all heavy atoms in the ligand with respect to the CHARMm minimized ligand taken from the LPDB.
The results of decoy pose selection for the 247 complexes are summarized in Table 1, which reports the average RMSD of the selected, ie: lowest energy scoring, decoy pose for each protocol, as well as the fraction of cases where the selected pose is ≤2 Å RMSD, a commonly used cutoff, and ≤1 Å RMSD, a stringent cutoff. The application of the ICDA algorithm clearly results in substantial improvements in these overall results, and applying it in the “unbiased” protocol results in only a small degradation of the results relative to the “biased” results. This is important, because it suggests that 1) it is not necessary to know the native pose, only to sample conformations near it, and 2) it is not necessary to apply ICDA to all residues in the protein, as we did in the “biased” results, only to the residues close to the ligand, as we did in the “unbiased” results, greatly reducing the computational expense. Our results when using the ICDA optimization are comparable to the best reported in the literature for this dataset; an assessment of scoring functions by Ferrara et al.  reports a success rate (≤2 Å RMSD) of 80–90% for a subset of 189 LPDB complexes. Our results also highlight the utility of an MM-GBSA energy function for pose selection, as also demonstrated in other work [14, 16, 28].
Figure 1 reports the breakdown of RMSD for all protocols. The effect of the ICDA optimization is particularly evident in the significant increase in cases with RMSD ≤1 Å and the decrease in failed predictions with RMSD >4 Å. Table 2 lists 41 cases for which ICDA optimization makes a large difference in pose prediction, defined as a difference in RMSD between the default and unbiased protocols of greater than 2 Å. For each case, the columns list the RMSD of the selected pose, ie: the lowest energy scoring decoy pose, for the “unbiased” ICDA and default protocols. To examine why receptor preparation makes a difference in scoring accuracy in these cases, the table lists the changes in protonation states (from defaults) in the binding site predicted by the ICDA algorithm. Out of 95 changes, 66 changes relate to the protonation of an Asp side chain, while 21 changes relate to a His residue, and 9 changes relate to protonation of a Glu side chain. The exceptionally high RMSD obtained for 1ivp in the default protocol is due to a decoy positioned outside of the protein, far from the binding site.
Table 3 lists the difference in energy between the native and selected pose for all cases listed in Table 2. Energy differences are calculated as Edecoy – Enative, with the native referring to the native ligand pose, subject to the same protocol as the set of decoys. The energy differences obtained from both the unbiased and default protocols show that the lowest energy decoys poses, in most cases, achieve energies lower than the optimized native. This indicates that the decoy sets provide an effective search of conformational space. The ICDA receptor preparation method also dramatically reduces the energy gaps between the native and selected poses in many cases.
Table 4 lists instances of changes in ambiguous Asn/Gln/His side chain orientations, i.e., “flipped” residues, within 3.6 Å of the binding site among the 41 complexes in Table 2. However, the flipped residues do not interact directly with the ligand, and as such likely do not contribute to the accuracy of pose selection.
Table 2 shows that, among cases in which ICDA optimization makes a difference, three receptors dominate: HIV protease (11 complexes), endothiapepsin (9 complexes), and penicillopepsin (8 complexes). These receptors are aspartyl proteases characterized by a catalytic dyad of aspartic acid residues in the active site: Asp 25 and Asp 25′ in HIV protease, Asp 32 and 215 in endothiapepsin, and Asp 33 and 213 in penicillopepsin. In general, experimental and theoretical studies of these complexes [68–72] have shown protonation of one of the Asp residues on the receptor, with the exception of complexes of phosphinate and phosphonate inhibitors in which the PO2 group is stabilized by protonation of both Asp residues .
Our results for HIV protease are consistent with those reported in the literature and show that in general the receptor is protonated on either Asp 25 or Asp 25′; for example complex 1aaq shown in Figure 2a. The only HIV protease complex in Table 2 containing a phosphinate inhibitor is 1hos shown in Figure 2b, and our results are consistent with studies showing that complexes of these inhibitors are protonated on both Asp residues .
For endothiapepsin and penicillopepsin, Table 2 shows the protonation of one dyad Asp residue (Asp 32 in endothiapepsin and Asp 33 in penicillopepsin) for the complexes that do not contain phosphonate or phosphinate inhibitors (1eed, 1epp, 2er7, 3er5, 4er1, 4er2, 1apu, 1apv, 1apw). For the complexes that contain phosphinate and phosphonate inhibitors (1ent, 2wea, 2web, 1ppk, 1ppl, 1ppm), both dyad Asp residues (Asp 32 and 215 in endothiapepsin, and Asp 33 and 213 in penicillopepsin) are protonated, for example 1ppk in Figure 2c. Our results show that the endothiapepsin complexes 3er3 and 5er2 are diprotonated. These complexes do not contain phosphonate or phosphinate inhibitors, or any other type of negatively charged ligand, and as such we have no simple explanation for these results. The other protonation state changes listed for the endothiapepsin and penicillopepsin receptors are for residues that do not interact directly with the ligand, with the exception of Asp 12 which is in close proximity to an acceptor oxygen on the ligands in complexes 2er7 and 3er5.
Table 5 reports average RMSD for complexes that contain a cofactor or ion within 3.6 Å of the ligand. These results indicate successful parameterization of the relevant cofactors in the all-atom OPLS-AA force field [53–54]. Some previous studies have noted challenges in pose selection when a metal cation is coordinated to the ligand in the binding site . While our results do not show higher RMSD when a manganese, magnesium or iron cation is in proximity to the ligand, the presence of a zinc ion does result, on average, in higher RMSD.
The resolution of the crystal complexes in the LPDB ranges from 1.0 to 3.2 Å; 130 complexes have a resolution of below 2 Å. The resolution of the crystal complex is shown to correlate slightly with the RMSD of pose selection. Using the “unbiased” protocol, native-like poses (RMSD ≤2 Å) are selected for 87% of complexes with resolution below or equal to 2 Å, and 82% of complexes with resolution poorer than 2 Å.
Our protocol assumes a pH of 7 in all calculations. We made this choice because pH information is inconsistently reported in PDB files. In reality, however, many complexes were obtained at acidic or basic pH values. This inconsistency may account for some prediction failures. Using the “biased” protocol, a total of 19 complexes result in RMSD > 4 Å. Among these failures, six are for complexes with pH values below 6.5 (1a30, 1afl, 1bai, 1cpi, 1eed, 1epo) and another six (1a07, 1avn, 1ba8, 1cnw, 1cnx, 1cny) are for complexes with pH values above 7.5.
For some complexes in the LPDB , the database contains alternative sets of decoys generated with a modified protonation state on either the receptor or ligand. Running the protocol using these modified structures gives insight into the effect of pH on pose selection. For example, for the protease 1bai which was crystallized at a pH of 6.2, an alternative decoy set was generated using a modified ligand, protonated on a carboxylic acid group. Using the modified ligand and alternative decoy set results in an RMSD of 2.0 Å, in contrast to the unmodified set which results in an RMSD of 16.5 Å.
Several of the failed predictions are for ligands that are significantly solvent exposed, for example several carbonic anhydrase II complexes (1avn, 1cnw, 1cnx, 1cny); see Figure 3 in which the 1cny ligand protrudes out of the receptor structure. In some cases, such as 1a07, the protruding ligand interacts with a second unit of a dimer structure , which is not present in the LPDB complex . In other cases the ligand may be interacting with neighboring molecules within the crystal unit cell.
To compare our results with a state of the art ligand docking program, we ran the Glide docking program (standard precision, version 50207) [76,77] on our test set using the score in place option. Of the 247 complexes we ran through our protocol, we were successful in running 227 complexes in Glide. Twelve cases failed to run because of the large size of the ligand, for example 1fkn (125 atoms), 1hhk (156 atoms), and 3er5 (185 atoms), due to a limit of 35 rotatable bonds in the ligands that can be docked by Glide. Another two cases (1b38, 1b39) failed due to the presence of a magnesium ion. The complexes eliminated from the original dataset include six histocompatability antigen receptors, four endotheapepepsin receptors, two cell division protein kinase receptors and eight other receptor types.
Table 6 shows results for the 227 cases, for Glide and the “unbiased ICDA” and default PLOP protocols. Overall, the results show that the MM-GBSA scoring procedure is only competitive with Glide scoring when the ICDA receptor preparation procedure is applied. Docking and scoring functions such as those used in Glide have been extensively parameterized and validated in a variety of tests, including the ability to reproduce ligand conformations in crystal structures. Many of the proteins in our test set (106 out of the 227 complexes described here) were used in Glide parameterization . In light of this, the fact that the force field based scoring works nearly as well, with no parameterization at all, is an achievement, and lends further evidence to support the potential utility of MM/GBSA scoring functions for inhibitor discovery. Nonetheless, the gap in performance probably does point to some limitations of the physics-based scoring function, most likely the implicit solvent model, for example in treating solvation in cavities . It is also clear that obtaining correct protonation states and hydrogen atom placements is absolutely critical for good performance with the MM/GBSA scoring function. Finally, we note that MM/GBSA scoring is slower than Glide and other empirical scoring functions. We do not anticipate it being used for very large scale in silico screening (except as a secondary “rescoring” protocol ), but the speed is appropriate for, e.g., congeneric series of hundreds or thousands of compounds.
We have demonstrated that optimizing protonation states and the positions of polar hydrogens in binding sites can play a significant role in selecting native-like ligand poses. We used a scoring function based on a molecular mechanics force field (OPLS) and continuum solvent (GB) for selecting poses, lending further evidence that physics-based scoring functions of this type can be usefully employed in docking protocols [28–29]. We believe that it will also be possible to incorporate the optimization of protonation states into protocols for “induced fit” docking that treat conformational changes in the protein due to ligand binding [44, 79–81]. Receptor preparation can in principle also be conducted in conjunction with treatment of the ligand protonation states using a program such as Epik  which uses empirial pKa prediction to generate possible protonation states of ligands under specified conditions. Finally, future studies will examine the application of the protocol to virtual screening and ranking ligands in congeneric series .
We thank Dr. Kai Zhu for help with the ICDA algorithm. This work was supported in part by the Sandler Program in Basic Sciences, and NIH grants GM071790 and AI035707 (to MPJ). MPJ is a consultant to Schrodinger LLC.