Abstract
G-protein-coupled receptors (GPCRs) mediate multiple signaling pathways in the cell, depending on the agonist that activates the receptor and multiple cellular factors. Agonists that show higher potency to specific signaling pathways over others are known as “biased agonists” and have been shown to have better therapeutic index. Although biased agonists are desirable, their design poses several challenges to date. The number of assays to identify biased agonists seems expensive and tedious. Therefore, computational methods that can reliably calculate the possible bias of various ligands ahead of experiments and provide guidance, will be both cost and time effective. In this work, using the mechanism of allosteric communication from the extracellular region to the intracellular transducer protein coupling region in GPCRs, we have developed a computational method to calculate ligand bias ahead of experiments. We have validated the method for several β-arrestin–biased agonists in β2-adrenergic receptor (β2AR), serotonin receptors 5-HT1B and 5-HT2B and for G-protein–biased agonists in the κ-opioid receptor. Using this computational method, we also performed a blind prediction followed by experimental testing and showed that the agonist carmoterol is β-arrestin–biased in β2AR. Additionally, we have identified amino acid residues in the biased agonist binding site in both β2AR and κ-opioid receptors that are involved in potentiating the ligand bias. We call these residues functional hotspots, and they can be used to derive pharmacophores to design biased agonists in GPCRs.
Introduction
G-protein-coupled receptors (GPCRs) are a superfamily of signaling proteins of critical importance to cellular function, and comprise the largest group of drug targets. The binding of endogenous agonists to their target GPCRs results in coupling to a variety of intracellular transducer proteins, including various heterotrimeric G-proteins and β-arrestins, that mediate different downstream signaling pathways. Certain GPCR agonists elicit differential responses to different signaling pathways and are known as biased agonists, and the phenomenon is called functional selectivity or biased signaling (Urban et al., 2007; Rajagopal et al., 2010; Rasmussen et al., 2011; Wisler et al., 2014).
The ability of biased agonists to cause differential signaling resulting in more focused cellular changes could translate to lowering the unfavorable side effects of a drug. This property of biased agonists can markedly improve their therapeutic index when compared with unbiased or balanced agonists (Mailman, 2007; Neve, 2009; Violin et al., 2014). Design of biased agonists remains a daunting challenge due to the complexity of the cellular and structural factors that contribute to the bias. Therefore, an understanding of the structural underpinnings of biased signaling of both G-protein bias as well as β-arrestin bias of GPCRs is needed to design biased agonists for GPCRs.
Linking GPCR Conformation to Cellular Function.
Several factors contribute to biased signaling by an agonist–GPCR pair in cells: 1) structural features of the agonist–GPCR–G-protein or agonist–GPCR–β-arrestin complex (Shukla et al., 2014) and 2) cell-specific and tissue-specific factors (Zhou and Bohn, 2014). Designing biased agonists is a challenge due to the lack of structural information on how an agonist–GPCR pair selectively couples to specific G-proteins or β-arrestins. Moreover, because GPCRs are dynamic, with various possible active-state conformations, the structural ensemble of the agonist–GPCR–transducer complex is one of the key determinants driving the biased signaling (Shukla et al., 2014). The selectivity of GPCRs to trigger distinct intracellular signaling profiles is very ligand dependent. Spectroscopic studies have now established that ligands differentially stabilize a range of GPCR conformations (Yao et al., 2006; Kahsai et al., 2011; Liu et al., 2012; Nygaard et al., 2013). Therefore, knowledge about how the agonist affects the agonist–GPCR–G-protein structural ensemble will significantly aid the design of biased agonists.
Challenges in Biased Ligand Identification and Design.
Given the huge challenges in crystallography and NMR in determining structural ensembles of agonist-GPCR-effector complexes, computational methods can greatly aid in delineating the structural factors that determine agonist bias and its outcome. Despite a wealth of biochemical and biophysical studies on inactive conformations, there is paucity of structural information on active conformations of GPCRs. More importantly the dynamics of the agonist-GPCR-effector complex is one of the major factors leading to a conformational ensemble, which governs the functional selectivity of biased signaling.
A better understanding of the structural basis of G protein/β-arrestin selection will provide for a more precise targeting of GPCRs with pharmaceuticals and also inform ongoing structure-based drug discovery efforts. Hence, our goal in this work is to identify structural elements in the GPCR that confer functional selectivity either to the G-protein signaling pathway or to the β-arrestin signaling pathway. Our unanswered question has been, how do the various biased agonists influence the role of these structural determinants in biasing the signaling to effect β-arrestin signaling or G-protein signaling? Here we have developed a computational method to calculate the ligand bias ahead of experiments, which would aid the design of biased agonists. The two-part goal in this work is: 1) to develop a method to compute the ligand bias ahead of experiments and 2) to identify the residues in the agonist binding site that potentiate the ligand mediated bias. We define these residues as functional hotspots. To this end we chose systems that cover a wider base of bias systems.
The human β2-adrenergic receptor (β2AR) and angiotensin receptor AT1R have been studied extensively for discovery of β-arrestin–biased ligands (Violin et al., 2014; Wisler et al., 2014). Therefore, robust experimental data on the quantitative ligand bias factor calculated using the Black-Leff model (Black and Leff, 1983) are available for these receptors. We chose β2AR for our study because the three-dimensional (3D) crystal structure of the fully active Gs-protein bound state is available and this constitutes a robust test case for the computational method. The experimental values of the ligand bias factor have been calculated using cAMP secondary messenger assays for assessing the G-protein coupling potency for the β2AR and using Tango or other β-arrestin recruitment assays for calculating the ligand potency toward β-arrestin–mediated signaling pathways (Rajagopal et al., 2011; Weiss et al., 2013; Zhou et al., 2013).
There has been a surge in efforts to develop G-protein–biased opioid receptor agonists that serve as better analgesics with reduced side effects coming from receptor desensitization (Soergel et al., 2014; Luttrell et al., 2015; Stahl et al., 2015). We tested three agonists, two of which showed different extents of G-protein bias for the κ-opioid receptor (κOR) (Zhou et al., 2013). We chose to study the κOR rather than the μ-opioid receptor because the active state structure of κOR is not available and we wanted to test our computational method using a homology model of the receptor that makes the testing more realistic.
The serotonin receptors 5-hydroxytryptamine 1B (5-HT1B) and 5-HT2B have been crystallized with the same agonist ergotamine, yet ergotamine shows β-arrestin bias in 5-HT2B and not in 5-HT1B (Wacker et al., 2013). We chose to study the system to find out which residues in 5-HT2B cause the bias in ergotamine. This example would provide a validation for the same ligand showing different biases in two related receptors. Another similar example is the bias behavior of epinephrine in the triple mutant T68F, Y132G, Y219A β2AR, denoted here as β2ARTYY. This mutant was designed using evolutionary trace analysis to identify residues that could cause bias in β2ARTYY (Shenoy et al., 2006). Thus, choosing the variety of different systems described here provides a test for the robustness of our computational method.
To summarize, we have studied β2AR and its biased mutant β2ARTYY with several agonists, 5HT1B and 5HT2B with the same agonist, and κOR with three agonists including two G-protein–biased agonists. These systems are also listed in Supplemental Table 1 of the Supplemental Information for convenience.
Materials and Methods
A New Computational Method for Calculating Ligand Bias
The experimentally measured ligand bias factor depends on the ligand affinity (KA−1) and its efficacy (τ) toward a specific signaling pathway. The bias factor is calculated with respect to a reference ligand using the Black and Leff operational model (Kenakin et al., 2012). Because macroscopic properties such as ligand potency and coupling efficacy are still intractable computationally, we developed a method to calculate the ligand bias from atomic-level properties of the structural ensemble of the agonist-GPCR complex that potentiates the agonist bias.
Agonist and G-protein binding to GPCRs leads to G-protein activation. Such activation at a distance is potentiated by the allosteric communication between the ligand binding site and the G-protein coupling site in the fully active state of GPCRs. We recently developed a computational method to calculate the strength of the allosteric communication pipelines from the extracellular region of the GPCR to the intracellular regions where the G-protein or β-arrestin couples to the receptor. In this work we hypothesize that the strength of the allosteric communication pipelines from the extracellular (EC) region passing through the residues in the ligand-binding site (ligand shown in blue spheres in Fig. 1) to the residues in the G-protein and β-arrestin coupling interface in the GPCRs is related to the ligand potency and its coupling efficacy.
The model used for calculating the ligand bias in GPCRs. The allosteric communication pipelines starting from the residues in the extracellular loop region and passing through the residues in the agonist binding site (agonist shown as blue spheres) and terminating in the residues that couple to the G-protein (GPI shown as pink surfaces) are shown as pink sticks. The allosteric communication pipelines that connect to the β-arrestin interface (BAI shown as green surfaces) are shown in green sticks. The ratio of the strength of these two pipelines for a test ligand with respect to a reference ligand is defined as the ligand bias.
Figure 1 shows the scheme used here for calculating the ligand bias in GPCRs: the allosteric communication pipelines that connect to the G-protein interface are shown as pink sticks, and those to the β-arrestin interface are shown as green sticks. We defined the ligand bias as the ratio of the strength of the allosteric communication pipelines to the G-protein pathway to that of the β-arrestin signaling pathway for a test ligand, to that of the same quantity for the reference ligand. We used a reference ligand in these calculations so that these quantities can be directly compared with the experimental bias factor.
Computational Methods
The computational methods for calculating the allosteric communication pipelines in GPCRs involve the following steps: 1) preparing the GPCR structure with agonists for the molecular dynamics (MD) simulations and 2) performing MD simulations in explicit membrane bilayer and water. The list of receptors used in this work is shown in Supplemental Table 1. The two-dimensional representation of the agonists used in this study is shown in Supplemental Fig. 1. The details of the receptor structure preparation, docking of the agonists, and MD simulations are provided in the Supplemental Information file. We performed brute force MD simulations with no constraints for wild-type β2AR with several ligands listed in Supplemental Table 1, β2ARTYY mutant with epinephrine, serotonin receptors 5-HT1B and 5-HT2B with ergotamine, and homology model of κOR with three different agonists bound.
Allosteric Pipeline Analysis
The MD simulation trajectories were used to calculate the allosteric communication pipelines for each agonist-GPCR pair using the Allosteer computational method (Bhattacharya and Vaidehi, 2014a,b; Bhattacharya et al., 2016; Vaidehi and Bhattacharya, 2016), which we have developed previously elsewhere. The details of the Allosteer computational method for calculating the allosteric communication pipelines in GPCRs are described by Bhattacharya and Vaidehi (2014a,b). Here we provide a brief description of the method as used in this work.
Using the MD trajectories, we first calculated the correlated movement between residues on the EC surface of the receptor and the residues in G-protein or β-arrestin coupling interface. Here we have used the trajectories from the five 200-nanosecond MD simulations collected as a total of 1-microsecond of MD simulations to calculate the mutual information in torsion angle space between all pairs of residues in each receptor-ligand complex. Using the mutual information, we then calculated the shortest pathway from the EC region to the G-protein or β-arrestin coupling surfaces that maximizes the correlated movement. This method is called Allosteer (Bhattacharya and Vaidehi, 2014a).
Using the Allosteer method we then calculated the allosteric communication pipelines starting from the residues in the EC loops (defined in the section Residue Group Definitions) passing through the residues in the agonist-binding site to the residues on the β-arrestin or G-protein interacting interfaces. The list of residues in each of these regions is given in Supplemental Table 2. For each MD run, the top 10% of allosteric pathways ranked by the total mutual information were used for further calculation of the ligand bias. The strength of an allosteric communication pipeline is the number of overlapping allosteric communication pathways contained in the pipeline (Bhattacharya and Vaidehi, 2014a; Vaidehi and Bhattacharya, 2016).
Residue Groups Definitions
Here we list the residues that define the various regions in the GPCRs. These residues were included in the calculation of the allosteric communication pipelines.
Ligand Binding Site Residues.
The residues located within 5 Å of the agonist-binding site in the GPCR and contacting the agonist in more than 40% of the snapshots from the MD simulation trajectories were chosen as the binding site residues for that agonist. The super set of binding site residues for all agonists to a particular receptor were listed by combining the residues list for all the agonists to the receptor studied. This superset of residues was used to calculate the allosteric communication pipeline from the EC loops going through the ligand binding site to the G-protein or β-arrestin coupling interface. The list of residues in the ligand binding site used in this study is shown in Supplemental Table 2.
β-Arrestin Interface Residues and G-Protein Interface Residues.
The list of residues in the GPCR interacting with the G-protein was obtained from Supplementary Table 1 of the crystal structure paper of β2AR-Gs protein complex (Rasmussen et al., 2011). The Protein Data Bank (PDB) ID of this crystal structure is 3SN6. The list of residues interacting with the β-arrestin interface was obtained by aligning the β2AR structure from PDB ID: 3SN6 to the crystal structure of rhodopsin bound to visual arrestin (PDB ID: 5DGY) (Zhou et al., 2016) and identifying the residues in the receptor that were within 5 Å of visual arrestin in the aligned structure. The κOR structure was aligned to that of β2AR, and the β-arrestin interface (BAI) and G-protein interface (GPI) residues were translated for the κOR based on the alignment of residues and corresponding Ballesteros-Weinstein numbers. The list of residues in the GPI and BAI used in this study are shown in Supplemental Table 2.
EC Region Residues.
Residues comprising the EC loops of the receptor in addition to those residues forming two turns in the EC region of each transmembrane helix formed the set of EC residues. The EC loop residues included in the calculation of allosteric communication pipelines for the GPCRs studied here are provided in Supplemental Table 2.
Definition of the Computational Ligand Bias
Using the MD trajectories, we calculated the allosteric communication pipelines of every agonist-GPCR complex to the G-protein and β-arrestin interface. We then calculated the ratio of the strength of the allosteric communication pipeline (PR) for the β-arrestin coupling interface to that of the G-protein coupling interface residues. We did the same calculation for all agonists, including the reference agonist, which is usually a balanced agonist, used in the experimental assays for the bias factor calculations. The computational ligand bias values is defined as:(1)where PRlig is the ratio of the strength of the allosteric communication pipeline for a given agonist-receptor pair for the β-arrestin to the G-protein coupling interfaces, and PRref is the same for a reference agonist.
Testing the Robustness of the Method to Calculate the Ligand Bias
The computational ligand bias shown in eq. 1 is a ratio that is obtained by calculating the strength of the allosteric communication pipelines to the β-arrestin and the G-protein coupling interfaces. The error in the calculation of the strength of the allosteric pipelines is dependent on: 1) the number of snapshots used from the MD trajectories and 2) the simulation time length of each MD simulation. The dependence on the number of snapshots stored in the MD simulations is because the state of the system has to be captured often enough to enable the recording of all significant, meaningful correlations. However, recording snapshots at extremely frequent intervals would increase the probability of registering spurious correlations and dependencies, which in turn will introduce noise into the bias calculations.
The second factor, the length of the MD simulations, is important because longer MD simulations capture the collective correlated motions of torsion angles in proteins. However, a tradeoff has to be made between the cost of computational resources and accuracy of ligand bias calculated here. Therefore, we examined the optimum number of snapshots to be stored and length of MD simulations to reduce the significant errors in the ligand bias calculations. Thus, we calculated the variation in the S.E. in the pipeline strength as a function of these two variables. The S.E. in these calculations was analyzed with the aim of choosing optimal values for both (Supplemental Figs. 2 and 3). Upon analysis of the S.E. values, storing snapshots from MD simulations at every 20-picosecond time step in a 200-nanosecond MD simulation (10,000 snapshots) was chosen as the optimum frame rate, and 200 nanoseconds was chosen as the simulation time length for the calculations. For more details on these tests we refer the reader to the Supplemental Information.
Experimental Methods.
β-Arrestin Recruitment Assay and cAMP Detection Assay.
To test the predictivity of the computational methods, three compounds (isoproterenol, formoterol, and carmoterol) were tested in a β-arrestin assay as well as a cAMP assay on human β2-adrenoceptors. Isoproterenol and formoterol were purchased from Sigma-Aldrich Chemie, (Hoxter, Germany), and carmoterol was synthesized using the protocol described in the US Patent US2010/0113790 A1 (Kankan et al., 2010). For β-arrestin assay we chose the PathHunter β-Arrestin Assay system (DiscoveRx, Fremont, CA); for cAMP, the AlphaScreen (PerkinElmer Life and Analytical Sciences, Waltham, MA) technology was employed.
β-Arrestin Recruitment Assay.
The PathHunter β-Arrestin Assay system (DiscoveRx) was used to determine the β-arrestin recruitment to human β2-adrenoceptors. Chinese hamster ovary K1 (CHO-K1) cells stably expressing β2-adrenoceptors fused to ProLink and the β-arrestin-enzyme acceptor were seeded onto opaque 384-well plates at 5000 cells per well in Dulbecco’s modified Eagle’s medium/Ham’s F-12 (1:1) without phenol red supplemented with 1% heat-inactivated fetal bovine serum. The cells were maintained at 37°C in a 5% CO2, humidified atmosphere overnight. The next day the cells were stimulated in growth medium with different concentrations of agonists at 37°C for 60, 30, 15, and 5 minutes. Arrestin recruitment resulting in β-galactosidase enzyme complementation was detected by use of the PathHunter Flash detection kit (DiscoveRx). Luminescence was read on an Envision plate reader (PerkinElmer Life and Analytical Sciences). We performed 13 experiments for isoproterenol, 9 for formoterol, and 4 for carmoterol.
cAMP Detection Assay.
Changes in intracellular cAMP levels were determined by AlphaScreen technology (PerkinElmer Life and Analytical Sciences) in CHO-K1 cells stably expressing β2-adrenoceptors following the manufacturer’s instructions. Cells in suspension (15,000 cells per well) were stimulated in opaque 384-well plates with respective agonists at different concentrations in Hank’s buffered saline solution supplemented with 5 mM HEPES, 0.1% bovine serum albumin, and 500 mM 3-isobutyl-1-methylxanthine for 30 minutes at room temperature. The cells were lysed by using Alphascreen reagents. After 2 hours, the plates were read on an Envision plate reader (PerkinElmer Life and Analytical Sciences). The concentration of cAMP in the samples was calculated from a standard curve. We performed 12 experiments each for isoproterenol, formoterol, and carmoterol.
Black and Leff Model Used for Calculating the Bias Factor
The concentration–response data for each agonist at the G protein and β-arrestin signaling assays were fitted to the Black-Leff model of agonism (Black and Leff, 1983):(2)where [A] is the agonist concentration, τ is equal to RT/KE (RT: receptor density; KE: intrinsic agonist efficacy), n is the transducer slope, and Em is the maximal response of the system. The transduction coefficient for each agonist for a given pathway was calculated as log(τ/KA) (Kenakin and Christopoulos, 2013). For fitting the data, the Black-Leff equation was recast to a different form according to Tschammer et al. (2011), and log(τ/KA) was directly obtained from the fit. Em was estimated as the maximum value of the signaling response for a given assay.
Setting the value of n = 1 gave a good fit for all the dose–response data. In practice, n was allowed to vary within a very narrow range (0.9–1.2) to account for statistical variability. The fitting was performed using the Genetic Algorithm (GA) module in MATLAB (MathWorks, Natick, MA) (the operational model fitting of GraphPad Prism did not always find a solution for all data sets). Instead of starting from a single initial guess of the solution, 10,000 initial guesses were randomly generated within a prescribed range. The provided range for KA was 10−15 to 1, while for log(τ/KA) it was 0–15 (the range for n was stated earlier). Using the different initial guesses, the algorithm converged to a solution within the provided tolerance limit of 10−8.
For each agonist and a given pathway, Δlog(τ/KA) was calculated as log(τ/KA) − log(τ/KA)ref, where log(τ/KA)ref is the value for the reference agonist isoproterenol. ΔΔlog(τ/KA) was then calculated as Δlog(τ/KA)arrestin − Δlog(τ/KA)G protein. Thus, the bias factor is given by 10ΔΔlog(τ/KA).
To estimate the variability in the calculated log(τ/KA) values, we calculated Sij2 for agonist i and pathway j as:(3)where k denotes the number of experiments and y corresponds to the log(τ/KA) for the agonist/pathway pair. The total variability of the estimates Spooled is given by:
(4)where dferror is the degree of freedom given by:
(5)The 95% confidence levels for the calculated log(τ/KA) values are given by:
(6)where S.E. is given by:
(7)and T corresponds to a two-tailed t test with 95% confidence. For calculating the confidence levels for ΔΔlog(τ/KA), we used the same method (eq. 6), except the S.E. is given by:
(8)where n refers to the number of experiments for the given pathway, and nref is the number of experiments for the reference ligand.
Results
Calculation of Ligand Bias Factor Using Black-Leff Operational Model for Carmoterol and Formoterol.
We measured the ligand bias factor of carmoterol and formoterol using isoproterenol as the reference ligand. Using the concentration-dependent curves measured using the cAMP assay for G-protein activation and PathHunter assay for β-arrestin recruitment, as detailed in the Materials and Methods section, we calculated the ligand bias factor for carmoterol and formoterol with reference to isoproterenol in β2AR. We used the Black-Leff model to obtain bias factors as described in the Materials and Methods section. The calculated bias factors are plotted in Fig. 2A and are also shown in Supplemental Table 3. Our results show that carmoterol has a 2-fold higher bias factor toward the β-arrestin signaling pathway compared with formoterol in β2AR.
Comparison of computational ligand bias (y axis) to the experimental ligand bias factors (x axis) for (A) the β-arrestin–biased agonists in β2AR. epi, epinephrine; iso, isoproterenol; sal, salbutamol; form, formoterol; BI, BI-167107; carm, carmoterol. (B) the G-protein–biased agonists in κOR. u69, u69593. (A) The bias values obtained using epinephrine as the reference agonist are shown as green triangles, and those obtained using isoproterenol as the reference agonist are shown as green circles. The respective reference agonists from both studies are marked at zero. The error bars for the bias factors of all β2AR agonists represent S.E. with 95% confidence limits except for BI-167107 which is 99% confidence limit. The number of experiments for these cases were n ≥ 3. (B) The bias values obtained for the G-protein–biased system κOR is shown using U69593 as the reference agonist. For the opioid receptor agonists, the experimental bias factor was calculated from at least three independent experiments. All computational ligand bias values are reported with S.E. and were obtained from five independent MD simulations.
Comparison of Computational Ligand Bias to Experimental Bias Factor.
Figure 2 shows the comparison of the computed ligand bias to the experimental ligand bias factor for various β-arrestin and G-protein–biased agonists for the β2AR and κOR, respectively. In one set of experiments epinephrine was used as the reference ligand to calculate the bias factor from experimental titration curves (Rajagopal et al., 2011), and isoproterenol was used as a reference ligand in another set of experimental measurements (Weiss et al., 2013). Supplemental Fig. 4, A and B, shows the same comparisons for 5-HT receptors and β2ARTYY mutant.
As seen in Fig. 2 the computed ligand bias recapitulates the rank ordering of the ligand bias factors for all the agonists in all the four receptors for both β-arrestin and G-protein–biased agonists. It is important to note that the ligand bias ordering has been captured for the G-protein–biased agonists where we used a homology model for κOR. The magnitude of computed ligand bias shown in Fig. 2A is an order of magnitude lower than the experimental ligand bias factor. This is because we have not computed the macroscopic property of the ligand bias factor as measured by experiments. Instead we have computed the relative strength of the allosteric communication pipelines for the test ligand to that of the reference ligand and this is a molecular factor that potentiates and leads to change in ligand potency and efficacy that manifests as ligand bias. Thus, the computational method provides a rapid method to identify biased ligands ahead of experiments, and is useful in designing other biased ligands.
The agonists epinephrine, isoproterenol, and salbutamol are considered balanced agonists to β2AR (Rajagopal et al., 2011) whereas formoterol, 5-hydroxy-8-[1-hydroxy-2-[[2-methyl-1-(2-methylphenyl)propan-2-yl]amino]ethyl]-4H-1,4-benzoxazin-3-one (BI-167107), and carvedilol are β-arrestin–biased, with carvedilol acting as an inverse agonist to G-protein (Wisler et al., 2007). In the study by Rajagopal et al. (2011), the bias factors of epinephrine, isoproterenol, formoterol and salbutamol was measured, and the investigators concluded that formoterol was a β-arrestin–biased agonist, while the remaining ligands failed to exhibit any considerable bias toward either arrestin or G-protein.
As shown in Fig. 2A, in agreement with experiments, formoterol is more biased toward β-arrestin compared with isoproterenol or salbutamol in β2AR. Similarly, the super agonist BI-167107 is more biased toward β-arrestin compared with isoproterenol. Carvedilol is a weak β-arrestin–biased agonist and G-protein antagonist, so experimental bias factor values could not be obtained. However, the computational ligand bias values indicate that carvedilol is indeed β-arrestin–biased. As a blind test to our method using MD simulations and calculation of ligand bias, we further predicted that carmoterol is a β-arrestin–biased ligand compared with isoproterenol. We then performed the cAMP and the PathHunter assays as described in the Results section and calculated the bias factors shown in Supplemental Table 3. As predicted, carmoterol is indeed β-arrestin–biased with reference to isoproterenol.
In κOR (Fig. 2B), we correctly predicted that the agonists 1.1 and 2.1 are G-protein–biased compared with N-methyl-2-phenyl-N-[(5R,7S,8S)-7-pyrrolidin-1-yl-1-oxaspiro[4.5]decan-8-yl]acetamide (U69593). The experimental results (Zhou et al., 2013) showed that when using the ligand U69593 as the reference agonist, probe 1.1 and probe 2.1 are G-protein–biased, with probe 1.1 being a stronger G-protein–biased agonist compared with 2.1. Theoretical ligand bias calculations were obtained for these three ligand-receptor systems, and the values were found to be in excellent agreement with the experimental bias factor values, not only indicating the G-protein bias of both probes with respect to U69593, but also reflecting the greater G-protein bias for probe 1.1 compared with probe 2.1.
Comparison of the computational ligand bias for ergotamine in 5-HT1B and 5-HT2B as well as epinephrine in β2AR and the β2ARTYY mutant showed good agreement with the experimental bias factors reported (Supplemental Fig. 4). The serotonin receptors are an example of a system with different receptor serotypes in complex with the same ligand, and one receptor system being neutral while the other is β-arrestin–biased. The computational ligand bias values were in agreement with experimental bias factor values. There is no quantitative experimental bias factor available for the β2ARTYY mutant.
These results for several systems and different types of agonists show the robustness and the reliability of our method in predicting the ligand bias ahead of experiments. As described in detail in the Materials and Methods section, the confidence level for the prediction of ligand bias is 90%. Thus, the theoretical ligand bias calculations were validated against experimental bias factor values, and they were found to be clearly indicative of whether a preference for either G-protein and/or β-arrestin exists.
Identifying Functional Hotspot Residues in the Ligand Binding Site That Contribute to the Bias.
Figure 3 shows the comparison of the shape of the ligand binding sites in epinephrine (neutral agonist) and formoterol (β-arrestin–biased agonist), respectively, bound to the fully active state of β2AR as modeled in this work. For fair comparison, this figure shows the same view in the same receptor orientation for both agonists. Unlike epinephrine, the binding site of formoterol extends to the extracellular loop 2 (ECL2) region of the receptor. The residues D1133.32, F193ECL2, S2035.42, and S2045.43—although present in both the ligand binding sites—showed different side chain conformations (Supplemental Fig. 5) in formoterol-bound β2AR compared with the epinephrine-bound β2AR. Thus, changes in the rotameric states of these residues reshape the common regions of the cavities forming the ligand binding site for both epinephrine and formoterol.
The shape of the agonist binding cavity in β2AR of (A) a neutral agonist epinephrine and (B) formoterol, a biased agonist. The allosteric communication pipelines to the G-protein and the β-arrestin coupling interfaces are also depicted. This figure shows β2AR surfaces sliced to reveal the shape of the ligand binding cavity.
The allosteric communication pipelines from the EC loops passing through the residues in the formoterol-binding site to the β-arrestin coupling interface are stronger in formoterol-β2AR complex (as shown in Fig. 3B) than in the epinephrine-β2AR complex. However, it is interesting to note that the allosteric communication pipelines to the G-protein interface are equally strong in both epinephrine and formoterol (Fig. 3A). Both epinephrine and formoterol show strong G-protein communication pipelines in the fully active state of the receptor, and the allosteric pipeline strength to the β-arrestin coupling interface is comparatively weaker in epinephrine than in formoterol. This is also true in the experimental assays, where Rajagopal et al. (2011) showed that the Emax for formoterol and epinephrine was similar for the cAMP response while that for β-arrestin recruitment was lower for epinephrine than for formoterol. To date there is no known biased agonist that shows a strong preference for β-arrestin and not for G-protein in β2AR. It should be noted that we have used the fully active G-protein–bound conformation of β2AR for all these calculations, and the allosteric pipeline calculations are sensitive to the conformational state of the receptor (Bhattacharya and Vaidehi, 2014a).
One of the challenges in designing biased agonists is the lack of structural details of biased agonist binding sites that could aid structure–activity relationship (SAR) studies. The binding site of biased agonists looks similar to that of balanced agonists from crystal structures (Warne et al., 2012). Therefore, it would be beneficial to SAR studies if we could delineate the residues in the biased agonist binding site that contribute to the ligand bias. Toward this goal, we identified the residues in the biased agonist binding site that showed a significant contribution to the allosteric communication pipelines of either the β-arrestin or G-protein coupling interface by calculating the differences in the strengths of the allosteric communication pipelines to the BAI and GPI for the biased and unbiased agonists.
Figure 4 show the residues within 5 Å of the formoterol- and epinephrine-binding sites, respectively. The larger circles show sustained contact with the ligand. There are many residues that are common to both the epinephrine- and formoterol-binding sites. As seen in Fig. 4A, the four residues T1103.29, F193ECL2, S2045.43, and Y3087.35 show a significant contribution to the allosteric communication pipeline to the β-arrestin interface compared with the G-protein pipelines with reference to epinephrine (shown in Fig. 4B). These residues are the functional hotspot residues that confer β-arrestin bias to the formoterol-β2AR pairing. Carmoterol-bound β2AR also shows T1103.29, F193ECL2, and S2045.43 as functional hotspots in its binding site. Thus, three functional hotspot residues (T1103.29, F193ECL2, and S2045.43) are common to both formoterol and carmoterol. This shows that GPCRs might have specific residues that show functional selectivity to signaling pathways in the receptor.
(A) The residues within 5Å of formoterol in β2AR. (B) The binding site residues within 5 Å of epinephrine in β2AR. The size of the circles is proportional to the percentage of snapshots from the MD simulations that show the contacts. Large circles: agonist-residue contacts present in more than 75% of the snapshots from MD simulations. Medium-size circles: between 60% and 75% of the snapshots. Small circles: between 40% and 60% of the snapshots. Residues shown in red text with yellow highlighting are the functional hotspot residues.
The 3D view showing the relative orientation of these residues in the formoterol-binding site is shown in the Supplemental Fig. 6. To identify which chemical moieties in formoterol contribute to bias, we identified the moieties in formoterol that contact the functional hotspot residues in its binding site. The functional hotspots T1103.29 and Y3087.35 showed strong contact with the 4-methoxyphenyl-propanyl group attached to the protonated amine (Supplemental Fig. 6). This moiety that contacts the functional hotspots was absent in epinephrine. Therefore, extending the agonist beyond the protonated amine group and involving contacts with T1103.29 and Y3087.35 could yield biased agonists in β2AR.
Identifying Functional Hotspot Residues in the G-Protein–Biased Agonist Binding Site in κOR.
Figure 5 shows the residues within 5 Å of the G-protein–biased κOR agonist probe 1.1 and that of the unbiased agonist U69593, respectively. U69593 is predominantly nested in the region between transmembrane helices TM2, TM3, and TM7 (Fig. 5B). On the other hand, the G-protein–biased agonist 1.1 is more flexible in the binding site and made contacts with residues in TM2, TM3, TM5, and TM6, as shown in Fig. 5A. The residues Y652.64, D843.32, K1735.39, I2406.55, and Y2587.35 are functional hotspot residues in the 1.1 biased agonist binding site that contribute significantly to the strength of the allosteric communication pipeline to the G-protein interacting interface compared with the β-arrestin interacting interface. The residues Y652.64, D843.32, I2406.55, and Y2587.35 are common functional hotspot residues to both the 1.1 and 2.1. The residues D843.32 and Y652.64 are common to both the 1.1 and U69593 binding sites, but the residues K1735.39, I2406.55, and Y2587.35 are unique to 1.1. The residues D843.32 and Y652.64 also show a difference in the side-chain rotamer conformations when bound to U69593 or probe 1.1 (Supplemental Fig. 7). A Y258A7.35 mutation in the κOR has a strong deleterious effect on the potencies of agonists U69593, dynorphin A, the octahydroisoquinolinone carboxamide 1xx, and salvinorin A (Vardy et al., 2013). Mutation of residues D843.32 to Ala/Asn and Y652.64 to Ala/Phe negatively affected the affinity and potency of the receptor in the presence of U69593, 1xx, and salvinorin A (Vardy et al., 2013). In summary, we found that the G-protein–biased agonists in κOR displayed a rather rigid side-chain behavior, whereas, the β-arrestin–biased agonists such as formoterol in β2AR displayed a more flexible side-chain rotamer distribution. Supplemental Fig. 8 shows the 3D orientation of the probe 1.1 bound to κOR.
(A) The residues within 5Å of the G-protein–biased agonist 1.1 in κOR. (B) The binding site residues within 5 Å of the balanced agonist U69593 in κOR. The size of the circle is proportionate to the percentage of snapshots from the MD simulations that show the contacts. Large circles: agonist-residue contacts present in more than 75% of the snapshots from MD simulations. Medium-size circles: between 60% and 75% of the snapshots. Small circles: between 40% and 60% of the snapshots. Residues shown in red text with yellow highlight are the functional hotspot residues.
Thus, the functional hotspot residues in the ligand binding site identified by our method can be used 1) for generating a pharmacophore for identifying newer β-arrestin–biased ligands for β2AR and 2) for SAR studies for enhancing or improving the bias factor.
Discussion
Using MD simulation trajectories and the Allosteer method to calculate the allosteric communication pipelines, we have developed a computational method to calculate the ligand bias in GPCRs. The ligand bias is defined as the relative strength of allosteric communication pipelines starting from the residues in the EC loops passing through the residues in the agonist binding site to the residues in the β-arrestin or G-protein interacting interface in the GPCR. We found that the calculated ligand bias shows the same ordering as the experimental bias factors. Using the computations, we correctly predicted that carmoterol is β-arrestin–biased in β2AR with respect to isoproterenol. We have identified the functional hotspot residues in both β2AR and κOR, as shown in Figs. 4A and 5A, respectively, that can be used in designing β-arrestin– or G-protein–biased ligands.
Even to date, the identification and optimization of biased ligands is serendipity-driven. Usually the starting points for ligand optimization are identified in finding campaigns, such as high-throughput screening, fragment-based screening, or virtual screening. The hit compounds are tested for their ligand bias experimentally, and the promising starting points for further optimization have to show the desired biased profile to begin. The knowledge of functional hotspots can revolutionize this canonical hit-finding process. Using our computational method can be beneficial in two ways: 1) to identify biased agonists from virtual screening and 2) to use the functional hotpots to modify any hit ligand as a starting point for bias optimization. This allows us to pick starting points for lead optimization based on crucial parameters such as pharmacokinetic properties and off-target profiles rather than the bias of the initial hit. The ligand bias can be built in by modifications to contact the functional hotspots.
Another important, sought-after tool in the pharmaceutical industry is the ability to assess whether a particular target GPCR is amenable to bias or “biasable,” similar to assessing the “druggability” of a target protein. Therefore, we asked the question, which of the two receptors, β2AR or κOR, is more “biasable”? Figure 6 shows a heat map of the residue contribution to the differential β-arrestin and G-protein allosteric communication pipelines in the formoterol-bound active state of β2AR and the 1.1-bound active state of κOR, respectively. The heat map is an indication of the difference in the strength of allosteric communication by each residue toward β-arrestin in β2AR and toward G-protein in κOR. The more red the residue and thicker the size of the backbone cartoon in Fig. 6, the greater is its contribution to the ligand bias. We observed that there are large clusters of residues on TM5 and TM6 with high strength of allosteric communication toward the G-protein coupling interface in κOR. Such strong functional hotspot residues are absent in the β2AR. This could mean that κOR is probably more readily biasable compared with β2AR.
Heat map depicting the extent of bias in β2AR and in κOR. The thickness of the cartoon is directly proportional to the contribution of that residue to arrestin or G-protein bias. (A) β2AR in complex with formoterol (β-arrestin–biased). Regions of the receptor that are strongly β-arrestin–biased are represented by shades of blue to red. (B) κOR in complex with probe 1.1 (G-protein–biased). Regions of the receptor that are strongly G-protein–biased are represented by shades of blue to red.
In summary, our computational method to calculate the ligand bias and identify functional hotspot residues in GPCRs is extensible to other class A GPCRs, which will greatly aid in biased ligand design for GPCRs.
Authorship Contributions
Participated in research design: Nivedha, Tautermann, Vaidehi.
Conducted experiments: Nivedha, Bhattacharya, Lee, Kollak, Kiechle.
Performed data analysis: Nivedha, Tautermann, Bhattacharya, Vaidehi.
Wrote or contributed to the writing of the manuscript: Nivedha, Tautermann.
Footnotes
- Received August 24, 2017.
- Accepted January 16, 2018.
This work was funded by Boehringer-Ingelheim and by National Institutes of Health National Institute of General Medical Sciences [Grant R01-GM097296].
↵
This article has supplemental material available at molpharm.aspetjournals.org.
Abbreviations
- 3D
- three-dimensional
- β2AR
- β 2-adrenergic receptor
- BAI
- β-arrestin interface
- BI-167107
- 5-hydroxy-8-[1-hydroxy-2-[[2-methyl-1-(2-methylphenyl)propan-2-yl]amino]ethyl]-4H-1,4-benzoxazin-3-one
- EC
- extracellular
- ECL2
- extracellular loop 2
- GPCR
- G-protein-coupled receptor
- GPI
- G-protein interface
- 5-HT1B
- 5-hydroxytryptamine (serotonin) receptor 1B
- 5-HT2B
- 5-hydroxytryptamine (serotonin) receptor 2B
- MD
- molecular dynamics
- κOR
- κ-opioid receptor
- PDB
- Protein Data Bank
- PR
- allosteric pipeline ratio
- SAR
- structure–activity relationship
- TM
- transmembrane domain
- U69593
- N-methyl-2-phenyl-N-[(5R,7S,8S)-7-pyrrolidin-1-yl-1-oxaspiro[4.5]decan-8-yl]acetamide.
- Copyright © 2018 by The American Society for Pharmacology and Experimental Therapeutics