Essential Dynamics Ensemble Docking for Structure-Based GPCR Drug Discovery

The lack of biologically relevant protein structures can hinder rational design of small molecules to target G protein-coupled receptors (GPCRs). While ensemble docking using multiple models of the protein target is a promising technique for structure-based drug discovery, model clustering and selection still need further investigations to achieve both high accuracy and efficiency. In this work, we have developed an original ensemble docking approach, which identifies the most relevant conformations based on the essential dynamics of the protein pocket. This approach is applied to the study of small-molecule antagonists for the PAC1 receptor, a class B GPCR and a regulator of stress. As few as four representative PAC1 models are selected from simulations of a homology model and then used to screen three million compounds from the ZINC database and 23 experimentally validated compounds for PAC1 targeting. Our essential dynamics ensemble docking (EDED) approach can effectively reduce the number of false negatives in virtual screening and improve the accuracy to seek potent compounds. Given the cost and difficulties to determine membrane protein structures for all the relevant states, our methodology can be useful for future discovery of small molecules to target more other GPCRs, either with or without experimental structures.

Introduction

Many G protein-coupled receptors (GPCRs) are being investigated as important therapeutic targets, but the success rate of structure-based drug design (SBDD) for GPCRs remains to be further improved (Hauser et al., 2017; Wootten et al., 2018; Odoemelam et al., 2020). One of the primary challenges is that the three-dimensional (3D) structures of most GPCRs have not been fully determined. Even with latest breakthroughs in protein structure prediction like AlphaFold (Jumper et al., 2021), the available structures may not represent the conformational states needed for accurate SBDD. The receptor (ADCYAP1R1, hereafter referred to as PAC1R) of the pituitary adenylate cyclase-activating peptide (PACAP), an emerging therapeutic target for stress-related disorders (Hammack et al., 2009; Ressler et al., 2011; Roman et al., 2014; Missig et al., 2017; Liao et al., 2019a), is a good example. Currently, the full-length PAC1R structures in the Protein Data Bank (PDB) are short isoforms (Uniprot ID: P41586-3) (Kobayashi et al., 2020; Liang et al., 2020), but the structures of the most prevalent long isoforms—PAC1null (Uniprot ID: P41586) or PAC1hop (Uniprot ID: P41586-2) — are still unavailable (Liao et al., 2019b). All the published structures of PAC1R are complexed with peptide agonists and a heterotrimeric G protein complex (Figure 1), and thus do not represent the inactive conformations required for antagonist development. So far, it is thought that over 40% of GPCRs have more than one isoform (Marti-Solano et al., 2020), and each GPCR can adopt multiple conformational states which can be stabilized upon interactions with binding partners (Li et al., 2013; Vardy and Roth, 2013). For accurate SBDD, it is important to employ conformations of the most medically relevant isoform, as it is to this ensemble of 3D pocket structures that the drug must show affinity. Here, we used PAC1R as a model system and investigated how to improve modeling accuracy and to gain predictive power for SBDD with limited 3D structural information, using the method of Essential Dynamics Ensemble Docking (EDED). With the proof of principle, this method can be readily generalized to develop new therapeutic targets to target a wider range of GPCRs.

FIGURE 1

FIGURE 1. Cartoon illustrations of the PACAP-bound PAC1R model (PDBID: 6M1I, PAC1very short) and our homology model (template PDBID: 4L6R, PAC1null, simulation snapshot at 500 ns). The PAC1null isoform is more biomedically relevant than the very short isoform. The PACAP peptide is shown as a helix cartoon (pale green); the 21-amino acid ECD insert (see the sequence in Supplementary Figure S1) is shown as a flexible coil (purple). This study focused on docking to the peptide-binding pocket.

PAC1R and its endogenous peptide hormone PACAP play an important role in neural development, calcium homeostasis, glucose metabolism, circadian rhythm, thermoregulation, inflammation, feeding behavior, pain modulation, as well as stress, and related endocrine responses (Harmar et al., 2012; Bortolato et al., 2014; Culhane et al., 2015). For example, increased levels of PACAP in the blood have been reported in women diagnosed with post-traumatic stress disorder (Ressler et al., 2011), implicating chronic activation of the PAC1R in the disorder. Other studies (Boehr et al., 2009; Missig et al., 2017) have suggested that PAC1R activation mediates the adverse emotional consequences of chronic pain via downstream MAPK/ERK activation. Thus, these prior studies indicate that PAC1R antagonism, especially with small-molecule antagonists, represents a new strategy to treat stress, chronic pain, and related disorders (Ressler et al., 2011). Similar to other class B GPCRs, PAC1R possesses a heptahelical transmembrane domain (7TM) and an extracellular domain (ECD) (Odoemelam et al., 2020). Most of the neural and peripheral tissues known to date contain the PAC1null or PAC1hop isoforms that includes a 21-amino acid insert in the ECD (Figure 1), which is missing in available PAC1R structures in the PDB (May and Parsons, 2017). This ECD insert was found highly dynamic in our previous modeling studies (Liao et al., 2017; Liao et al., 2021), but its role in regulating PAC1R remains unknown. While PAC1R antagonists are being developed as potential treatments for stress-related disorders, the agonist-bound cryo-EM structures are not directly applicable to computational design or screening of PAC1R antagonists. GPCRs spontaneously adapt active and inactive signaling states, each of which are characterized by broad conformational ensembles. In the confirmational selection view, agonists and antagonists stabilize GPCR conformations of the active and inactive ensembles, respectively (Boehr et al., 2009; Abrol et al., 2013). It is now well accepted that to accurately design GPCR ligands as drug candidates, one should use active conformations for agonist design and inactive conformations for antagonist design. With the transition between active and inactive GPCR conformations occurring on the millisecond timescale (Vilardaga, 2010; Heyden et al., 2013; Weis and Kobilka, 2014; Scherer et al., 2015), it is computationally demanding to obtain the inactive PAC1R conformations from the agonist-bound cryo-EM structures via molecular dynamics (MD) simulations. Instead, we seek to use a homology model in this work and test with the EDED method.

