Design of Drug Efficacy Guided by Free Energy Simulations of the β 2 -Adrenoceptor

,


Introduction
The G-protein-coupled receptor (GPCR) superfamily constitutes the largest group of eukaryotic membrane proteins and plays essential roles in cellular communication.GPCRs share a common architecture composed of seven transmembrane (TM) helices and have extracellular (orthosteric) binding sites that recognize signaling molecules. [1]Development of GPCR ligands has had a tremendous impact on human health, and 34 % of the currently approved drugs interact with members of this receptor family. [2]Access to chemical probes has also been crucial for characterizing signaling pathways and revealed that GPCR activation is more complex than initially anticipated. [3]PCRs exist in an equilibrium between different conformations that can be modulated by ligand binding. [4,5] thosteric ligands are traditionally classified as agonists, antagonists, or inverse agonists.The binding of an agonist shifts the equilibrium towards active receptor conformations, which are able to interact with intracellular partners (e.g.G-proteins and β-arrestins) and thereby trigger signaling cascades.The receptor response is agonist-specific and compounds can produce the maximal stimulation (full agonist) or lower efficacy (partial agonists).Antagonists also bind to the orthosteric site but do not influence the conformational equilibrium between active and inactive states.Inverse agonists will stabilize inactive receptor conformations that are unable to stimulate intracellular partners.8] Accurate predictions of how ligand binding influences the signaling of a receptor would be valuable in drug design.
Activation of the β 2 adrenergic receptor (β 2 R) has been extensively characterized by biophysical and pharmacological methods. [6,7,9,10] Hi[13][14] In the intracellular region, activation involves a large (14 Å) outward displacement of TM6, leading to the formation of the G-protein binding site.In contrast, conformational changes in the orthosteric binding site are relatively small.The main difference between active and inactive conformations is a 2 Å inward bulge of TM5 (Figure 1A), leading to hydrogen bonding between adrenaline and Ser203 5.42 and Ser207 5.46 (superscripts refer to Ballesteros-Weinstein numbering). [3,14,15] Tere is a strong allosteric coupling between the orthosteric and G-protein binding sites.[18] These results have important implications for structure-based drug design.If structures of the active and inactive receptors are available, the efficacy of a ligand could be predicted by identifying compounds with a difference in affinity for these conformations.By considering the conformational selectivity in the hit-to-lead generation phase of drug discovery, the efficacy profile could be optimized to achieve maximal therapeutic effect with minimal side effects.
[21][22][23][24][25][26][27] Although it should in theory be possible to use brute force simulations to predict ligand efficacy, such calculations are prohibitively expensive to carry out because activation occurs on millisecond time scales.Enhanced sampling simulations have enabled detailed studies of the β 2 R activation mechanism, but these techniques are limited to studying small sets of ligands.We explored an alternative approach that does not require simulations of the full activation process.[30][31][32] We first assessed if free energy simulations could reproduce experimentally determined ligand affinities for the β 2 R stabilized in active or inactive conformations.Based on these results, we also evaluated if the calculated conformational selectivity of a ligand was correlated with its efficacy in functional assays.In a second step, we computationally designed 15 compounds and experimental evaluation of these led to the discovery of several agonist scaffolds with nanomolar potencies.Our results illustrate how leads with tailored signaling efficacy can be designed using structure-based methods.

Simulations of one receptor state predict ligand affinity, but not efficacy
MD simulations of agonists and antagonists bound to the β 2 R were performed to assess if free energy calculations could be used to predict ligand affinity and efficacy.The calculated binding free energies were compared to experimental data for 25 β 2 R ligands from three studies, which measured ligand affinities to active and inactive receptor conformations stabilized by nanobodies [18] or mutants. [33,34] ereas traditional measurements of binding affinities will reflect binding to both active and inactive receptors, this data enabled a comparison between experiments and simulations for a specific conformational state.The simulations were initiated from crystal structures of the active and inactive receptor embedded in a hydrated lipid bilayer.The overall receptor conformation was maintained in an active-or inactive-like conformation in the simulations by using spherical boundary conditions.The spherical system was centered on the binding site and was treated as fully flexible, whereas atoms outside a radius of 25 Å were tightly restrained.For each ligand and receptor state, binding affinities were calculated relative to a reference ligand (adrenaline or isoprenaline).Free energies were calculated by combining MD simulations with the FEP technique. [35]ach compound was alchemically transformed into the reference ligand in complex with the receptor and in aqueous solution, which can be used to calculate relative ligand affinities using the FEP method.Large agonists and antagonists (e.g., salmeterol and carvedilol) that extended into secondary binding pockets were not considered because transformations involving a large number of atoms lead to difficulties to obtain converged free energy estimates.The free energies were calculated based on a perturbation network, [36] which was optimized to achieve converged simulation results (Supporting Information Figure S1-S3).Several transformation pathways were used to calculate the free energy for each compound pair and convergence was assessed using the cycle closure error (average error for the inactive and active conformation = 0.7 and 0.6 kcal mol À 1 , respectively).Predicted and experimental relative binding free energies for the active (R*) and inactive (R) states (ΔΔG bind,R* and ΔΔG bind,R , respectively) were in good agreement with a mean unsigned error (MUE) of 1.5 kcal mol À 1 and a coefficient of determination (R 2 ) of 0.72 (Figure 1B, Supporting Information Tables S1 and S2).It should be noted that some antagonists have higher affinities for the active receptor conformation than potent agonists, illustrating that calculations based on a single receptor state cannot predict ligand efficacy.

