|Home | About | Journals | Submit | Contact Us | Français|
Several studies have suggested that disrupting interactions of the G protein βγ subunits with downstream binding partners might be a valuable study for pharmaceutical development. Recently, small molecules have been found which bind to Gβγ with high apparent affinity in an enzyme-linked immunosorbent assay (ELISA), have demonstrated selective inhibition of interactions of Gβγ with downstream signaling partners, and have been shown to increase antinociceptive effects of morphine and inhibit inflammation in vivo. In this paper we examine several docking and scoring protocols for estimating binding affinities for a set of 830 ligands from the NCI diversity set to the Gβ1γ2 subunit and compared these with IC50s measured in a competition ELISA with a high-affinity peptidic ligand. The best-performing docking protocol used a consensus score and ensemble docking and resulted in a 6-fold enrichment of high-affinity compounds in the top-ranked 5% of the ligand data set.
G protein-coupled receptors (GPCRs) trigger cellular signaling cascades through interactions with heterotrimeric guanine nucleotide binding proteins (G proteins) and are involved in many diverse physiological processes such as sensory perception, modulation of cardiac rhythm, neurotransmission, attraction of motile cells by chemotaxis, and regulation of mitosis.1,2 GPCRs are common targets for pharmaceutical intervention; around half of all drugs in clinical use are GPCR agonists or antagonists.3 In the standard model for signaling, the receptor acts as a catalyst for the exchange of GDP for GTP on the Gα subunit, which leads to dissociation from the Gβγ subunits. Gα and Gβγ are then free to interact with downstream signaling partners. Five distinct β subunits and twelve distinct γ subunits have been identified in the human and mouse genomes; these subunits can pair to form many Gβγ different combinations. Evidence suggests that the different Gβγ isoforms interact with different GPCRs, although this selectivity is not completely understood.4
Several studies have suggested that disruption of interactions of Gβγ with downstream binding partners might be a valuable strategy for pharmaceutical development.4–21 Recently, small molecules have been found which bind to Gβ1γ2, as demonstrated by the ability to compete with a high-affinity peptide ligand (SIGK) in an enzyme-linked immunosorbent assay (ELISA). An X-ray crystal structure of SIGK bound to Gβ1γ2 has been determined and shows that SIGK binds to the same region of the β subunit as the switch II region of Gα, in the center of the β-propeller.22 This region has been postulated to be a hotspot or a small area on a protein–protein interaction surface particularly important for mediating binding. The small molecules were discovered with the help of a small virtual screen of ~2000 compounds from the National Cancer Institute (NCI) diversity library, using the FlexX software package.23 One compound, M119 (NSC119910), was subsequently shown to inhibit Gβγ-dependent activation of phospholipase C and phosphoinositide 3-kinase γ (PI3K) and to increase the antinociceptive effects of morphine in mice.21 A very similar compound, gallein (which is commercially available at high purity), was found by surface plasmon resonance to bind Gβ1γ2 with Kd ~ 400 nM and was shown to inhibit PI3K activation and chemoattractant-dependent neutrophil migration and lessen carrageenan-induced inflammation in mice.20 In general, disruption of protein–protein interactions has the potential to be an important strategy for developing therapeutics. Although contact surfaces are often large in size and relatively flat,24–26 studies involving point mutations show that only a small subset of residues contribute most of the free energy of binding; these “hotspots” are often located at the center of the interface.27–33 Recent reviews of the principles of protein–protein interactions and their significance for drug design have been published.34,35
Virtual screening by docking compounds from chemical libraries has become widely used in drug discovery, for finding new lead compounds when a high-resolution structure of the target is available.24–26,36–40 Many docking programs have been developed over the past few decades, and good success in screening applications has been reported in some cases. However, accurate prediction of binding affinities remains challenging.41–46 All recently available programs provide false negatives (active compounds that are not given a high-ranking score by docking programs) and false positives (inactive compounds that are given a high-ranking score). Individual programs often show better ability to predict high-affinity compounds for particular receptors or compound classes. It is difficult to say in advance whether a particular program or scoring function will be successful for a certain case without any supporting experimental knowledge about the target.
In this paper, we examined several docking and scoring protocols for estimating binding affinities for a subset of 830 ligands from the NCI diversity set to the Gβ1γ2 subunit and compared these with IC50s measured in competition ELISA with the high-affinity SIGK peptidic ligand. The DOCK6 and GLIDE software packages were used, and the influences of the assignment of protonation states and partial charges of ligands and receptors were evaluated. Enrichment factors from both software programs (the number of high-affinity ligands found as a fraction of the amount of the database screened) were compared depending on variations. The influences of the assignment of protonation states and partial charges of ligands and protein receptor were evaluated. Solvation effects were introduced using the Poisson–Boltzmann continuum method. The effects of structural flexibility were explored in detail by molecular dynamics (MD) simulation. The aim of this study is to obtain the best-optimized docking method to apply to a larger database.
830 compounds taken from the NCI diversity set (comprising a total of 1990 compounds) were tested for competition with the SIGK peptide in ELISA as described in refs 16 and 21. Compounds with less than 300 μM of median inhibitory concentration (IC50) value were assumed to be high-affinity compounds; a total of 36 such high-affinity compounds were found in this assay (see Table 1).
The crystal structure of Gβγ complexed with the SIGK peptidic ligand (PDB ID 1XHM) was separated into protein and ligand PDB files. The protein was processed to add hydrogen atoms and remove water molecules and minimized to alleviate bad steric contacts with the Protein Preparation Wizard utility in Maestro (Schrödinger, Inc.). Two different sets of partial charges were employed for the receptor: AMBER charges,47 assigned using the CHIMERA program,48 and OPLS/AA charges,49 which were automatically assigned by Glide version 4.5 (Schrödinger, Inc.) in the receptor grid calculation. The protein–ligand file was subsequently used to select spheres in the hotspot of the receptor in the process of receptor preparation by the DOCK6 program.50
Coordinates for molecules in the NCI diversity set were downloaded from NCI. Metal ions were removed and hydrogens were added using SYBYL version 6.8 (Tripos, Inc.). Tautomers were generated using ionization states estimated for pH 7 with the LigPrep program (Schrödinger, Inc.). Four different sets of partial charges (CM1,51 CM2,52 AM1-BCC,53 and Mulliken54) were assigned to each ligand using the Antechamber program in the AMBER software package. OPLS ligand charges were also obtained using Glide by the default setting.
Using the DMS program in the DOCK6 software package, the solvent-accessible molecular surface of the protein binding site was calculated using a probe radius of 1.4 Å. Receptor spheres were generated using the program SPHGEN. Spheres covering the hotspot were selected within 10 Å from the positions of the heavy atoms of the SIGK peptide ligand. The grid box enclosing the selected spheres was generated with an extra 5 Å added in each dimension. Ligand flexibility was employed during the docking process. When performing docking with Glide, the inner box, which defines the range of motion for the center of each ligand, was set to 10 Å on each side. The outer (enclosing) box, which includes the entire hole at the center of Gβγ, was set to 11 Å longer than the inner box on each side. Docking calculations were performed using both standard-precision (SP) and extraprecision (XP) modes.
The Poisson–Boltzmann continuum method implemented in MEAD55 was used to introduce solvation effects. The best poses of ligand molecules from DOCK6 were used for solvation free energy calculations. A grid spacing of 1.0 Å was applied, and a focused grid of 0.25 Å was utilized. The differences of solvation free energies of ligands in protein and in water were calculated.
A relatively simple way to take in to account the effects of receptor flexibility is to dock to an ensemble, rather than a single structure.56–58 To create an ensemble of receptor conformations for docking, molecular dynamics (MD) simulations of the unbound receptor were performed with the SANDER module of the AMBER9.0 software package. The parm99 force field parameters were used, and solvent effects were incorporated using the pairwise generalized Born model of Case and co-workers.59–61 Prior to a production run, the receptor was first relaxed by steepest descent energy minimization, followed by conjugate gradient energy minimization. The system was then heated over 25 ps from 0 to 300 K. An MD simulation of 1 ns was performed to produce a conformational ensemble of the receptor at 300 K. A nonbonded cutoff distance of 12 Å and a time step of 2 fs was used. Constant temperature was maintained using the Langevin thermostat, and all bonds involving hydrogen atoms were constrained using the SHAKE algorithm. Coordinates were saved every 1000 steps, for a total of 500 snapshots. The snapshots were classified into ten groups by the fuzzy c-means (FCM) clustering algorithm using the GA Fuzzy Clustering program.62 From each cluster, the structure closest to the center was chosen for docking. A docking score was obtained for each of the ten conformations. These scores were then averaged to rank compounds.
Enrichment factors (i.e., the number of high-affinity compounds found as a function of the amount of the ranked database examined) were calculated using both unrefined original forms of NCI database ligands and ligands with protonation states corrected for pH 7, using LigPrep. AM1-BCC charges were used for the ligand, with AMBER charges for the receptor. When the structures in the unrefined original ligand database were inspected, most of the titratable functional groups such as carboxyls were found to be protonated. The charges on the acid functional groups were also significantly different from those on the refined ligands. However, most of the basic functional groups such as amines were deprotonated in both the refined and unrefined ligand structures. Docking was performed with DOCK6 using the Grid scoring function as described above. Enrichment factors are shown in Figure 1. With unrefined, original protonation states, DOCK6 could not effectively identify active ligand compounds, since the enrichment curve was only slightly better than random selection. After the ligand protonation state was corrected, the percent known active ligands found in the top tenth of the ranked database was increased from 13% to 33%, suggesting that the protonation state of ligands has a significant effect on molecular recognition by the Gβγ active site.
To investigate the role of partial charges in docking and scoring accuracy, charges for Gβγ were assigned using both AMBER99 and OPLS/AA force fields; four different sets of charges (CM1, CM2, AM1-BCC, and Mulliken) were assigned to the ligands. The DOCK6 program was used to generate poses and obtain scores. Results are shown in Figure 2 and suggest that overall, more accurate results are obtained using AMBER receptor charges. Using AM1-BCC, 28% of the known high-affinity ligands were found in the top 10% of the database using OPLS receptor charges, while 33% of known ligands were found using AMBER receptor charges. In general, AM1-BCC charges resulted in the greatest enrichment, followed by Mulliken, CM2, and CM1. This trend held for both sets of receptor charges. The influence of the receptor charges was less significant than that from the ligand charges.
Success in enrichment and in accurately reproducing correct binding poses has been reported for the Glide docking program.63–65 We computed enrichment factors with Glide; results are shown in Figure 3. The best enrichment (33% of high-affinity ligands found in the top 10% of the database) was slightly higher with DOCK6 than with Glide (21% of high-affinity ligands found in the top 10% of the database). Combinations of AMBER and OPLS receptor charges and AM1-BCC and OPLS ligand charges were used. Similar to results with DOCK6, the AMBER receptor charges were generally superior to the OPLS receptor charges when using Glide. When OPLS charges were used for the ligands, Glide SP with AMBER receptor charges found 21% of known ligands in the top 10% of the ranked database, while Glide SP with OPLS receptor charges found only 12% of known ligands. When AM1-BCC charges were used for the ligand and AMBER charges were used for the receptor, the fraction of actives found in the top 10% dropped from 21% to 9%. Results indicate that the AMBER/OPLS receptor/ligand charge combination is the most effective to identify known active compounds for this system with Glide. However, when the AMBER/OPLS receptor/ligand charge combination was applied to DOCK6, the enrichment in the top 10% of the database was less than that from the AMBER/BCC receptor/ligand charge combination, as shown in Figure 2a.
The Glide extended-precision (XP) mode has previously been reported to be more effective for finding known active compounds for a large and diverse set of ligands and receptors than Glide SP.66 However, it was designed as a refinement tool for use with relatively good ligand poses from Glide SP because it demands more intensive computational expense. We took the best-scoring poses from Glide SP mode using the AMBER/OPLS combination of charges and docked again using Glide XP mode. This actually provided degradation of the fraction of actives found in the top 10% of the database from 21% to 12%. However, in the top 30% of the database, enrichment was essentially the same.
To account for solvation, we used the Poisson–Boltzmann continuum method implemented in MEAD. Because PB energies are known to be very sensitive to receptor–ligand geometry, only the best poses from previous docking with DOCK6 were considered. The ligand charge set used in this method was AM1-BCC, which gave the best enrichment out of the charge sets we considered. The solvation free energy differences of the ligand in protein and in water were added to the DOCK6 Grid scoring function. (Solvation free energy calculations were not included for Glide, because the Glidescore already includes a solvation term.) Calculated solvation free energy differences were small and had little effect on enrichment factors.
We constructed a simple consensus score by adding together the grid score from DOCK6 using AMBER/AM1-BCC receptor/ligand charges and the Glidescore from Glide SP with AMBER/OPLS charges. In this case, the fraction of actives in the top 5% of the database increased to 28%, and the fraction of actives found in the top 10% was increased to 37%. When the correction for solvation free energies was added to the consensus scoring function, the enrichment factors in the top 5% and 10% of the database were the same (Figure 4).
Ensemble docking—that is, docking to ensembles of conformations of receptors taken from molecular dynamics simulations—was performed to account for receptor flexibility. MD simulations did not significantly change the overall receptor structure but produced fluctuations (Figure 5). The receptor and ligand charge combinations used were AMBER/AM1-BCC for DOCK6 and AMBER/OPLS for GlideSP, which resulted in the best enrichment as described above. In general ensemble docking improved enrichment but only by a slight amount (Figure 6). For DOCK6, the fraction of actives found in the top 10% of the database was slightly improved (from 33% to 38%). For Glide SP, ensemble docking resulted in fewer actives found in the top 10%. However, consensus scoring with a combination of the DOCK6 grid score, solvation correction from MEAD, and the Glidescore resulted in an increase of the number of actives found in the top 10% of the database to 38%. The result shows that the performance with ensemble docking is better for the top fraction of the database, which is most valuable as it is this fraction that is likely to be tested in practice.67 Several studies have shown that docking scores are often highly correlated with molecular weight, indicating that nonspecific interactions are the most significant contributors to scores.67 Therefore, for comparison, we calculated enrichment simply ranking compounds by molecular weight. Results indicate that a consensus score with ensemble docking provides superior enrichment in the top fraction of the database; however, in considering larger fractions (greater than the top fifth or so), a simple ranking by molecular weight does just as well.
In this work we examined the ability of docking and scoring methods to predict the binding of small molecules to a protein–protein interface, the hotspot of the G protein β subunit. The influence of variations in the ligand protonation state, initial structures of proteins and ligands, partial charges, scoring functions, and docking protocols was investigated using the DOCK6 and Glide programs. Docking to a broad protein–protein interaction surface is a difficult problem, but it was possible to obtain significant enrichment in the number of active compounds found in the top-ranking portion of the database. The use of a consensus score and taking receptor flexibility into account by docking to an ensemble of configurations from molecular dynamics simulation resulted in greater enrichment. Limitations remain, such as the neglect of changes of conformational entropy upon binding. Furthermore, docking programs could not directly rank compounds according to affinities, and there was little correlation observed between scores and measured inhibition constants. For example, the highest-affinity compound M119 (IC50 = 200 nM) was filtered out of the top 10% of the database during this screening process, although it appeared in the top 30% for all methods. In agreement with previous studies, enrichment provided by docking is comparable to ranking compounds simply by molecular weight, although ensemble docking using a consensus score was better at finding high-affinity compounds in the top-ranked fraction of the database.
Somewhat surprisingly, the binding area of small molecules observed in docking does not always overlap with the binding of the high-affinity SIGK peptide ligand (Figure 7). In general, small compounds were observed to dock inside the hole of the Gβγ hotspot, rather than on the surface. While this might be an artifact of the modeling, if poses are reasonable accurate, it suggests that small molecules can disrupt the interaction without completely occluding the binding surface.
Future work will involve applying the methods investigated here to larger databases of millions of compounds. It is expected that docking/scoring methods will be able to make a contribution toward discovery of new high-affinity ligands with the potential for selective inhibition of interactions of Gβγ with downstream signaling partners as well as to a better understanding of the structural basis for signaling specificity.
We thank Dr. Orr Ravitz and Cen Gao for their assistance with computational problems. This work was supported in part with funds from the Camille and Henry Dreyfus Foundation.