Ensemble docking utilizes multiple receptor models for pocket sampling, obtained from clustering the conformations sampled by MD simulations for molecular docking, and displays noted improvement at identifying GPCR ligands when compared to docking against a single experimental structure (Lin et al., 2002; Huang and Zou, 2007; Amaro et al., 2018; Velazquez et al., 2018; Li et al., 2019; Acharya et al., 2020; Bhattarai et al., 2020; Chandak et al., 2020; Jukič et al., 2020; Patel et al., 2021; Li et al., 2022). EDED is distinct from prior ensemble docking approaches, mainly in clustering and selection of receptor models. Global root mean square deviation (RMSD) is convenient to cluster similar structures, but the highly dynamic extracellular and intracellular loops (ECLs and ICLs) of GPCRs can significantly compromise the otherwise good similarity between the 7TM structures. Thus, clustering based on global RMSD can generate many models that, while representative of global changes, are irrelevant to the intricate differences within the local binding pocket of the GPCR. This additional overhead ultimately lowers both the efficiency and accuracy of ensemble docking when using the global RMSD approach for clustering. EDED avoids this issue by focusing on both local similarity and essential dynamics of the binding pocket. Although computational power is more accessible than ever, streamlined workflows which expend computational resources only on worthwhile calculations are always desirable. Herein, we applied EDED to PAC1R with as few as four receptor models, whose results show a reduced false negative rate and a good correlation between the small molecule efficacy and the predicted score. Our results provide the evidence for initial success to develop small-molecule antagonists for PAC1R and pave the way for future structure-based GPCR drug discovery.

Results and Discussion

Inactive Conformations of PAC1null and Key Interactions With Small Molecules

Towards discovery of novel PAC1R antagonists, the inactive state conformational ensemble of PAC1R was estimated using an all-atom MD simulation of a ligand-bound PAC1R model from homology modeling (Figure 1). Our reference ligand is an analog of known PAC1 antagonists (Beebe et al., 2008) that were discovered previously using structure activity relationships. We created the antagonist-bound model by docking the reference compound into the PAC1R homology model. This complex model was simulated in the POPC membrane for 500 ns, and for the entire length of the simulation the ligand remained bound in roughly the starting conformation (Supplementary Figure S2). Other features like the closed ECD and straight transmembrane helix six (TM6), as well as short separation between TM3, TM5 and TM6, are consistent with a deactivated structure of a class B GPCR (Wu et al., 2020).

Despite the overall appearance of an inactivated receptor, there were critical changes within the orthosteric pocket during the MD simulations. Using EDED, four members of the inactive conformational ensemble (states S0, S1, S2, and S3 ordered by observed population) were extracted and reveal distinct conformations of the 7TM helices and different side chain orientations within the binding pocket (Figure 2). For one, bending of TM1 was observed to follow S0 < S2 < S3 < S1, where the most populated state (S0) was the most straightened helix. This correlated with local changes to residues Y161, L165, and S164 on TM1, and most significantly the stiffened TM1 in the S0 state enabled both π-stacking (with Y161) and hydrogen bonding (with S164) interactions. On the other hand, displacement of TM7 in the S1 state relative to S0 caused replacement of the hydrogen bond with S164 in the S0 state with a new hydrogen bond with S390. The interactions between the indole on the ligand and V368, L388, and E392 were modulated between the different receptor states with generally tighter interactions in the S1 and S3 states, in comparison with the S0 and S2 states. In addition, changes in TM5 affected the ability of K310 to form the stable interactions with the electron rich substituents on the ligand in states S0 and S3 which were diminished in the S1 and S2 states. Ultimately, this analysis reveals how EDED is able capture the subtle changes in pocket structure that are highly relevant for accurate modeling of ligand-receptor interactions when performing SBDD.

FIGURE 2

FIGURE 2. Four representative PAC1R conformations using in EDED reveal important changes in the binding pocket. The protein is shown in cartoon and the reference ligand is shown as sticks. The histogram of all trajectory frames projected onto the first two principal components of residues within the ligand-binding pocket of PAC1R. Black dots labelled with numbers from 0 to 3 are the representative structures (S0, S1, S2, and S3) determined by the K-means clustering algorithm.

Comparison of Docking to a Single Receptor Model and to the Conformational Representatives

Compared with docking to the ligand-free homology model and ligand-free cryo-EM structure, EDED significantly improved the identification of candidate compounds (Figure 3). The average binding score of the top 350 (approximately 2.5%) of compounds docked to the ligand-free homology model improved from −5.9 to −9.4 kcal/mol when docked against the ensemble. Likewise, it improved from an average of −5.8 to −9.4 kcal/mol when compared to the PACAP-bound model (PDBID: 6M1I). This gives an average 3.6 kcal/mol improvement in average docking score of the top selected compounds. Additionally, EDED identified six compounds predicted to bind to PAC1R with comparable binding score (−11.2 kcal/mol) as our reference ligand.

FIGURE 3

FIGURE 3. Violin plots of the docking score distribution of the top 350 compounds to different receptor models. The dash line shows the –9.0 kcal/mol cutoff used to prioritize compounds for synthesis.