Calculated affinity shifts upon receptor activation predict ligand efficacy
The ability of a compound to trigger receptor activation has been shown to correlate with the difference in ligand affinity for the active and inactive conformations. [18]Based on this observation, we assessed if ligand efficacy could be predicted by comparing the calculated relative binding free energies for the two receptor states (the relative affinity shift, ΔG shift = ΔG bind,R* À ΔG bind,R ) (Figure 2A).The predicted efficacy is determined by the calculated relative affinity shift and the choice of reference ligand.The agonists adrenaline and isoprenaline were selected as references.These agonists were selected as references because the compounds fit well into both the active and inactive binding sites.Larger compounds (e.g., the antagonist carazolol) were more challenging to model into the active state and longer simulations would be required to obtain converged free energy calculations.In this case, the relative affinity shift of an agonist is expected to be close to zero as these compounds have similar efficacy profiles.In contrast, ΔG shift would be expected to be positive for inverse agonists and antagonists because these compounds have an altered balance between the affinities to the two conformational states.Finally, partial agonists should display intermediate relative affinity shifts (Figure 2B).To assess our computational approach, we performed retrospective binding free energy simulations for a set of 31 β 2 R ligands, which were divided into three groups based on assays measuring G-protein mediated signaling (agonists, partial agonists, and non-agonists: weak partial agonists, antagonists, and inverse agonists) (Supporting Information Table S3).Predicted relative affinity shifts are shown in Figure 2C (Supporting Information Table S4) and were converged based on analysis of cycle closure errors (average error = 0.7 kcal mol À 1 , Supporting Information Figure S3).In agreement with experimentally observed efficacy profiles, the agonists typically showed the smallest relative affinity shifts (median ΔG shift of 1.3 kcal mol À 1 ) whereas the non-agonists had the largest shifts (median ΔG shift of 6.6 kcal mol À 1 ).Partial agonists showed intermediate values with a median ΔG shift of 3.3 kcal mol À 1 .The statistical significance of these predictions was confirmed by using a Kruskal-Wallis test (pvalue < 10 À 4 ) followed by a Dunn's post hoc test.The average affinity shifts were significantly different between agonists and the two other groups (partial agonists and nonagonists with p-values < 0.05 and < 0.001, respectively).For comparison, ranking compounds by their calculated relative affinity to the active state did not show any significant difference between the three groups (Kruskal-Wallis test: pvalue > 0.05, Supporting Information Figure S4).

Computational design of non-catechol β 2 R agonists
A major challenge for drug discovery targeting adrenergic receptors has been to replace the catechol group and maintain agonist activity.To identify novel non-catechol agonists, two virtual chemical libraries were generated by replacing the catechol group in the β 2 R agonist colterol with other building blocks (Figure 3A).Building blocks containing either five-or six-membered aromatic rings with at least one hydrogen bond donor or acceptor were identified, resulting in two libraries with a total of � 34 000 unique compounds.Each library was first docked to the active β 2 R conformation and top-ranked molecules with a binding mode similar to adrenaline were visually inspected.Fifty-six compounds were selected based on interactions with key residues in the binding site (Asp113 3.32 , Asn293 6.55 , Ser203 5.42 , and Ser207 5.46 ), followed by free energy simulations in the active and inactive conformations (Supporting Information Figure S5).Based on the calculated affinity shifts and synthetic feasibility, seven compounds were selected from the two libraries (1-7, Scheme 1 and Table 1).Five compounds were based on replacing the catechol with sixmembered aromatic rings and, among these, one was predicted to be an antagonist (3, ΔG shift = 7.4 kcal mol À 1 ), one to be a weak partial agonist (2, ΔG shift = 5.0 kcal mol À 1 ), and the remaining corresponding to full agonism (compounds 1 and 4-5, ΔG shift = 1.2, 1.1, and À 1.9 kcal mol À 1 , respectively).The two selected compounds with five-membered aromatic rings (6 and 7) were predicted to be full agonists (ΔG shift = À 0.8 and 1.6 kcal mol À 1 , respectively, Supporting Information Table S5).

