|Home | About | Journals | Submit | Contact Us | Français|
With the recent breakthroughs in G protein-coupled receptor structure, one can now compare experimentally determined structures with the most recent modeling and docking methods. A community-wide blind prediction experiment (GPCR Dock 2008) was conducted in coordination with the publication of the human adenosine A2A receptor bound to the ligand ZM241385 crystal structure (Science 322, 1211 (2008)). Twenty-nine participating groups submitted 206 models that were evaluated for the accuracy of the ligand binding mode and the overall receptor model. Several new insights emerged including the critical importance of disulfide bonds in the extracellular loops, helix residue registry, and domain knowledge.
Molecular modeling has a pertinent role in rational drug discovery and design1, 2. Reliable three-dimensional models can provide valuable insights into basic principles of molecular recognition and aid in the structure-based drug design approach to lead finding and optimization3. G protein-coupled receptors (GPCRs) are membrane proteins involved in signal transduction pathways and are important therapeutic targets for numerous diseases4, 5. As such, significant structure prediction efforts using methods ranging from de novo to homology-based approaches have been applied to members of the GPCR family6, 7. Until recently, much of GPCR homology modeling efforts have been based on the templates of bovine rhodopsin and bacteriorhodopsin, with refinement of the models by molecular dynamics simulations, ligand docking, and incorporation of additional biochemical and biophysical data8–12. The refinement step is necessary in building accurate models, especially around the ligand-binding site, due to the expected structural differences both among members of the family resulting from the generally low sequence identity and the large diversity of ligands accommodated within the family7, 13–15, and among various conformational states associated with the different ligand efficacy16–18.
The most recently solved GPCR structure is the 2.6 Å crystal structure of the human adenosine A2A receptor bound to an antagonist19. Adenosine receptors belong to the class A rhodopsin-like GPCR family, and are implicated as promising therapeutic targets in a wide range of conditions, including cerebral and cardiac ischaemic diseases, sleep disorders, immune and inflammatory disorders, and cancer20. The A2A structure shows an overall seven transmembrane helix structural similarity to the rhodopsin and adrenergic receptor structures, with shifts in the positions and orientations of the helices and a markedly different structure of the extracellular loops19.
To evaluate the current progress in GPCR structure prediction and docking, we carried out a community-wide blind prediction experiment (GPCR Dock 2008) in coordination with the release of the human adenosine A2A receptor structure in October 200819. GPCR Dock 2008 was organized in a similar manner as the CASP and CAPRI experiments21, 22, and was aimed to assess the current status of GPCR structure modeling and docking, as well as to highlight areas for future efforts in method development. In this paper, we report the results of the assessment together with our analysis of the current status of GPCR structure and ligand docking predictions.
In August 2008, prior to the publication of the human adenosine A2A structure in October 200819, participating predictors were asked to blindly predict and submit up to ten ranked models of the human A2A receptor in complex with the ligand ZM241385, starting from the amino acid sequence of the receptor and a 2D structure of the ligand. A total of 63 different individuals initially registered, with 206 models submitted by 29 different individuals in the final data set. Note that 37 of the 206 submitted models were either missing the ligand or had incorrect bond connectivity for the ligand. We assessed the remaining 169 models for the prediction accuracy of the ligand binding mode, and all 206 models for the prediction accuracy of the receptor alone.
RMSD (root mean square deviation) is used as a quantitative measure of the similarity between two superimposed atomic coordinates, and is calculated by
where δ is the distance between N pairs of equivalent atoms from the two coordinates. RMSD values (units of Å) can be calculated for any type and subset of atoms, e.g. Cα atoms of proteins (Cα RMSD) for all residues, for residues in the transmembrane helices or the loops; heavy atoms of small-molecule ligands (ligand RMSD).
Z-score is a standard dimensionless score that normalizes a value with respect to the sample mean and standard deviation, and is defined by
where x is the raw value, μ is the mean and σ is the standard deviation.
Assessment criteria are dependent on the purpose of the generated models. Given the primary value of the GPCR structural models in expanding our knowledge in basic molecular recognition and their potential use in design and development of new small molecules, the quality of the models was primarily assessed by the accuracy of the ligand binding mode. Particular attention and care was made towards the fact that the crystal structure is a static structure with positional errors, and the value of modeling is ultimately to guide drug discovery and biological insight. Our numerical measure of accuracy for the ligand binding mode was based on two metrics, ligand RMSD and the number of correct receptor-ligand contacts. Neither metric alone was sufficient to capture the accuracy of prediction around the ligand-binding site; hence, both were used and combined into a z-score scoring scheme to rank the models.
The ligand RMSD between model and crystal structure is calculated as the coordinate root mean-square deviation (RMSD) for the 25 non-hydrogen atoms of ZM241385 after superimposing the protein Cα atoms of the model and the crystal structure. In addition, the ligand RMSD is also calculated excluding the phenoxy group of ZM241385 that has high B-factor values. The number of correct contacts is counted as the number of correctly predicted native contacts observed between protein atoms and the ligand. A native contact is defined as any inter-atomic distance within 4 Å of the ligand in the crystal structure. There are 75 such receptor-ligand contacts, and an additional 15 contacts formed with water.
The models were ranked by assigning a combined mixed z-score to each model. The combined z-score was calculated as the average of z-scores for ligand RMSD and the number of correct contacts: Zcombined=(−ZLigandRMSD+ ZN _ CorrectContacts)/2. The z-scores for ligand RMSD and the number of correct contacts were computed in two passes as follows: i) assign a z-score to each model using the average and standard deviation values from all models, ii) re-compute the average and standard deviation excluding models with z-scores more than two standard deviations above (for ligand RMSD) and below (for the number of correct contacts) the average, iii) re-assign a z-score to each model using the revised average and standard deviation values obtained in step ii. The best model, i.e. the model with the highest combined z-score, from each group was analyzed.
The submitted models show a wide distribution in prediction accuracy of the ligand binding mode, with average values of 9.5 Å (s.d. 3.8 Å) for ligand RMSD and 4 (s.d. 7) for the number of correct contacts (Figure 1A).. These statistics indicate that the majority of the submitted models do not predict the ligand position and the binding interactions very accurately. Furthermore, the lack of strong correlation between ligand RMSD and binding site RMSD (Figure 1B), e.g. models with <4.0 Å binding site RMSD have a range of 2.8 to 17.2 Å ligand RMSD, suggests that the performance of some ligand docking methods is poor and can be improved.
Very few models score well in both ligand RMSD and the number of correct contacts (only 13 out of the 169 total receptor-ligand models have >1 combined z-score, compared to 40 models that score well solely in ligand RMSD, ZligandRMSD <−1.0). For models with relatively low ligand RMSD values but small number of correct contacts, the inaccuracy in binding interactions may be attributed to errors in the sidechain placement of the ligand-binding residues. Although nearly a third of the models capture the hydrogen bonding interaction between the N2536.55 sidechain and the exocyclic N15 atom of the ligand (44 out of 169 models have <4 A N253 OD1 – ZM241385 N15 interaction distance), other key receptor-ligand interactions such as the aromatic stacking interaction between the F1685.29 sidechain and the bicyclic ring of the ligand are not captured well in most models (Figure 2).
While the overall outcome clearly shows that there are remaining challenges in accurately predicting the ligand binding mode, the quality of the predictions for the receptor alone appear to be much higher: 4.2 ± 0.9 Å for the receptor Cα RMSD, and 2.8 ± 0.5 Å for the transmembrane (TM) helices Cα RMSD; loop regions, with the exception of ICL1 (intracellular loop 1), are not modeled accurately in the majority of the models (Figure 3A,B and Figure 4). It is notable that some groups that accurately predict the TM region of the receptor do not predict the ligand binding mode very well (TM Cα RMSD: 2.0 Å for Pogozheva/Lomize, and 2.1 Å for Horst/Roy), indicating that the methods for modeling the receptor and docking the ligand can be generally considered as distinct steps in the generation of models for the receptor-ligand complex.
Despite the apparent difficulty in accurately predicting the receptor-ligand interactions, some models had consistent features with the crystal structure, although model ranking continues to be one of the most challenging areas of development. Here, we focus on the predictions from the top ten groups, ranked according to the combined z-score, and assess the model quality at greater details (Figure 5). Note that, with predictions for only one target, the statistical significance of the group ranking cannot be judged, as is typically done in CASP experiments by a head-to-head comparison on common targets between the top groups23. To further support our selection of the best predictions, we ranked all models using an alternative metric, binding site contactRMSD, which treats all ligand-binding residues with equal weight, and is an RMSD of receptor-ligand contact distance over all ligand-binding residues. We found that both the z-score ranking and the contactRMSD ranking agree on the selection of the best model, and the majority of other top predictions.
The best model (Costanzi) among all submitted models has a ligand RMSD of 2.8 Å RMSD and 34 of 75 correct contacts (Figure 6A and Table I). The ligand is modeled in a native-like binding pose, with an extended conformation and a nearly perpendicular orientation to the membrane plane. The model accurately predicts some of the key receptor-ligand interactions: it captures the hydrogen bonding interaction between the N2536.55 sidechain and the exocyclic amino group (N15 atom) of the ligand, and the aromatic stacking interaction between the F1685.29 sidechain and the bicyclic triazolotriazine core of the ligand. Compared to the crystal structure, the ligand in the model is positioned deeper in the binding pocket, bringing the furan ring closer to TM helices III and V. The inaccuracy in the ligand position is most likely due to errors in the sidechain positions of the two crucial ligand-binding residues (F1685.29 and E1695.30) in ECL2 (extracellular loop 2) and the sidechain orientation of M1775.38 at the extracellular end of TM helix V: the aromatic ring of F1685.29, which interacts with the bicyclic ring, is positioned too deeply; the adjacent E1695.30 forms a hydrogen bonding interaction with the hydroxyl group in the phenolic substituent instead of the exocyclic N15 atom near the bicyclic ring; the sidechain of M1775.38 is not oriented toward the binding cavity. In addition, the family conserved disulfide bond C773.25 – C1665.27 is predicted accurately, but the disulfide bond in ECL3 C2596.61 – C2626.64 is not, presumably contributing to the inaccuracy in the sidechain orientation of H2646.66, which is not pointed toward the binding site.
The best predictions from the top six groups (Costanzi, Katritch/Abagyan, Lam/Abagyan, Davis/Barth/Baker, Maigret, Jurkowski/Elofsson) highlight the successes and difficulties in accurately predicting the ligand binding pose and receptor-ligand interactions (Figures 6B,C,D and Table I). The extended ligand conformation is accurately predicted in all six models, and the nearly perpendicular orientation is captured in four of the six models. The hydrogen bonding interaction between the N2536.55 sidechain and the exocyclic N15 atom of the ligand is correctly modeled in four models; however, in one of the four, the ligand makes no interaction with residues in ECL2. The aromatic stacking interaction between the F1685.29 sidechain and the bicyclic ring of the ligand is correctly modeled in four models; however, in all four models, the ligand is positioned too deeply in the binding pocket, and the M1775.38 sidechain is not oriented toward the binding cavity. There is one model that does not accurately capture either the hydrogen bonding interaction with N2536.55 or the aromatic stacking interaction with F1685.29, while five of the six models accurately predict the family conserved disulfide bond C773.25 – C1665.27. None of the six models capture the hydrogen bonding interaction between E1695.30 in ECL2 and the exocyclic N15 atom of the ligand.
Other models that ranked near the top (Kanou, Goddard, Bologa, Olson) are slightly less accurate but show similar trends as the top six models in their ability to accurately predict the ligand binding mode (Table I). The ligand is modeled in a native-like extended conformation in three of the four models. The hydrogen bonding interaction between the N2536.55 sidechain and the exocyclic N15 atom of the ligand is modeled accurately in three of the four models; whereas, the aromatic stacking interaction between the F1685.29 sidechain and the bicyclic ring of the ligand is modeled accurately in only one of the four models. The family conserved disulfide bond C773.25 – C1665.27 is captured in two models. Remarkably, one of the models (Goddard) accurately places the E1695.30 sidechain proximal to the exocyclic N15 atom of the ligand, and almost captures the hydrogen bonding interaction, even though the overall ECL2 conformation is inaccurate.
The best predictions were generally not ranked as the best models by the predictors at the time of model submission prior to the release of the crystal structure (Table I). Only two of the six best models were ranked first, and three of the six groups show a weak correlation between their model ranking and the model quality as assessed by the combined z-score for the accuracy around the ligand-binding site. Furthermore, the additional models submitted by the six groups are generally of lower quality than the best predictions (Table I). Only one of the six best models has a z-score that is within one standard deviation of the group average z-score.
The assessment of the submitted models shows that the best participating methods have the ability to predict close native-like ligand binding, but with limitations in capturing all of the key receptor-ligand interactions and correctly estimating model quality by ranking. The majority of the submitted models are quite far from achieving a native-like ligand binding pose. The most challenging aspect of GPCR structure prediction highlighted in this assessment seems to be in accurately modeling the ligand interactions with residues in the extracellular loop regions. This result is not surprising given the lack of structural homology in the loops among the known GPCR structures24, and the generally difficult task of modeling loops25, 26.
The most successful prediction methods relied on homology modeling approaches based on the template structures of β-adrenergic receptors, and in some cases with the additional template structures of rhodopsin, (PDB IDs: 2RH1 (β2AR), 2VT4 (β1AR), 1U19 (bovine Rhod), 2Z73 (squid Rhod)) to generate models of the receptor, followed by docking of the ligand to one or more receptor models using small-molecule docking programs such as Glide27, ICM28, GOLD29, and AutoDock30 (see Supplementary Information for description of prediction methods). The alignment of the human A2A sequence to the template structure seemed to have been straightforward, given the family conserved motifs and residues in the TM helices31. The ECL2 was modeled by de novo approaches in many of the top predictions (Katritch/Abagyan, Lam/Abagyan, Davis/Barth/Baker, Jurkowski/Elofsson, Goddard), but only partially modeled in the best prediction (Costanzi) for a short segment of eight residues, located N-terminus to TM helix V, that includes the disulfide-bond-forming C1665.27. Some of the criteria used to select and rank the final receptor-ligand complex models were: docking scores, conformational energy of the complex, agreement with mutagenesis and structure-activity relationship data, and binding selectivity studied by virtual ligand screening or by modeling other subtypes of adenosine receptor.
The reliability of the homology modeling approach depends on the availability of suitable templates32. The results of the current assessment show that the structures of β-adrenergic receptors alone or together with rhodopsin were suitable transmembrane templates in predicting the general structure of the adenosine A2A receptor. However, given the expected structural diversity in class A GPCRs, it is unclear whether the current set of techniques applied to the structure prediction of the A2A–ZM241385 complex would result in a similar level of accuracy for the prediction of other GPCRs, especially for those belonging to subfamilies that are phylogenetically distant from the amine and the opsin receptor clusters33. We believe the database of GPCR structures needs to expand further to provide suitable templates for accurate modeling of those other receptors.
The inaccuracies in homology models can arise from errors in side-chain packing, main-chain shifts in aligned regions, errors in unaligned loop regions, misalignments, and incorrect templates34. These errors relate to the issue of “adding value” to the template structure, which was addressed in the recent CASP experiment35, and also seems to be applicable to GPCR modeling. Indeed, ligand interactions with residues located in structurally divergent regions from the templates are consistently not modeled accurately in all of the six best predictions: hydrogen bonding interaction between E1695.30 in ECL2 and the exocyclic N15 atom of the ligand is not captured, and the sidechains of H2646.66 in ECL3 and M1775.38 in the extended bulge structure unique to A2A at the extracellular end of TM helix V are not oriented toward the binding site. An exception is the aromatic stacking interaction between F1685.29 in ECL2 and the bicyclic ring of the ligand, which is correctly modeled in some of the predictions. F1685.29 is located in the loop, but it is structurally homologous to F1935.32 which interacts with the carbazole heterocycle of the ligand carazolol in the β2AR structure, hence modeling of this interaction may have been guided by homology. Interestingly, F1685.29 is modeled more accurately than E1695.30 even though mutagenesis data shows mutation of E1695.30 to alanine reduces the affinity for both antagonists and agonists36, and no data is available for F1685.29. The inaccuracy in the orientation of the ligand binding pose, e.g. the parallel orientation with the phenolic substituent positioned close to TM helices II and III, may in part be due to the inaccurate modeling of the helical shifts in TM helices I, II, and III. The helical shifts alter the location of the binding pocket and redefine the pocket size and shape19; thus, it is expected that accurately modeling the helical shifts would contribute to a better prediction of the ligand binding pose. The helical shifts were most accurately modeled by an effective use of multiple template structures of rhodopsin and β-adrenergic receptors (Pogozheva/Lomize), or an all-atom refinement approach implemented by the ROSETTA program using a physically realistic model that recapitulates protein interatomic and protein-solvent interactions in the membrane environment37 (Davis/Barth/Baker).
Other sources of error include not modeling the water molecules that are either structurally important or directly involved in ligand binding interactions3. The ligand binding cavity in the A2A–ZM241385 structure has four ordered water molecules19, yet none of the submitted predictions included water molecules. We tried re-docking the ligand to the crystal structure using ICM28 and found that a native-like binding pose (within 1 Å heavy atom RMSD for the bicyclic ring and the furanyl substituent of the ligand, and <3.0 Å overall ligand RMSD) can be recovered without any water molecules, which suggests that water may not be critical for accurately predicting the ligand interactions. However, modeling water molecules together with the ligand may contribute in general to a better prediction of the ligand binding pose or affinity. Additional re-docking studies with the docking protocols used by the participating methods would help assess the effect of the water molecules, and the accuracy of the docking methods separately from that of the receptor modeling methods.
Lastly, it is interesting that the best model was from the Costanzi group who worked previously with adenosine receptor modeling and docking. Their domain knowledge on the adenosine receptor may have helped with the interpretation of the mutagenesis and ligand interaction data and such domain knowledge is almost certainly critical.
Accurate prediction of GPCR structure and ligand interactions clearly remains a challenge despite the increase in the number of experimentally solved GPCR structures within the last year. Assessment of these predictions highlights similar issues addressed by the CASP predictions for template-based modeling targets, i.e. the difficulty in loop modeling, refinement and improvement over the best available template, and model ranking. Accurate modeling of the structurally divergent region such as the extracellular loops, disulfide bond formation affecting helix residue registry, and the helical shifts in the TM region seems to be particularly critical for accurately predicting the key ligand interactions in GPCRs, and this area is perhaps the most in need of technological development. Progress in GPCR modeling and docking will require further improvements in the current prediction methods to “add value” to the best available templates and generate models that will be more useful for applications in structure-based drug design.
We thank the participating predictors for making this assessment possible. We thank Mike Hanson, V.-P. Jaakola, Chris Roth, and Vadim Cherezov for help with the analysis and comments on the manuscript, and Katya Kadyshevskaya and Vadim Cherezov for figure preparation. We are grateful to the Goddard group for providing the script to calculate the binding site contactRMSD. We thank Angela Walker for data tracking and assistance with the manuscript, and Josh Kunken for IT help during the assessment. This work was supported in part by the NIH Roadmap grant P50 GM073197 (JCIMPT) and Protein Structure Initiative grant U54 GM074961 (ATCG3D).
Arthur Olson (Department of Molecular Biology, The Scripps Research Institute)
Wiktor Jurkowski and Arne Elofsson (Center of Biomembrane Research, Department of Biochemistry & Biophysics, Stockholm University)
Slawomir Filipek (Laboratory of Biomodelling, International Institute of Molecular and Cell Biology)
Irina Pogozheva and Andrei Lomize (Peptide Synthesis and Molecular Recognition Laboratory, University of Michigan)
Bernard Maigret (Orpailleur team, LORIA, Nancy University)
Jeremy Horst1, Ambrish Roy 2, Brady Bernard1, Shyamala Iyer1, Yang Zhang2, and Ram Samudrala1 (1 Computational Biology Group, University of Washington; 2 Department of Molecular Biosciences, Center for Bioinformatics, University of Kansas)
Osman Ugur Sezerman (Sabanici University)
Gregory V. Nikiforovich1 and Christina M. Taylor2 (1 MolLife Design LLC; 2 Department of Biochemistry and Molecular Biophysics, Washington University)
Stefano Costanzi (Laboratory of Biological Modeling, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health)
Y. Vorobjev2, N. Bakulina2, and V. Solovyev1,2 (1 Department of Computer Science, Royal Holloway, University of London; 2 Softberry Inc.)
Kazuhiko Kanou1, Daisuke Takaya1,2, Genki Terashi1, Mayuko Takeda-Shitaka1,2 and Hideaki Umeyama1,2 (1 School of Pharmacy, Kitasato University; 2 RIKEN Systems and Structural Biology Center)
William A. Goddard III, Youyong Li, Soo-Kyung Kim, Bartosz Trzaskowski, Ravinder Abrol, and Adam Griffith (Materials and Process Simulation Center, California Institute of Technology)
Vsevolod Katritch, Manuel Rueda and Ruben Abagyan (Molsoft LLC)
Ian Davis, Patrick Barth, David Baker (Department of Biochemistry, University of Washington)
Michael Feig (Department of Biochemistry & Molecular Biology, Michigan State University)
Michal Brylinski, Hongyi Zhou, Seung Yup Lee and Jeffrey Skolnick (Center for the Study of Systems Biology, Georgia Institute of Technology)
Liliana Ostopovici-Halip and Cristian Bologa (Division of Biocomputing, University of New Mexico)
Polo Lam and Ruben Abagyan (Department of Molecular Biology, The Scripps Research Institute)
Eric S. Dawson, Kristian Kaufmann, Nils Woetzel, and Jens Meiler (Center for Structural Biology, Vanderbilt University)
Feng Ding, Adrian Serohijos, Shuangye Yin, and Nikolay V. Dokholyan (Department of Biochemistry and Biophysics, University of North Carolina at Chapel Hill)
David Rodriguez and Hugo Gutiérrez-de-Teràn (Fundación Pública Galega de Medicina Xenómica. Complejo Hospitalario Universitario de Santiago, A Choupana, s/n. Santiago de Compostela)
Henri Xhaard (Center for Drug Research, Faculty of Pharmacy, University of Helsinki)