To gain physical insight into the improvement of the docking scores, the binding pose(s) of the top compounds from both methodologies were examined. We have previously reported the key role of R199 in PACAP-induced activation of PAC1R (Liang et al., 2020). This is further corroborated by strong cation-pi interactions with the residue in our models. Interaction with R199 across all the ensemble conformations became a critical determining factor for which top ensemble docking compounds should be prioritized for synthesis and/or computational optimization. Examining the compounds which have ensemble docking scores close to or better than our reference ligand, this interaction is present for all six top scoring ligands in at least one of the docking poses. This is in contrast with the homology model and the PACAP-bound model where only relatively few of the top compounds from this methodology were able to engage in this key interaction. Also, of note are induced fit effects where the MD simulation of our reference ligand in the pocket may affect the binding pocket through subtle shifts in the backbone and the rotation of side chains. In the rigid receptor docking to the homology model, the 7TM helical bundle is closer together, defining a more compact orthosteric pocket. Thus, it is only accessible for small ligands to bind deep into the pocket below R199. In contrast, the conformations in the ensemble docking are more open, better allowing ligands to access the pocket. This can be seen by where most ligands found their best pose. Although both datasets were docked against a grid centered on R199, the ensemble docking results have the majority of top ligands below the residue, low in the pocket. When docked against the homology model, the top ligands are higher in the pocket at lowest in line with R199.

The new ligands examined within the orthosteric pocket showcased the ability of ensemble docking to provide integral confirmations omitted by static modelling, with the ensemble approach providing key ligand poses corresponding to interactions with new side chains revealed in the ensemble. Aside from R199, several key contacts were discovered from study of the top ligands bound to each receptor in the ensemble (Supplementary Figure S3). These contacts expand the understanding of the orthosteric pocket dynamics and can be exploited in small molecule rational design. In comparison with consistent interactions to the ligand-binding pocket of the homology model, these results suggests that EDED may reveal new crucial ligand-receptor interactions even from a rigid template.

A thermodynamically driven approach to scoring the binding poses of a given compound to multiple receptor structures was used to assess the binding affinity of the docked ligands. This approach quantitatively captures various physical phenomena that are often considered when computing overall docking scores: 1) the relative likelihood of the receptor obtaining the different conformations are explicitly included, and 2) the binding of the ligand to the receptor changes the energies of the complex differentially in the distinct conformations. Importantly, this model properly handles confounding cases that other approaches, such as a simple direct averaging of different docking scores, would not describe well. For instance, for any given ligand, a protein is hypothetically able to adopt an unlikely conformation ( $\mathrm{\Delta }G{\_}_{\mathrm{conf}1,i}\gg 0$, i.e., much higher relative free energy than the structure with lowest relative free energy) where the binding of the ligand to the protein could be quite favorable ($-\mathrm{\Delta }{G}_{bin{d}_i}$ approximately equal to 10 kT). Simply including this state in an average of docking scores would treat it as equivalently important as conformations that are far more relevant to the signaling states of the protein. Our approach avoids such errors, by including the energetics of the receptor confirmations, assuring that the overall energy of these rare states is indeed still relatively high and do not contribute significantly to the final score in Eq. 1. In sum, our docking score considers the difference in overall energies of the bound receptor conformations and is appropriate for comparison with a physical experiment that is unlikely to be able to distinguish between different bound conformations (Eq. 1).

$e^{\mathrm{\Delta }G/kT}=e^{\left(G_{bound}-G_{unbound}\right)/kT}=\frac{P_{unbound}}{P_{bound}}=\frac{{{\displaystyle \sum }}_{i=1}^n{P}_{cluster i}}{{{\displaystyle \sum }}_{i=1}^n{P}_{complex i}}(1)$

Where $G_{bound}$ and $G_{unbound}$ are the energies associated with the ligand being bound or unbound to any receptor conformation, respectively, $P_{bound}$ and $P_{unbound}$ are the total probabilities of the ligand being bound or unbound to any receptor conformation in the ensemble, respectively, $P_{cluster i}$ is the probability of a specific receptor conformation (calculated from the MD, see SI for more information), and $P_{\mathit{complex} i}$ is the probability of the ligand being bound to that specific ensemble conformation. We note that our model is still more appropriate than equal weighting for cases where one does not trust the relative energies of the different conformations obtained directly from the MD simulations. In such cases setting the $\mathrm{\Delta }G{\_}_{\mathrm{conf}1,i}$ to 0 for each conformation (i.e., each conformation is equally likely) reduces Eqs 1–2.

$\mathrm{\Delta }{\text{G}}_{bind,equal weighting}=\ln \left(\frac{n}{{{\displaystyle \sum }}_{i=1}^n{e}^{-\mathrm{\Delta }{G}_{bin{d}_{i/kT}}}}\right)kT(2)$

Clearly, Eq. 2 is not a simple weighted average of the different binding scores, however to our knowledge this analysis is lacking in the literature.

Evaluation of EDED Predictions

Additional to testing EDED with compounds from ZINC, we also tested 23 small-molecule compounds which were classified as strong, moderate, and weak antagonists in PAC1R activity assays (unpublished data from Prof. Victor May). The design of these small molecules was based on previously published work outlining the structure-activity relationship between small molecules and the PAC1 receptor (Beebe et al., 2008). Ligand-based virtual screening was then performed and yielded the 23 compounds which were experimentally tested. Docking each analog against all four conformations in the ensemble and scoring them as previous described (Eq. 1) shows modest correlation to experimental results (Figure 4). The strong experimental antagonist had the highest predicted binding affinities with an average −10.4 kcal/mol, while the moderate and weak antagonists both had worse predicted binding affinities −9.8 kcal/mol and −8.5 kcal/mol, respectively.

FIGURE 4

FIGURE 4. Ensemble weighted glide scores ($\text{Δ}$ G_bind) of 23 experimentally tested compounds. Compounds with strong, modest, and poor ERK inhibitive activity are depicted in green, blue, and red, respectively. Corresponding colored lines represent the average ensemble weighted glide score for that category. A cutoff of –9 kcal/mol was applied for predicted antagonists to be compared to their experimental results showing either strong or medium inhibition (active) or weak inhibition (inactive).