Biological activity of designed ligands
Compounds 1-7 were first evaluated in radioligand binding assays (Table 1) and five of these (1, 2, 5-7) were ligands of the β 2 R. Compounds 1 and 2 had K i values of 17 and 58 nM, respectively, and these ligands were predicted to be full and weak partial agonists, respectively (Table S5).The compounds were based on indazole (1) and indole (2) moieties that mimic the catechol of adrenaline (Figure 3B).We do not consider compounds 1 and 2 to represent new scaffolds because indoles have previously been explored as high affinity β 2 R ligands (e.g., in carazolol and pindolol). [39,40] wever, it should be noted that an indole-like moiety does not necessarily result in a specific efficacy because this scaffold is part of both agonists and antagonists.Measurements of G s -mediated receptor activation in response to the compounds revealed that both compounds were partial agonists with maximal efficacies of 61 % (EC 50 = 8.5 nM) and 48 % (EC 50 = 55 nM) for 1 and 2 relative to noradrenaline, respectively (Table 1 and Figure 3C-E).Both compounds 1 and 2 showed low levels of β-arrestin recruitment with a maximal efficacy relative to noradrenaline of 19 % (1, EC 50 = 16.0 nM) and 7 % (2, EC 50 = 120 nM).The enantiopure form of compound 1 displayed very high affinity (compound 1 a, K i = 9.3 nM).Subsequent to our discovery of compound 1, this structure was described in a patent, which confirmed that the most active enantiomer is the (R)-form, in agreement with our predictions. [41]In contrast to the indole scaffold, the indazole is able to mimic two key hydrogen bonds of the catechol group (Ser203 5.42 and Asn293 6.55 ), which explains the higher efficacy and potency of this agonist (Figure 3B).

Fine-tuning signaling by perturbing receptor-ligand interactions
Five analogs of 1 and 2 were synthesized to evaluate the ability of free energy simulations to predict the effect of small structural modifications on ligand efficacy.Compounds were designed to form an additional hydrogen bond with Ser207 5.46  (11) or to disrupt interactions with Asn293 6.55 and Ser203 5.42 (12-15) (Figure 4A).Detailed synthetic procedures for 11-15 are described in the Supporting Information.
Predictions from the free energy simulations were compared with results from binding and functional assays.The addition of an exocyclic nitrogen to the indazole moiety to form a hydrogen bond with Ser207 5.46  (11) decreased the calculated affinity for both the active and inactive states compared to the reference agonists isoprenaline and compound 1, which agreed with the experimentally determined loss of affinity (K i = 490 nM) (Supporting Information Table S5).The predicted ΔG shift value of 0.6 kcal mol À 1 was compatible with agonist activity, which was confirmed by the functional assays (E max = 66 % in the G-protein signaling assay).The addition of a 2-methyl group to the indole moiety to probe interactions with Asn293 6.55  (12) resulted in a ΔG shift of 3.9 kcal mol À 1 , corresponding to weak partial agonist activity.Experimentally, compound 12 resulted in a K i of 4.0 nM and partial agonist activity at the G s proteinmediated signaling (E max = 38 %), in agreement with the prediction that this compound would show lower efficacy than compound 1.Finally, a benzothiophene (13) and benzofuran (14 and 15) group were introduced to probe interactions with Ser203 5.42 and Asn293 6.55 , which resulted in affinity shifts between 6.4 and 6.7 kcal mol À 1 , corresponding to values observed for antagonists.This result also agreed well with the low levels of receptor activation observed experimentally (E max = 16, 19, and 18 % for compounds 13, 14, and 15, respectively, for G-protein signaling).To assess the ability of the free energy simulations to predict differences in ligand efficacy profiles in congeneric series of ligands, the affinity shifts of these compounds, as well as four reference ligands, were compared to the response in assays measuring G s -mediated signaling (Figure 4B).There was a strong correlation between the predicted affinity shifts and receptor activation rate (Spearman correlation coefficient = À 0.85, p-value < 0.05).