It is worth noting that our EDED method is best used to identify potential antagonists from a collection of compounds, but the dockings scores (like Glide SP, XP, and our EDED score) to estimate binding energies should be interpreted with caution (Elokely and Doerksen, 2013; Pantsar and Poso, 2018; Pinzi and Rastelli, 2019). While we successfully reduced the false negative rate (FNR) with EDED, there is still a high false positive rate (FPR). A delicate balance between ensemble size and the FPR has previously been reported, inspiring us to select a relatively small ensemble for analysis (Mohammadi et al., 2022). Our FPR is comparable to prior studies employing both ensemble and static methods for virtual screening (Ferreira et al., 2010; Deng et al., 2015; Hou et al., 2018). Additionally, the experimental assays provided here are a measure of antagonistic ability, and not binding affinity. As quantitative binding assays remain to be performed, it is possible some of the false positives (compounds with poor experimental results but high ensemble docking scores) bind tightly but are not effective antagonists, i.e., they do not stabilize the inactive conformations or prevent cognate ligand binding in other ways. With the extended view provided by EDED, we envision that the chance of obtaining a false negative prediction is likely reduced in our model when compared with static Glide docking. This added width within the sampled energy landscape (from the new side chain confirmations) allows our EDED method to achieve more accurate sampling of potential ligand-receptor interactions, thus increasing the chances of finding a hit compound otherwise overlooked in the static model. Overall, EDED displayed an accuracy of 57% in predicted binding affinity when compared to our experimental results, an increase when compared with Glide’s empirical scoring function (Adeshina et al., 2020). Combined with the overall low variance in EDED docking scores for the top 350 compounds analyzed (Figure 3), we believe our methodology represents a robust route for the recognition of small molecules with high receptor affinity.

Conclusion

In conclusion, we have developed and implemented EDED, an ensemble docking inspired methodology for SBDD. By focusing on the essential dynamics of the ligand binding pocket, our method is distinct from many prior studies that built receptor clusters solely based on the root mean square deviation (RMSD) of the entire protein backbone (Kufareva and Abagyan, 2012). Further, the use of clustering within this reduced dimensionality conformational space directly considers the local structural similarity of the ligand-binding pocket. We demonstrate that EDED captures the critical changes in the 3D structure of the binding pocket that are known to correlate strongly with binding affinity of ligands. Our approach is partially based on the assumption that differences in the binding pocket itself (as opposed to the protein as a whole) predominately give rise to the different binding poses and energies that are the goal of any ensemble docking workflow. Using the EDED derived representative structures, we screened a large dataset of compounds and successfully identified novel small molecule antagonists of the PAC1 receptor. However, EDED is not specific to a single GPCR and will likely accelerate the design of small molecule drugs that target other GPCRs with currently unknown conformational states.

Methods and Models

Receptor Model Preparation in EDED

One key idea of EDED is to obtain chemically relevant receptor models for docking. Instead of using the agonist-bound PAC1R structure, we generated a homology model of inactive PAC1R (with the canonical variant sequence, Uniprot ID: P41586) with a template of the glucagon receptor (PDBID: 4L6R, ∼40% similarity) (Boehr et al., 2009). This PAC1R model incorporated the inactive features of class B GPCRs such as a continuous helix along TM6 and a closed ECD. A small-molecule PAC1R antagonist, our reference ligand, was placed in the orthosteric pocket via molecular docking (Glide, Schrödinger Inc.). The complex model was later simulated to sample the inactive conformational ensemble.

Receptor Model Sampling in EDED

To sample inactive conformations for docking, the ligand-bound PAC1R model was simulated with the OPLS3 (Harder et al., 2016) force field in explicit SPC solvent in the NPT ensemble (300K, 1 atm, Martyna-Tuckerman Klein coupling scheme) using classical MD simulations. A POPC membrane was place around the 7TM using the Orientations of Proteins in Membrane (OPM) database (Lomize et al., 2012). The simulation was performed in the Maestro-Desmond program (Bowers et al., 2006) (GPU version 5.4) with a timestep of 2 fs, recording interval of 4.8 ps, and a total simulation time of 500 ns The Ewald technique was used for the electrostatic calculations. The van der Waals and short-range electrostatic interactions were cut off at 9 Å. Hydrogen atoms were constrained using the SHAKE algorithm. Two extended simulations were also examined to confirm the ligand poses and receptor confirmations. Once again, a POPC membrane was placed around the 7TM bundle using OPM. NAMD 2.11 was used as the simulation package for these replicates (Phillips et al., 2020). The CHARMM36 forcefield (Lee et al., 2016) was used with a TIP3 solvent model in a NPT ensemble (310 K, 1 atm) Force switching was utilized at the range of 10–12 Å to approximate the LJ interactions. Langevin piston/Nose-Hoover (Martyna et al., 1994; Feller et al., 1995) methods were utilized for the pressure control with a piston period of 50 fs and a decay time of 25 fs Langevin coupling of these simulations with a dampening coefficient of 1 ps⁻¹ was also utilized. Long range electrostatic interactions were modeled with the particle mesh Ewald method (Essmann et al., 1995). These simulations were run with a 2 fs timestep and combined for 350 ns of data. MD trajectories were analyzed using in-house Python and TCL scripts as well as Visual Molecular Dynamics (VMD). (Humphrey et al., 1996).

Receptor Ensemble Selection in EDED