Discovered agonists stabilize active-like receptor conformations
To further characterize the identified agonists, enhanced sampling MD simulations of complexes with compounds 1, 7, and 11 were carried out.In these simulations, the receptor was fully flexible and the enhanced sampling protocol ensured that the active-like free energy basin was exhaustively sampled.These simulations thereby allowed us to study ligand-induced conformational changes.The results were compared with previously performed simulations of reference ligands (three agonists, two antagonists, and one inverse agonist) and in the absence of ligand (Supporting Information Figure S6). [26]In agreement with our experimental data, the conformational states stabilized by compounds 1, 7, and 11 were more similar to agonists than to non-agonists (Supporting Information Figure S6A).In the connector region, which is located close to the orthosteric site, compounds 1, 7, and 11 stabilized an active-like conformation of Ile121 3.40 and Phe282 6.44 , which was also observed in simulations of the reference agonists.(Supporting Information Figure S6B).

Strategies to design drugs with tailored efficacy profiles
Three observations that are relevant for the design of GPCR agonists emerge from this study.First, a strength of using a physics-based approach is that no prior knowledge regarding the pharmacophore features that are important for agonism is required.In theory, the free energy calculations only rely on access to structures of the active and inactive receptor conformations and a reference agonist.
It should be noted that the β 2 R represents a favorable case because the binding modes could be modelled with high confidence and key interactions for agonism have been identified.The same approach can be applied to the 25 GPCRs for which there are experimental structures of inactive and active conformations (Supporting Information Table S6).
A second interesting observation from our virtual screens was that the top-ranked compounds rarely were able to form the complete set of hydrogen bonds formed by adrenaline, which may be important to achieve full activation of β adrenergic receptors [42] and suggests that larger libraries may improve the hit rates.For example, molecular docking screens have recently identified GPCR agonists in libraries with > 100 million compounds. [43]Due to the computational cost of free energy simulations, our approach is limited to predictions of relative binding free energies for small sets of compounds and is not suitable for screens of diverse libraries.In order to ensure facile synthesis, only compounds that could be made from commercially available building blocks were part of our library with 34 000 compounds.If we instead had generated the library based on all theoretically possible building blocks up to 13 heavy atoms, [44] > 1 million additional compounds would be accessible (a 40-fold larger library).As this library expansion will come at the expense of synthesis efforts, it will become important to select compounds using reliable computational models.Previous structure-based virtual screens have identified β 2 R agonists, [45][46][47] but a large fraction of the compounds predicted by the approximate docking scoring functions have been inactive.By combining docking screens with rigorous free energy simulations of top-ranked compounds, the balance between speed and accuracy can be optimized.
A third important observation is that the simulations were able to predict structural modifications resulting in a spectrum of efficacy profiles for a congeneric series of compounds, which could be used to tune ligand efficiency in lead optimization.Notably, potent agonists designed from our simulations showed stronger activation of G-proteinmediated signaling compared to β-arrestin recruitment, suggesting that the simulated active conformation may represent a biased signaling state.Based on our experimental data, the functional selectivity of compound 1 a is comparable with the drug salbutamol, which has been described as a biased agonist. [6]Interestingly, the compounds have slightly different profiles.Whereas salbutamol showed higher ligand efficacy, compound 1 a is more potent with EC 50 values in the single-digit nanomolar range.The structural basis of the observed functional selectivity is unclear and is challenging to characterize without performing computationally demanding MD simulations of Gprotein and β-arrestin bound states.A limitation of our protocol is that the receptor is maintained in the active or inactive conformations that were captured by crystallography, whereas multiple distinct states are accessible in the experiments.Additional conformations could be identified by MD simulations of functionally selective ligands combined with the characterization of allosteric communication pipelines using the methods developed by Vaidehi and coworkers. [48,49] y extending the free energy simulations to include structures representing several different active conformations, it may be possible to rationally design functionally selective ligands.

Conclusion
Identification of lead compounds with a specific efficacy profile is one of the major challenges in GPCR drug discovery.Here, we show that free energy simulations can be used to predict the functional activity and demonstrated the efficiency of this approach by designing β 2 R agonists.The large number of available GPCR structures will make it possible to apply the same approach to other targets, which can accelerate the development of drug candidates with tailored signaling efficacy.