We first aligned the 7TM of PAC1R (residues 156–405) to the homology model to reduce noise due to translational movement. Next, the coordinates of the centers of mass for any residue whose side chain was within 3 Å of any ligand atom in the static model were collected and parsed using in-house designed TCL and python scripts. A dimension reduction based on principal component analysis (PCA) was used to determine which collective motions (termed principal components, PCs) contributed most to variations in the overall conformations of the binding pocket. The first fifteen PCs (accounting for 90% of the cumulative variance) were clustered using a K-means clustering algorithm implemented by PyEmma (Liao et al., 2021). Based on inspection of the first two PCs (Figure 5), four cluster centers were identified. As these cluster centers are not precise frames within the trajectory but are instead points in the PC space, the cluster centers’ PC coordinates were approximately projected back to the original Cartesian coordinates. Frames from the trajectory which had PC values closest to the centers based on a RMSD measurement, were then selected as the ensemble docking receptor structures. This approach allowed a minimum of representative frames to capture the most variance of the binding pocket as opposed to other methodologies which often have many structures. Also, our physics-based approach is transferrable to other GPCRs and expanded clustering. In fact, our focus on the relevant receptor models likely requires less sampling in MD simulations and fewer clusters for subsequent docking, a practical advantage for large-scale screening.

FIGURE 5

FIGURE 5. Overview of computational workflow for development of PAC1R antagonists. Right Column: selection of input ligands from a structure database (in this example the ZINC15 (Sterling and Irwin, 2015) database). Custom filters were used to select raw structures with desirable properties (molecular weight, logP, etc.). These structures are then prepared using Schrödinger’s ligprep software program. Left column): the PAC1null homology model is constructed from the protein’s sequence, simulated for 500 ns, and the raw coordinates are analyzed. The representative structures are used in ensemble docking. Hit compounds are selected based on visual inspection of the results.

Docking and Scoring of Potential PAC1R Antagonists

Receptor grid models were generated using the three-dimensional structures selected as detailed above with R199 selected as the center of the docking box with an 18-Å cutoff. Docking was carried out using Schrödinger Virtual Screening Workflow (Friesner et al., 2006) (VSW) at three consecutive levels of precision, both for small molecules docked to the static homology model and to the conformation ensemble. Small molecules docked to our PAC1null ensembled were given an overall score, Ensemble $\mathrm{\Delta }$ G_bind, based on Eq. 3.

$\text{Ensemble }\mathrm{\Delta }{\text{G}}_{bind}=\ln \left(\frac{1+{{\displaystyle \sum }}_{i=2}^n{e}^{-\mathrm{\Delta }{G}_{con{f}_{1,i}}}}{{{\displaystyle \sum }}_{i=1}^n{e}^{-\mathrm{\Delta }{G}_{con{f}_{(1,i)}}-\mathrm{\Delta }{G}_{bin{d}_{(1,i)/kT}}}}\right)kT(3)$

In Eq. 3, $\mathrm{\Delta }{G}_{\mathrm{conf}1,i}$ is the difference in energy (in units of kT) between the lowest energy receptor conformation and each subsequent conformation calculated using the clustered trajectory, and $-\mathrm{\Delta }{G}_{bin{d}_i}$ is the corresponding Glide XP docking score to that same conformation. While $\mathrm{\Delta }{G}_{\mathrm{conf}1,i}$ is representative of the apo receptor free energy, it is worth noting that simulation data used to generate these confirmations included the ligand bound within the pocket.

Docking was carried out against compounds 1) pseudo-randomly selected from the ZINC15 (Sterling and Irwin, 2015) database, 2) as analogs of known antagonists to the static ligand-free homology model, the cryo-EM structure, and the conformational ensemble. In total, a small test set of 10,000 drug-like compounds were selected and download from the ZINC database and docked using Schrödinger’s VSW as described previously.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

KM, JR, SS, and JL contributed to the conception and design of the study. KM performed the model preparation, simulation setup, docking, and docking analysis. KM and NH performed the simulation analysis. JR derived the EDED scoring function. KM, JR, and JL wrote the first draft of the manuscript. KM, JR, NH, SS, and JL wrote sections of the manuscript. All authors contributed to the revision, read, and approved the submitted version.

Funding

The work was mainly supported by the NIH grant R01-GM129431 to JL. SS. was partially supported by the Army Research Office (Grant 71015-CH-YIP). STS was supported by the U.S. Army Research Office (Grant 71015-CH-YIP). Part of the computational facilities was also supported by an NSF CAREER award (Grant CHE-1848444 awarded to STS).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Drs. Victor May and Matthias Brewer (UVM) for sharing the experimental data and for helpful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.879212/full#supplementary-material

References

Abrol, R., Kim, S.-K., Bray, J. K., Trzaskowski, B., and Goddard, W. A. (2013). “Chapter Two - Conformational Ensemble View of G Protein-Coupled Receptors and the Effect of Mutations and Ligand Binding,” in Methods in Enzymology. Editor P. M. Conn (Academic Press), 520, 31–48. doi:10.1016/b978-0-12-391861-1.00002-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Acharya, A., Agarwal, R., Baker, M. B., Baudry, J., Bhowmik, D., Boehm, S., et al. (2020). Supercomputer-Based Ensemble Docking Drug Discovery Pipeline with Application to Covid-19. J. Chem. Inf. Model. 60 (12), 5832–5852. doi:10.1021/acs.jcim.0c01010

PubMed Abstract | CrossRef Full Text | Google Scholar

Adeshina, Y. O., Deeds, E. J., and Karanicolas, J. (2020). Machine Learning Classification Can Reduce False Positives in Structure-Based Virtual Screening. Proc. Natl. Acad. Sci. U.S.A. 117 (31), 18477–18488. doi:10.1073/pnas.2000585117

PubMed Abstract | CrossRef Full Text | Google Scholar

Amaro, R. E., Baudry, J., Chodera, J., Demir, Ö., McCammon, J. A., Miao, Y., et al. (2018). Ensemble Docking in Drug Discovery. Biophys. J. 114 (10), 2271–2278. doi:10.1016/j.bpj.2018.02.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Beebe, X., Darczak, D., Davis-Taber, R. A., Uchic, M. E., Scott, V. E., Jarvis, M. F., et al. (2008). Discovery and SAR of Hydrazide Antagonists of the Pituitary Adenylate Cyclase-Activating Polypeptide (PACAP) Receptor Type 1 (PAC1-R). Bioorg. Med. Chem. Lett. 18 (6), 2162–2166. doi:10.1016/j.bmcl.2008.01.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattarai, A., Wang, J., and Miao, Y. (2020). Retrospective Ensemble Docking of Allosteric Modulators in an Adenosine G-Protein-Coupled Receptor. Biochim. Biophys. Acta (BBA) - General Subj. 1864 (8), 129615. doi:10.1016/j.bbagen.2020.129615