Figure 1 .
Figure 1.Molecular basis of β 2 R activation and relative binding free energy calculations in active and inactive receptor conformations.A) Structures of the active and inactive conformations of the β 2 R. Active (PDB accession code: 4LDO) and inactive (PDB accession code: 2RH1) structures are shown as blue and red cartoons, respectively.Adrenaline and key binding site residues are shown as sticks.Hydrogen bonds are depicted as yellow dashed lines.B) Correlation between calculated and experimental relative binding free energies of β 2 R ligands.Relative binding free energies of ligands (kcal mol À 1 ) were calculated relative to adrenaline (circle) or isoprenaline (triangle) in the active (blue) and inactive (red) conformations of the receptor.The free energies were calculated from three independent simulations and error bars correspond to the standard deviation.The mean unsigned error (MUE, kcal mol À 1 ) and the coefficient of determination (R 2 ) are shown in the top left corner.

Figure 2 .
Figure 2. Relative affinity shift calculations for β 2 R ligands.A) Thermodynamic cycle used to compute the relative affinity shift from free energy simulations.Alchemical transformations of the ligands A and B were performed in the active and inactive states of the receptor.B) Illustration of the expected affinity shifts for ligands with different efficacy profiles with a full agonist as reference compound.C) Calculated relative affinity shifts and experimental efficacy profiles of β 2 R ligands.Affinity shifts were calculated as the difference between the relative binding free energy of each ligand in the active and inactive state of the receptor with an agonist as the reference.Ligands are colored according to their efficacy profiles (Gprotein-mediated signaling): agonists (green), partial agonists (orange), and non-agonists (red).Error bars correspond to the standard deviation.The boxplot in the top left corner shows the affinity shift distribution of the different ligand groups.Boxes represent the 50th percentile of the data and the black bands show the median value.The whiskers correspond to 1.5 times the interquartile range.n.s.: not significant, * p � 0.05, ** p � 0.01, and *** p � 0.001 according to a Kruskal-Wallis test followed by a Dunn's test.

Figure 3 .
Figure 3. Computational design of non-catechol β 2 R ligands.A) A virtual chemical library was constructed by replacing the catechol group in the β 2 R agonist colterol with either five-and six-membered aromatic rings.B) Predicted binding modes of the designed compounds.Structures correspond to snapshots from an MD simulation of the β 2 R active conformation.The receptor is shown as a gray cartoon and key residues are shown as sticks.Ligands are represented as sticks and hydrogen bonds are depicted as yellow dashed lines.C) Radioligand displacement curves and binding affinities (K i ) of selected ligands at β 2 R derived from four (5, 6), five (2, 7), or six (1) experiments.D), E) Functional activity for Gprotein signaling (D) and β-arrestin recruitment (E).Graphs show mean curves � s.e.m. derived from three (5, 6 in E), four (1, 2, 5, 6 in D, 2 in E), or five (7 in D, 1, 7 in E) individual experiments.
[a] IP-One ® assay (Cisbio) with HEK293T cells co-transfected with β 2 R and the hybrid G-protein Gα qs .[b] PathHunter ® assay (Eurofins/DiscoverX) for measuring arrestin recruitment in EA-tagged arrestin-2 expressing HEK cells co-transfected with the ProLink tagged β 2 R-PK1.[c] K i values in [nM � s.e.m.] are the means of three to eight individual experiments each performed in triplicate with the number of experiments listed in parentheses (n).[d] EC 50 values for G-protein signaling in [nM � s.e.m.] indicating mean potencies derived from three to eight individual experiments each performed in duplicate.[e] E max values in [% � s.e.m.] indicating maximum efficacy relative to the full effect of noradrenaline (= 100 %) derived from (n) experiments.[f ] EC 50 values for β-arrestin recruitment in [nM � s.e.m.] derived from three to twelve individual experiments each performed in duplicate.n.d.= not determined.* Separated enantiomers.

Figure 4 .
Figure 4. Fine-tuning receptor signaling.A) Predicted binding modes of compounds 11-15.Structures correspond to snapshots after an MD simulation of the β 2 R active conformation.The β 2 R is shown as a gray cartoon, and key residues are shown as sticks.All ligands are represented as sticks with orange carbon atoms.Hydrogen bonds are depicted as yellow dashed lines.B) Correlation between predicted affinity shifts and G smediated signaling relative to noradrenaline.Error bars correspond to the standard deviation.The Spearman correlation coefficient (ρ) is shown in the top right corner.Abbreviations: Salbutamol (SAL), Noradrenaline (NORA), Colterol (COL), and Halostachine (HAL).

Table 1 :
Ligand binding and functional assays results for the tested compounds.