PubMed Abstract | CrossRef Full Text | Google Scholar

Boehr, D. D., Nussinov, R., and Wright, P. E. (2009). The Role of Dynamic Conformational Ensembles in Biomolecular Recognition. Nat. Chem. Biol. 5 (11), 789–796. doi:10.1038/nchembio.232

PubMed Abstract | CrossRef Full Text | Google Scholar

Bortolato, A., Doré, A. S., Hollenstein, K., Tehan, B. G., Mason, J. S., and Marshall, F. H. (2014). Structure of Class B GPCRs: New Horizons for Drug Discovery. Br. J. Pharmacol. 171 (13), 3132–3145. doi:10.1111/bph.12689

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowers, K. J., Chow, D. E., Xu, H., Dror, R. O., Eastwood, M. P., Gregersen, B. A., et al. (2006). “Scalable Algorithms For Molecular Dynamics Simulations On Commodity Clusters, SC '06,” in Proceedings of the 2006 ACM/IEEE Conference on Supercomputing, 11-17 Nov. 2006, 43.

Google Scholar

Chandak, T., Mayginnes, J. P., Mayes, H., and Wong, C. F. (2020). Using Machine Learning to Improve Ensemble Docking for Drug Discovery. Proteins 88 (10), 1263–1270. doi:10.1002/prot.25899

PubMed Abstract | CrossRef Full Text | Google Scholar

Culhane, K. J., Liu, Y., Cai, Y., and Yan, E. C. (2015). Transmembrane Signal Transduction by Peptide Hormones via Family B G Protein-Coupled Receptors. Front. Pharmacol. 6, 264. doi:10.3389/fphar.2015.00264

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, N., Forli, S., He, P., Perryman, A., Wickstrom, L., Vijayan, R. S. K., et al. (2015). Distinguishing Binders from False Positives by Free Energy Calculations: Fragment Screening Against the Flap Site of HIV Protease. J. Phys. Chem. B 119 (3), 976–988. doi:10.1021/jp506376z

PubMed Abstract | CrossRef Full Text | Google Scholar

Elokely, K. M., and Doerksen, R. J. (2013). Docking Challenge: Protein Sampling and Molecular Docking Performance. J. Chem. Inf. Model. 53 (8), 1934–1945. doi:10.1021/ci400040d

PubMed Abstract | CrossRef Full Text | Google Scholar

Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 103 (19), 8577–8593. doi:10.1063/1.470117

CrossRef Full Text | Google Scholar

Feller, S. E., Zhang, Y., Pastor, R. W., and Brooks, B. R. (1995). Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 103 (11), 4613–4621. doi:10.1063/1.470648

CrossRef Full Text | Google Scholar

Ferreira, R. S., Simeonov, A., Jadhav, A., Eidam, O., Mott, B. T., Keiser, M. J., et al. (2010). Complementarity Between a Docking and a High-Throughput Screen in Discovering New Cruzain Inhibitors. J. Med. Chem. 53 (13), 4891–4905. doi:10.1021/jm100488w

PubMed Abstract | CrossRef Full Text | Google Scholar

Friesner, R. A., Murphy, R. B., Repasky, M. P., Frye, L. L., Greenwood, J. R., Halgren, T. A., et al. (2006). Extra Precision Glide: Docking and Scoring Incorporating a Model of Hydrophobic Enclosure for Protein−Ligand Complexes. J. Med. Chem. 49 (21), 6177–6196. doi:10.1021/jm051256o

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammack, S. E., Cheung, J., Rhodes, K. M., Schutz, K. C., Falls, W. A., Braas, K. M., et al. (2009). Chronic Stress Increases Pituitary Adenylate Cyclase-Activating Peptide (PACAP) and Brain-Derived Neurotrophic Factor (BDNF) mRNA Expression in the Bed Nucleus of the Stria Terminalis (BNST): Roles for PACAP in Anxiety-Like Behavior. Psychoneuroendocrinology 34 (6), 833–843. doi:10.1016/j.psyneuen.2008.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Harder, E., Damm, W., Maple, J., Wu, C., Reboul, M., Xiang, J. Y., et al. (2016). OPLS3: A Force Field Providing Broad Coverage of Drug-Like Small Molecules and Proteins. J. Chem. Theory Comput. 12 (1), 281–296. doi:10.1021/acs.jctc.5b00864

PubMed Abstract | CrossRef Full Text | Google Scholar

Harmar, A. J., Fahrenkrug, J., Gozes, I., Laburthe, M., May, V., Pisegna, J. R., et al. (2012). Pharmacology and Functions of Receptors for Vasoactive Intestinal Peptide and Pituitary Adenylate Cyclase-Activating Polypeptide: IUPHAR Review 1. Br. J. Pharmacol. 166 (1), 4–17. doi:10.1111/j.1476-5381.2012.01871.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hauser, A. S., Attwood, M. M., Rask-Andersen, M., Schiöth, H. B., and Gloriam, D. E. (2017). Trends in GPCR Drug Discovery: New Agents, Targets and Indications. Nat. Rev. Drug Discov. 16 (12), 829–842. doi:10.1038/nrd.2017.178

PubMed Abstract | CrossRef Full Text | Google Scholar

Heyden, M., Jardon-Valadez, H. E., Bondar, A.-N., and Tobias, D. J. (2013). GPCR Activation on the Microsecond Timescale in MD Simulations. Biophysical J. 104 (2Suppl. 1), 115a. doi:10.1016/j.bpj.2012.11.666

CrossRef Full Text | Google Scholar

Hou, X., Rooklin, D., Yang, D., Liang, X., Li, K., Lu, J., et al. (2018). Computational Strategy for Bound State Structure Prediction in Structure-Based Virtual Screening: A Case Study of Protein Tyrosine Phosphatase Receptor Type O Inhibitors. J. Chem. Inf. Model. 58 (11), 2331–2342. doi:10.1021/acs.jcim.8b00548

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, S. Y., and Zou, X. (2007). Ensemble Docking of Multiple Protein Structures: Considering Protein Structural Variations in Molecular Docking. Proteins 66 (2), 399–421. doi:10.1002/prot.21214

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: Visual Molecular Dynamics. J. Mol. Graph. 14 (1), 33–38. doi:10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Jukič, M., Janežič, D., and Bren, U. (2020). Ensemble Docking Coupled to Linear Interaction Energy Calculations for Identification of Coronavirus Main Protease (3CLpro) Non-Covalent Small-Molecule Inhibitors. Molecules 25 (24). doi:10.3390/molecules25245808

CrossRef Full Text | Google Scholar

Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., et al. (2021). Highly Accurate Protein Structure Prediction with AlphaFold. Nature 596 (7873), 583–589. doi:10.1038/s41586-021-03819-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, K., Shihoya, W., Nishizawa, T., Kadji, F. M. N., Aoki, J., Inoue, A., et al. (2020). Cryo-EM Structure of the Human PAC1 Receptor Coupled to an Engineered Heterotrimeric G Protein. Nat. Struct. Mol. Biol. 27 (3), 274–280. doi:10.1038/s41594-020-0386-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kufareva, I., and Abagyan, R. (2012). Methods of Protein Structure Comparison. Methods (Mol. Biol. Clift. N.J.) 857, 231–257. doi:10.1007/978-1-61779-588-6_10

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J., Cheng, X., Swails, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., et al. (2016). CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field. J. Chem. Theory Comput. 12 (1), 405–413. doi:10.1021/acs.jctc.5b00935

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, D., Jiang, K., Teng, D., Wu, Z., Li, W., Tang, Y., et al. (2022). Discovery of New Estrogen-Related Receptor α Agonists via a Combination Strategy Based on Shape Screening and Ensemble Docking. J. Chem. Inf. Model. 62 (3), 486–497. doi:10.1021/acs.jcim.1c00662

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Jonsson, A. L., Beuming, T., Shelley, J. C., and Voth, G. A. (2013). Ligand-Dependent Activation and Deactivation of the Human Adenosine A2A Receptor. J. Am. Chem. Soc. 135 (23), 8749–8759. doi:10.1021/ja404391q

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Zhang, X. X., Lin, Y. X., Xu, X. M., Li, L., and Yang, J. B. (2019). Virtual Screening Based on Ensemble Docking Targeting Wild-Type P53 for Anticancer Drug Discovery. Chem. Biodivers. 16 (7), e1900170. doi:10.1002/cbdv.201900170

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, Y.-L., Belousoff, M. J., Zhao, P., Koole, C., Fletcher, M. M., Truong, T. T., et al. (2020). Toward a Structural Understanding of Class B GPCR Peptide Binding and Activation. Mol. Cell 77 (3), 656–668.e5. doi:10.1016/j.molcel.2020.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, C., de Molliens, M. P., Schneebeli, S. T., Brewer, M., Song, G., Chatenet, D., et al. (2019). Targeting the PAC1 Receptor for Neurological and Metabolic Disorders. Curr. Top. Med. Chem. 19, 1399–1417. doi:10.2174/1568026619666190709092647

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, C., Remington, J. M., May, V., and Li, J. (2021). Molecular Basis of Class B GPCR Selectivity for the Neuropeptides PACAP and VIP. Front. Mol. Biosci. 8, 644644. doi:10.3389/fmolb.2021.644644

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, C., May, V., and Li, J. (2019). PAC1 Receptors: Shapeshifters in Motion. J. Mol. Neurosci. 68 (3), 331–339. doi:10.1007/s12031-018-1132-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, C., Zhao, X., Brewer, M., May, V., and Li, J. (2017). Conformational Transitions of the Pituitary Adenylate Cyclase-Activating Polypeptide Receptor, a Human Class B GPCR. Sci. Rep. 7 (1), 5427. doi:10.1038/s41598-017-05815-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J.-H., Perryman, A. L., Schames, J. R., and McCammon, J. A. (2002). Computational Drug Design Accommodating Receptor Flexibility: The Relaxed Complex Scheme. J. Am. Chem. Soc. 124 (20), 5632–5633. doi:10.1021/ja0260162

PubMed Abstract | CrossRef Full Text | Google Scholar

Lomize, M. A., Pogozheva, I. D., Joo, H., Mosberg, H. I., and Lomize, A. L. (2012). OPM Database and PPM Web Server: Resources for Positioning of Proteins in Membranes. Nucleic Acids Res. 40, D370–D376. doi:10.1093/nar/gkr703

PubMed Abstract | CrossRef Full Text | Google Scholar

Marti-Solano, M., Crilly, S. E., Malinverni, D., Munk, C., Harris, M., Pearce, A., et al. (2020). Combinatorial Expression of GPCR Isoforms Affects Signalling and Drug Responses. Nature 587 (7835), 650–656. doi:10.1038/s41586-020-2888-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Martyna, G. J., Tobias, D. J., and Klein, M. L. (1994). Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 101 (5), 4177–4189. doi:10.1063/1.467468

CrossRef Full Text | Google Scholar

May, V., and Parsons, R. L. (2017). G Protein-Coupled Receptor Endosomal Signaling and Regulation of Neuronal Excitability and Stress Responses: Signaling Options and Lessons from the PAC1 Receptor. J. Cell. Physiol. 232 (4), 698–706. doi:10.1002/jcp.25615

PubMed Abstract | CrossRef Full Text | Google Scholar

Missig, G., Mei, L., Vizzard, M. A., Braas, K. M., Waschek, J. A., Ressler, K. J., et al. (2017). Parabrachial Pituitary Adenylate Cyclase-Activating Polypeptide Activation of Amygdala Endosomal Extracellular Signal-Regulated Kinase Signaling Regulates the Emotional Component of Pain. Biol. Psychiatry 81 (8), 671–682. doi:10.1016/j.biopsych.2016.08.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohammadi, S., Narimani, Z., Ashouri, M., Firouzi, R., and Karimi‐Jafari, M. H. (2022). Ensemble Learning from Ensemble Docking: Revisiting the Optimum Ensemble Size Problem. Sci. Rep. 12 (1), 410. doi:10.1038/s41598-021-04448-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Odoemelam, C. S., Percival, B., Wallis, H., Chang, M.-W., Ahmad, Z., Scholey, D., et al. (2020). G-Protein Coupled Receptors: Structure and Function in Drug Discovery. RSC Adv. 10 (60), 36337–36348. doi:10.1039/d0ra08003a

PubMed Abstract | CrossRef Full Text | Google Scholar

Pantsar, T., and Poso, A. (2018). Binding Affinity via Docking: Fact and Fiction. Molecules 23 (8), 1899. doi:10.3390/molecules23081899

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, D., Athar, M., and Jha, P. C. (2021). Exploring Ruthenium‐Based Organometallic Inhibitors Against Plasmodium Falciparum Calcium Dependent Kinase 2 (PfCDPK2): A Combined Ensemble Docking, QM/MM and Molecular Dynamics Study. ChemistrySelect 6 (32), 8189–8199. doi:10.1002/slct.202101801

CrossRef Full Text | Google Scholar

Phillips, J. C., Hardy, D. J., Maia, J. D. C., Stone, J. E., Ribeiro, J. V., Bernardi, R. C., et al. (2020). Scalable Molecular Dynamics on CPU and GPU Architectures with NAMD. J. Chem. Phys. 153 (4), 044130. doi:10.1063/5.0014475

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinzi, L., and Rastelli, G. (2019). Molecular Docking: Shifting Paradigms in Drug Discovery. Int. J. Mol. Sci. 20 (18), 4331. doi:10.3390/ijms20184331

PubMed Abstract | CrossRef Full Text | Google Scholar

Ressler, K. J., Mercer, K. B., Bradley, B., Jovanovic, T., Mahan, A., Kerley, K., et al. (2011). Post-Traumatic Stress Disorder Is Associated with PACAP and the PAC1 Receptor. Nature 470 (7335), 492–497. doi:10.1038/nature09856

PubMed Abstract | CrossRef Full Text | Google Scholar

Roman, C. W., Lezak, K. R., Hartsock, M. J., Falls, W. A., Braas, K. M., Howard, A. B., et al. (2014). PAC1 Receptor Antagonism in the Bed Nucleus of the Stria Terminalis (BNST) Attenuates the Endocrine and Behavioral Consequences of Chronic Stress. Psychoneuroendocrinology 47, 151–165. doi:10.1016/j.psyneuen.2014.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Scherer, M. K., Trendelkamp-Schroer, B., Paul, F., Pérez-Hernández, G., Hoffmann, M., Plattner, N., et al. (2015). PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. J. Chem. Theory Comput. 11 (11), 5525–5542. doi:10.1021/acs.jctc.5b00743

PubMed Abstract | CrossRef Full Text | Google Scholar

Sterling, T., and Irwin, J. J. (2015). ZINC 15 - Ligand Discovery for Everyone. J. Chem. Inf. Model. 55 (11), 2324–2337. doi:10.1021/acs.jcim.5b00559

PubMed Abstract | CrossRef Full Text | Google Scholar

Vardy, E., and Roth, B. L. (2013). Conformational Ensembles in GPCR Activation. Cell 152 (3), 385–386. doi:10.1016/j.cell.2013.01.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Velazquez, H. A., Riccardi, D., Xiao, Z., Quarles, L. D., Yates, C. R., Baudry, J., et al. (2018). Ensemble Docking to Difficult Targets in Early-Stage Drug Discovery: Methodology and Application to Fibroblast Growth Factor 23. Chem. Biol. Drug Des. 91 (2), 491–504. doi:10.1111/cbdd.13110

PubMed Abstract | CrossRef Full Text | Google Scholar

Vilardaga, J.-P. (2010). Theme and Variations on Kinetics of GPCR Activation/Deactivation. J. Recept. Signal Transduct. 30 (5), 304–312. doi:10.3109/10799893.2010.509728

PubMed Abstract | CrossRef Full Text | Google Scholar

Weis, W. I., and Kobilka, B. K. (2014). The Molecular Basis of G Protein-Coupled Receptor Activation. Annu. Rev. Biochem. 87 (1), 897–919. doi:10.1146/annurev-biochem-060614-033910

PubMed Abstract | CrossRef Full Text | Google Scholar

Wootten, D., Christopoulos, A., Marti-Solano, M., Babu, M. M., and Sexton, P. M. (2018). Mechanisms of Signalling and Biased Agonism in G Protein-Coupled Receptors. Nat. Rev. Mol. Cell Biol. 19 (10), 638–653. doi:10.1038/s41580-018-0049-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Yang, L., Hang, K., Laursen, M., Wu, L., Han, G. W., et al. (2020). Full-Length Human GLP-1 Receptor Structure Without Orthosteric Ligands. Nat. Commun. 11 (1), 1272. doi:10.1038/s41467-020-14934-5

PubMed Abstract | CrossRef Full Text | Google Scholar