Ligand Binding Ensembles Determine Graded Agonist Efficacies at a G Protein-coupled Receptor*

G protein-coupled receptors constitute the largest family of membrane receptors and modulate almost every physiological process in humans. Binding of agonists to G protein-coupled receptors induces a shift from inactive to active receptor conformations. Biophysical studies of the dynamic equilibrium of receptors suggest that a portion of receptors can remain in inactive states even in the presence of saturating concentrations of agonist and G protein mimetic. However, the molecular details of agonist-bound inactive receptors are poorly understood. Here we use the model of bitopic orthosteric/allosteric (i.e. dualsteric) agonists for muscarinic M2 receptors to demonstrate the existence and function of such inactive agonist·receptor complexes on a molecular level. Using all-atom molecular dynamics simulations, dynophores (i.e. a combination of static three-dimensional pharmacophores and molecular dynamics-based conformational sampling), ligand design, and receptor mutagenesis, we show that inactive agonist·receptor complexes can result from agonist binding to the allosteric vestibule alone, whereas the dualsteric binding mode produces active receptors. Each agonist forms a distinct ligand binding ensemble, and different agonist efficacies depend on the fraction of purely allosteric (i.e. inactive) versus dualsteric (i.e. active) binding modes. We propose that this concept may explain why agonist·receptor complexes can be inactive and that adopting multiple binding modes may be generalized also to small agonists where binding modes will be only subtly different and confined to only one binding site.

G protein-coupled receptors constitute the largest family of membrane receptors and modulate almost every physiological process in humans. Binding of agonists to G protein-coupled receptors induces a shift from inactive to active receptor conformations. Biophysical studies of the dynamic equilibrium of receptors suggest that a portion of receptors can remain in inactive states even in the presence of saturating concentrations of agonist and G protein mimetic. However, the molecular details of agonist-bound inactive receptors are poorly understood.
Here we use the model of bitopic orthosteric/allosteric (i.e. dualsteric) agonists for muscarinic M 2 receptors to demonstrate the existence and function of such inactive agonist⅐receptor complexes on a molecular level. Using all-atom molecular dynamics simulations, dynophores (i.e. a combination of static three-dimensional pharmacophores and molecular dynamics-based conformational sampling), ligand design, and receptor mutagenesis, we show that inactive agonist⅐receptor complexes can result from agonist binding to the allosteric vestibule alone, whereas the dualsteric binding mode produces active receptors. Each agonist forms a distinct ligand binding ensemble, and different agonist efficacies depend on the fraction of purely allosteric (i.e. inactive) versus dualsteric (i.e. active) binding modes. We propose that this concept may explain why agonist⅐receptor complexes can be inactive and that adopting multiple binding modes may be generalized also to small agonists where binding modes will be only subtly different and confined to only one binding site.
Specific and coordinated cell-to-cell communication regulates the flow of information between cells, and proper information processing ensures physiological functions of biological systems. G protein-coupled receptors (GPCRs), 8 constituting the largest class of membrane proteins in mammals, are essential mediators of chemically and light-encoded information (1)(2)(3)(4). GPCRs sense a great variety of extracellular stimuli, e.g. neurotransmitters and hormones, and subsequently translate this information into an intracellular response via G proteins, ␤-arrestins, and possibly GPCR-interacting proteins (2)(3)(4)(5). Because of their abundance and relevance in regulating the majority of (patho-)physiological processes in humans, GPCRs have for a long time represented the most important drug targets being addressed by at least a third of all currently marketed drugs (6,7).
Agonist binding leads to receptor activation, which is followed by intracellular G protein recruitment and subsequent cell signaling. Breakthroughs in GPCR crystallography have led to inactive and active crystal structures of the same receptor protein. Among these are rhodopsin (8 -10) and more recently the ␤ 2 -adrenergic (11)(12)(13)(14), M 2 muscarinic (15,16), and -opioid receptors (17,18). These structures most likely represent energetically favored, relatively stable inactive and active receptor⅐ligand complexes. Despite the diversity of crystallized receptors, a common mechanism of receptor activation can be inferred from these structures in conjunction with a wealth of older biochemical data. Agonist-mediated receptor activation leads to an outward tilt of the intracellular parts of transmembrane helices V and VI, facilitating intracellular G protein binding, which in many cases is accompanied by an inward movement of the extracellular parts of transmembrane helices V and VI. This results in a contraction of the ligand binding site, closing it off toward the extracellular space.
Despite the wealth of GPCR structures, however, little is known about the structural dynamics of their inactive-to-active transitions and different agonist efficacies. Recent biophysical studies on various receptors (Class A and Class C) have shown that GPCRs reside in a dynamic equilibrium of distinct receptor conformations comprising inactive and active receptor states (19 -22). These studies suggest a common mechanism for agonist efficacy: agonists shift the preexisting equilibrium of different receptor conformations toward more active states (21,22).
Interestingly, the dynamic equilibrium of receptors remains heterogeneous even in the presence of saturating concentrations of full agonists and always contains a fraction of receptors in inactive states (19,(21)(22)(23). In fact, agonists alone are not sufficient to stabilize the fully active receptor state as seen in the crystal structures (19 -24). The fully active state is only reached upon addition of both saturating concentrations of agonist and a G protein or nanobody (19,21,23). However, even under these conditions, e.g. in the presence of the full agonist isoproterenol and the nanobody Nb80, a significant percentage of ␤ 2 -adrenergic receptors still remains in inactive states (21).
In line with this concept, partial agonists would stabilize fewer receptors in active states, and concomitantly a greater fraction of partial agonist-bound receptors would remain in inactive states (22,24). However, it is puzzling how agonistbound receptors can adopt active and inactive states. Based on the above mentioned evidence that receptors are "floppy," it is reasonable to hypothesize that agonists may have multiple binding modes; some of these agonist binding modes might then stabilize active agonist⅐receptor complexes, whereas other agonist binding modes might stabilize inactive agonist⅐receptor complexes.
However, because purely orthosteric partial agonists are likely to show only slightly different binding modes, it is technically challenging to investigate this phenomenon with current methods. In this study, we make use of a special case of agonism to identify multiple binding topographies of agonists and to provide a proof of principle that ligand binding ensembles can form the molecular basis of different agonist efficacies. This special case is dualsteric partial agonists for the muscarinic M 2 receptor (M 2 AChR) that comprise two pharmacophores targeting the orthosteric and allosteric binding sites, respectively (25). These dualsteric agonists have been suggested to possess two pharmacologically distinct binding modes (26 -28), a feature defined as dynamic ligand binding (28). The "extreme" molecular nature of the dualsteric agonists, i.e. spanning two binding sites, is anticipated to allow the identification of distinct ligand binding topographies.
By combining pharmacological methods, molecular docking, and all-atom molecular dynamics simulations based on M 2 AChR crystal structures, we here identify two distinct binding topographies of a set of dualsteric partial agonists. One binding mode resembles that of the co-crystallized orthosteric agonist iperoxo (16) and forms an active agonist⅐receptor complex. The other binding mode is purely allosteric and forms an inactive agonist⅐receptor complex by adopting binding topographies similar to those of common allosteric modulators (29). We show that both types of agonist⅐receptor complexes are present in a receptor population at the same time. Manipulation of the dynamic equilibrium of active and inactive agonist⅐receptor complexes by ligand design and by mutations reveals that the fraction of inactive agonist⅐receptor complexes decreases overall agonist efficacy.
Our findings suggest that an ensemble of active and inactive agonist⅐receptor complexes in which agonists adopt multiple binding topographies is a molecular basis for different agonist efficacies. Albeit having studied a highly specific case of GPCR agonism, our findings provide general evidence for the concept that agonists can principally stabilize inactive agonist⅐receptor complexes; this may also apply to agonists that target only one binding site of the receptor.

Results
Direct Pharmacological Evidence for Multiple Binding Modes of a Dualsteric Agonist-Dualsteric (also termed bitopic orthosteric/allosteric) ligands (30) for muscarinic receptors are composed of orthosteric and allosteric moieties covalently connected via linkers of different chemical nature (31)(32)(33)(34). Recent studies suggested that at least some of these bipharmacophoric ligands may have more than one binding mode (28,34). The experimental results are compatible with a theoretical scenario, i.e. dynamic ligand binding, in which such a ligand would adopt a dualsteric pose (Fig. 1a) and a purely allosteric binding mode (Fig. 1a). In this study, we provide direct evidence that iper-6naph indeed displays two distinct binding modes at M 2 AChRs. Iper-6-naph is a prototypical dualsteric ligand (Fig. 1a) that is composed of the orthosteric agonist iperoxo linked via a hexamethylene chain to an allosteric modulator derived from naphmethonium. Equilibrium binding of iper-6-naph using N-[ 3 H]methylscopolamine ([ 3 H]NMS) as the orthosteric radiotracer shows partial displacement of [ 3 H]NMS in a HEPES buffer (10 mM HEPES, 10 mM MgCl 2 , 100 mM NaCl) (Fig. 1b).
[ 3 H]NMS displacement indicates ligand competition between iper-6-naph and the radioligand at the orthosteric site. However, the incomplete [ 3 H]NMS displacement at high concentrations of iper-6-naph is indicative of the formation of ternary complexes consisting of [ 3 H]NMS and iper-6-naph bound to the receptor's orthosteric and allosteric sites, respectively. The binding of iper-6-naph to the allosteric site hampers [ 3 H]NMS dissociation, which is detected as incomplete displacement (Fig. 1b). To directly prove that iper-6-naph can indeed adopt a purely allosteric binding pose, we conducted [ 3 H]NMS displacement experiments in a 5 mM Na ϩ , K ϩ , P i buffer (for exact composition see "Experimental Procedures"). This buffer has been shown to increase the affinity of various allosteric modulators, most likely due to the absence of Mg 2ϩ cations, which have been suggested to compete with allosteric modulators for the allosteric binding site at muscarinic M 2 receptors (35). Under these conditions, iper-6-naph leads to a concentrationdependent increase of [ 3 H]NMS binding, which is direct evidence for an allosteric binding topography of iper-6-naph (Fig.  1b). On the molecular level, enhancement of [ 3 H]NMS binding is due to positive cooperativity of the allosteric moiety 6-naph with [ 3 H]NMS (Fig. 1c). The 5 mM Na ϩ , K ϩ , P i buffer enhances both allosteric affinity and efficacy of 6-naph (Fig. 1c). However, the 5 mM Na ϩ , K ϩ , P i buffer has little to no effect on the binding of the orthosteric agonist iperoxo (Fig. 1d). These data demonstrate that iper-6-naph adopts two distinct binding modes at M 2 AChRs that reside in dynamic equilibrium.
Structural Characteristics of a Ligand Binding Ensemble-The concept that ligands can adopt multiple distinct binding modes at the same GPCR has been defined as dynamic ligand binding (28) (Fig. 1a). However, the structural basis for this phenomenon remains unclear. To demonstrate that a GPCR is indeed capable of accommodating a ligand in multiple binding modes, i.e. forming a ligand binding ensemble, we here apply an interdisciplinary approach combining molecular modeling techniques guided by pharmacological experiments and vice versa. Our modeling strategy consists of a combination of allatom molecular dynamics (MD) simulations and three-dimensional pharmacophores (36,37).
The M 2 AChR displays strong allosteric coupling between the ligand binding site and the intracellular G protein-binding interface (16). Structurally, this is represented by a pronounced contraction of the orthosteric binding site (16). In fact, iperoxo binding leads to a complete closure of the orthosteric binding site toward the extracellular space (16). The so-called tyrosine lid formed by Tyr-104 3.33 , Tyr-403 6.51 , and Tyr-426 7.39 separates the orthosteric binding site from the allosteric binding site in the active M 2 AChR crystal structures (Fig. 2a). However, the side chain conformations of these tyrosine residues appear very flexible in MD simulations, and the receptor tolerates agonist binding within an orthosteric binding site that is not completely closed. Based on the available crystal structures of the M 2 AChR in active (Protein Data Bank codes 4MQS and 4MQT) (16) and inactive conformations (Protein Data Bank code 3UON) (15), we first generated receptor⅐ligand complexes by docking experiments. To dock iper-6-naph into the active M 2 AChR structure (Protein Data Bank code 4MQT), the tyrosine lid was remodeled by side chain sampling based on MD simulations of the iperoxo-bound crystal structure (Fig. 2, b and c). MD-based side chain sampling detects an open tyrosine lid suitable for docking bitopic agonists. The only differences between the crystal structure and the obtained active-like receptor model are the side chain conformations of Tyr-104 3.33 , Tyr-403 6.51 , and Tyr-426 7.39 .
Docking and subsequent all-atom MD simulations initiated from an active receptor⅐ligand complex reveal a dualsteric binding topography of iper-6-naph (Fig. 3a). The iperoxo moiety of iper-6-naph binds to the orthosteric binding site, whereas the 6-naph moiety protrudes toward extracellular domains, thereby engaging residues of the common allosteric binding site (Fig. 3a). The dualsteric binding pose of iper-6-naph stabilizes an active ligand⅐receptor complex because the agonistic properties of iper-6-naph are mediated by the orthosteric site (25) and hence are only compatible with a binding pose in which the iperoxo moiety targets the orthosteric binding site. Based on the MD simulations, additional molecular descriptors can be inferred that argue for an active ligand⅐receptor complex. First, the orientation of the iperoxo moiety of iper-6-naph is almost identical to that observed in the iperoxo-bound M 2 AChR crystal structure (Fig. 3, a and b). In both iperoxo and iper-6-naph MD simulations of the active M 2 AChR complexes, the positively charged nitrogen of iperoxo interacts with Asp-103 3.32 and displays -cation interactions with Tyr-104 3.33 , Tyr-403 6.51 , and Tyr-426 7.39 (Fig. 3, a and b). The triple bond exhibits hydrophobic contacts with Tyr-104 3.33 , Trp-155 4.57 , and Trp-400 6.48 . Additionally, Asn-404 6.52 forms a hydrogen bond with the oxygen of the 4,5-dihydroisoxazole moiety of iperoxo (Fig. 3a). Second, the volume of the orthosteric binding site of iper-6-naph⅐receptor complexes (249.6 Å 3 ) is similar to the orthosteric volume of the active iperoxo⅐receptor structure (222.3 Å 3 ) and that obtained by MD simulation of active iperoxo⅐receptor complexes (242.7 Å 3 ). In line with this, the distance between the C␣ atoms of Tyr-104 3.33 and Asn-404 6.52 as an indicator for the extent of contraction of the orthosteric binding site upon activation is similar in the iper-6-naph (13.1 Å) and iperoxo MD simulations (12.4 Å). In contrast, the crystal structure of the inactive M 2 AChR shows a larger volume of the orthosteric binding site (383.6 Å 3 ) and a larger distance between Tyr-104 3.33 and Asn-404 6.52 (13.6 Å).
Docking of iper-6-naph into the inactive M 2 AChR crystal structure (Protein Data Bank code 3UON) and subsequent MD simulations reveal a second binding mode of iper-6-naph that is purely allosteric (Fig. 3c). Overall, the purely allosteric binding pose of iper-6-naph is similar to that of the allosteric modulator naphmethonium (Fig. 3d). Iper-6-naph in its purely allosteric binding mode interacts with two centers formed by aromatic residues in the allosteric vestibule (Fig. 3c). The first center is lined by Tyr-177 EL2 and Trp-422 7.35 , and the second is lined by Tyr-80 2.60 and Tyr-83 2.63 as reported previously for prototypical allosteric modulators (29). In line with this, the ammonium group of the iperoxo moiety forms cationinteractions with Tyr-80 2.60 and Tyr-83 2.63 , and the other ammonium group displays cation-interactions with Trp-422 7.35 (Fig. 3, c and d). In addition, the triple bond is located opposite to Tyr-177 EL2 , and the two methyl groups of the small allosteric linker chain show lipophilic contacts to Ala-191 5.43 , Tyr-403 6.51 , and Tyr-426 7.39 (Fig. 3c). The naphthalene ring is embedded in a pocket built by Thr-187 5.39 , Phe-188 5.40 , and Val-407 6.55 . Moreover, the ammonium group of iperoxo forms a charge interaction with Glu-175 EL2 , and the other ammonium group features cationinteractions with Trp-99 3.28 and Tyr-426 7.39 (Fig. 3, c and d).
Functionally, the allosteric binding pose of iper-6-naph stabilizes an inactive ligand⅐receptor complex as its topography is comparable with the allosteric ligand naphmethonium (Fig.  3d). Naphmethonium and, in particular, its building block 6-naph are allosteric inverse agonists (28). Our MD simulations indicate that iper-6-naph is able to adopt two completely distinct binding modes at the same receptor: the dualsteric iper-6-naph⅐M 2 AChR complex resembles an active receptor, whereas the purely allosteric iper-6-naph⅐M 2 AChR complex is inactive.
We use pharmacological experiments to quantify the fractions of receptors that are bound to iper-6-naph in either the dualsteric or the allosteric binding pose (Fig. 3e)   over, the global analysis yields the dynamic efficacy, dyn , of iper-6-naph (for fitting details see "Experimental Procedures"). This new global analysis is superior to previous analyses (25,26,28,34) as it yields the allosteric affinity K B of a dualsteric ligand that has not been experimentally accessible before. Iper-6-naph displays partial M 2 AChR activation in relation to acetylcholine (ACh), which defines the maximal effect of the system (Fig. 3, e and f). Global analysis yields a significantly higher affinity of iper-6-naph for the dualsteric than for the allosteric binding mode (pK A and pK B are 7.40 Ϯ 0.06 and 6.84 Ϯ 0.06, respectively; p Ͻ 0.001, unpaired t test). Consequently, the fraction of active iper-6-naph⅐receptor complexes, f active (Fig. 3f), is significantly higher than the inactive fraction (79 versus 21%, p Ͻ 0.001).
Decrease of Agonist Efficacy Is Due to an Increase of Inactive Agonist-Receptor Complexes-We provide structural evidence that the ligand's preference for the purely allosteric binding mode determines its efficacy. We apply isox-6-naph, which differs from iper-6-naph only in its orthosteric moiety (for the molecular structures of all compounds see Fig. 4). Isox, the orthosteric building block of isox-6-naph, is a derivative of iperoxo displaying equal efficacy at muscarinic receptors albeit having a 100-fold lower affinity. Docking and MD simulations indicate that the topographies of iperoxo and isox in the orthosteric binding site of active M 2 AChR complexes are highly similar (Fig. 5, a and b). However, different hydrogen bond strengths were reported for furan and isoxazole compared with non-aromatic ring systems (38,39). Therefore, we expected a higher occurrence of geometries related to hydrogen bonds in the simulation of the iperoxo-bound M 2 AChR complex. To unveil these differences, we generated three-dimensional interaction point densities (dynophores) by MD-based sampling (Fig. 5c). Following this approach, we find that iperoxo forms a hydrogen bond to the orthosteric epitope Asn-404 6.52 in nearly all frames of the trajectory (98%), whereas isox shows this interaction in only 87% of the frames. Further analyses of hydrogen bond distances and angle frequencies confirm this observation, showing nearly ideal geometry for iperoxo, whereas for isox more deviations are observed (Fig. 5c). All other interactions show no difference between iperoxo and isox (Fig. 5b). Docking and MD simulations of M 2 AChR complexes bound to isox-6-naph in the dualsteric mode reveal overall homology to iper-6-naph (Fig. 5d). However, the dynophore analysis of iper-6-naph and isox-6-naph in their dualsteric binding modes reveals a similar occurrence of hydrogen bonds (98 versus 90%) as for the orthosteric building blocks iperoxo and isox (98 versus 87%). Hence, the topographical differences between iperoxo and isox (Fig. 5c) are preserved in the topographies of their dualsteric ligands iper-6-naph and isox-6-naph, respectively. In line with these structural data, global analysis of [ 3 H]NMS equilibrium binding and [ 35 S]GTP␥S binding data reveals a 100-fold lower affinity of isox-6-naph than iper-6-naph in their dualsteric binding modes as reflected by pK A (isox-6-naph) ϭ 5.12 Ϯ 0.15 and pK A (iper-6-naph) ϭ 7.40 Ϯ 0.06, respectively (Fig. 5f).
Isox-6-naph has recently been shown to be a less efficacious partial agonist than iper-6-naph (28). MD simulations of isox-6-naph in its allosteric binding mode reveal an unexpected intramolecular interaction. The isox moiety forms ainteraction with the allosteric moiety 6-naph that is not possible in the purely allosteric binding mode of iper-6-naph. This stabilizes an intramolecular conformation of isox-6-naph that is likely to bind with high affinity to the allosteric vestibule (Fig.  5e). In fact, the ammonium group of the isox moiety forms a charge interaction with Glu-175 EL2 . The other ammonium group displays cationinteractions with Tyr-80 2.60 and Tyr-426 7.39 . The triple bond lies opposite Tyr-177 EL2 , and the two methyl groups of the small allosteric linker chain show lipophilic contacts to Ile-178 EL2 and Tyr-426 7.39 . The naphthalene ring is embedded in a pocket lined by Phe-181 EL2 , Tyr-177 EL2 , and Thr-187 5.39 . Additionally, Ile-178 EL2 is able to form a hydrogen bond to one of the carbonyl groups of the allosteric moiety (Fig. 5e). This allosteric network suggests an optimal fit of isox-6-naph to the allosteric binding site. Indeed, global pharmacological analysis (Fig. 5f) yields a significantly higher affinity of isox-6-naph for the inactive than for the active pose (pK B and pK A are 7.11 Ϯ 0.15 and 5.12 Ϯ 0.15, respectively; p Ͻ 0.001, unpaired t test). This results in an almost exclusive allosteric binding topography of isox-6-naph (99 versus 1% for the inactive and active fraction, respectively; p Ͻ 0.001, unpaired t test).
Overall, isox-6-naph binds to the allosteric binding site with higher affinity than to the orthosteric binding site. Isox-6naph's preference for the purely allosteric binding mode explains its lower efficacy compared with iper-6-naph (Fig. 5f).
Increase of Agonist Efficacy Is Due to a Decrease of Inactive Agonist-Receptor Complexes-Iper-8-naph is a dualsteric ligand in which the orthosteric and allosteric building blocks are connected via an octamethylene linker. MD simulations of iper-6-naph (Fig. 3c) and isox-6-naph (Fig. 5e) in their purely allosteric binding modes imply that a longer linker would hamper a purely allosteric binding mode of iper-8-naph. Indeed, MD simulations of allosterically bound iper-8-naph⅐M 2 AChR complexes show a suboptimal topography of iper-8-naph in the allosteric binding site (Fig. 6a). Although iper-8-naph shares some key interactions with iper-6-naph (Fig. 3c) and isox-6-naph (Fig. 5e), its overall, purely allosteric topography is distinct (Fig. 6a). Moreover, its orientation in the allosteric binding site appears to be variable. Heavy atoms of iper-8-naph show a pronounced deviation of its docking pose during the simulation time. In contrast, iper-6-naph and even more so isox-6-naph remain in their orientation (Fig. 6b). However, MD simulations of dualsterically bound iper-8-naph⅐M 2 AChR complexes show homology to the dualsteric poses of iper-6-naph and isox-6naph. Global analysis of pharmacological data reveals that first a higher affinity of iper-8-naph for the dualsteric, active than for the allosteric, inactive pose (pK A and pK B are 8.50 Ϯ 0.06 and 6.87 Ϯ 0.06, respectively; p Ͻ 0.001, unpaired t test), second the fraction of active receptors is significantly greater than the fraction of inactive receptors (98 versus 2%; p Ͻ 0.001, unpaired t test), and third the efficacy of iper-8-naph resembles that of a full agonist (Fig. 6, c and d). Overall, a dualsteric ligand's preference for the dualsteric binding mode increases its efficacy. . Decrease of agonist efficacy is due to the preference for the purely allosteric binding mode. a, comparison of iperoxo in its co-crystallized conformation (dark gray) with isox from an MD simulation (light gray). The yellow spheres indicate lipophilic contacts, the red sphere indicates a hydrogen bond acceptor, and the positively charged center is shown as a blue star. b, similar pharmacophoric features of iperoxo (above) and isox (below) with key residues for ligand binding. The yellow spheres indicate lipophilic contacts, the red sphere indicates a hydrogen bond acceptor, and the positively charged center is shown as a blue star. c, the combination of three-dimensional pharmacophores and MD simulations led to the new concept of dynophores (dynamic pharmacophores) that are able to reflect time-dependent changes in the interaction pattern of ligands. The yellow cloud indicates lipophilic contacts, the red cloud indicates a hydrogen bond acceptor, and the positively charged center is shown as a blue cloud. Whereas the positively charged center and the lipophilic contacts could be observed in all frames of the MD simulation, the hydrogen bond acceptor feature was present in almost all frames (98%) for iperoxo and iper-6-naph but only in 87% of all frames for isox and 90% for isox-6-naph. Below, relative frequencies of distance and angle of the hydrogen bond acceptor. Curves for iperoxo and iper-6-naph are shown in black; curves for isox and isox-6-naph are shown in gray. d, the dualsteric binding modes of iper-6-naph (dark gray) and isox-6-naph (light gray) are highly similar. e, three-dimensional pharmacophore analyses of isox-6-naph in the purely allosteric binding mode. Note that isox-6-naph forms an intramolecularstacking interaction between the isoxazole moiety and the allosteric ring system. f, effects of isox-6-naph on [ 3 H]NMS equilibrium binding in HEPES buffer (orange curve) and G protein activation (green curve).
Rational Design of a Full Agonist with Exclusive Dualsteric Binding Topography-Dynamic ligand binding to M 2 AChRs is characterized by a binding mode ensemble comprising active and inactive ligand⅐receptor complexes bound to the same ligand in a dualsteric and an allosteric binding mode, respectively (Fig. 3). Both ligand⅐receptor complexes are in equilibrium, the position of which is sensitive to the affinity to both the orthosteric and allosteric binding sites (Figs. 1, 5, 6, and 7). This implies that it should be possible to design a putatively dual- steric ligand that does not bind to the receptor in the purely allosteric mode anymore and hence should display full agonism. Based on our structural model, we estimated the very extreme distance within the allosteric vestibule to be 15.7 Å, measured between Thr-84 2.65 and Thr-187 5.39 (Fig. 8b). Hereupon, we designed and synthesized the ligand iper-rigid-naph (for chemical synthesis, see supplemental Fig. 1). Iper-rigidnaph is a derivative of iper-6-naph comprising a linker that is longer (16.0 Å) and more rigid than the hexamethylene and octamethylene spacers (Fig. 8b). Attempts to dock iper-rigidnaph only into the allosteric vestibule of inactive M 2 AChRs failed due to multiple clashes with allosteric residues (Fig. 8b). In contrast, molecular docking and MD simulations of active M 2 AChR complexes indicate that iper-rigid-naph keeps the dualsteric binding mode, resulting in full receptor activation (Fig. 8a). The suggested loss of an allosteric binding pose of iper-rigid-naph is corroborated by pharmacological experiments (Fig. 8c). Global analysis of [ 3 H]NMS and [ 35 S]GTP␥S binding data indicates that iper-rigid-naph exclusively stabilizes active M 2 AChRs (Fig. 8d). In line with this, iper-rigid-naph displays full agonism (Fig. 8d).

Discussion
GPCRs are highly flexible membrane proteins (43,44), adopting a multitude of distinct conformations (inactive and active) even in the absence of ligands (19 -22, 24). Agonist binding leads to a shift of the receptor's conformational equilibrium toward active states. The equilibrium shift, however, is not quantitative. Even in the presence of a G protein mimetic, agonist-occupied receptors can adopt inactive states (21).
We have studied the existence, molecular details, and function of such inactive agonist⅐receptor complexes. Bitopic partial agonists for M 2 AChRs adopt a ligand binding ensemble consisting of a dualsteric and a purely allosteric binding mode (Fig.  3). Dualsterically bound agonist⅐receptor complexes induce G protein activation, whereas purely allosterically bound agonist⅐receptor complexes are inactive (Fig. 3). Multiple lines of evidence imply that both states reside in dynamic equilibrium (Figs. 1, 5, 6, and 7). First, small changes in the chemical structure of agonists (e.g. the orthosteric building block and the linker chain) strongly influence the affinities for either the active or the inactive receptor conformation (Figs. 5 and 6). Second, an allosteric triple mutant disrupts the purely allosteric binding topography of the agonists and thus fully rescues agonist efficacy (Fig. 7). Third, a putatively dualsteric agonist (iperrigid-naph), which was designed based on our findings, adopts only the dualsteric binding pose and hence displays full agonism (Fig. 8). Our data show that the population size of inactive agonist⅐receptor complexes determines overall agonist efficacy.
In addition to these experimental observations, Onaran and Costa (45) and others (46,47) have devised probabilistic models to theoretically describe a multitude of different receptor states called "receptor ensembles." Using probability partition functions, it can be shown that ligand binding to receptors changes the frequency distribution of receptor states. This theory directly implies that a ligand would have different affinities for different receptor conformations (47). Interestingly, probabilistic models do not distinguish between the chemical natures FIGURE 8. Rational design of a full agonist with exclusive dualsteric binding topography. a, transmembrane view of a representative conformation of the M 2 AChR in complex with iper-rigid-naph in the dualsteric binding mode taken from a 50-ns MD simulation. Agonistic moieties are shown with a green surface; antagonistic moieties are shown with a cyan surface. b, iper-rigid-naph cannot bind in a purely allosteric mode like iper-6-naph (Fig. 3c) because of the rigidified linker that is not able to enter the allosteric binding site. The maximum dilatation of the allosteric vestibule is shown as the distance between Thr-84 2.65 and Thr-187 5.39 (lower green line in b). This conformational constraint allows iper-rigid-naph to only bind in the dualsteric mode. c, effects of iper-rigid-naph on [ 3 H]NMS equilibrium binding in HEPES buffer (orange curve) and iper-rigid-naph-induced G protein activation (green curve). [ 3 H]NMS and [ 35 S]GTP␥S binding are plotted on the left and right ordinates, respectively. Data are mean Ϯ S.E. from five to nine independent experiments conducted in triplicate. d, the fractional population of active M 2 AChRs (f active ) stabilized by iper-rigid-naph (orange bar) and its maximal efficacy (E max ; green bar) relative to ACh were plotted on the left and right ordinates, respectively. f active was retrieved by fitting all data points in c globally to an operational model of agonism for dynamic ligands (see "Experimental Procedures"). Error bars represent S.E. of ligands or the binding sites that these ligands target preferentially. Hence, in principal, every ligand for GPCRs should have a multitude of microaffinities for different conformational states of the receptor. However, it is not known whether different agonist affinities go along with different agonist binding modes. To address this question experimentally, we have studied a special case of GPCR agonism because the rather extreme molecular nature of dualsteric agonists allows identifying multiple agonist binding modes with currently available techniques. Identifying multiple binding topographies of purely orthosteric agonists by crystallography appears to be technically challenging as less preferred or energetically less stable agonist⅐receptor complexes may not be retrieved easily in crystallization trials. Moreover, subtle differences in the agonist topography might exceed the resolution currently achieved. Although structural data of agonists bound to inactive state receptors have been reported (48,49), it is not known whether the agonists in these structures would adopt a different binding mode in the active state receptors. Furthermore, conventional molecular docking experiments seem unapt to describe such dynamic phenomena because they aim at predicting a single, most favorable conformation of the protein⅐ligand complex.
However, it is tempting to hypothesize that ligand binding ensembles indeed represent a more general mechanism for partial agonism. In agreement with probabilistic models of receptor ensembles (45,46), the conclusions drawn from this study should not be limited to bitopic agonists. In particular, our data suggest that it is of worth to study the existence of multiple binding modes of agonists that target only the orthosteric binding site. Our mechanistic model is supported by previously published data: the endogenous muscarinic AChR agonist ACh has been suggested to adopt both "productive" (i.e. signalingcompetent) and "non-productive" binding modes in the orthosteric site of M 1 AChRs (50). Moreover, for ligand-activated nuclear receptors, ligand binding ensembles have been crystallized (51)(52)(53)(54)(55) and might hence appear as a general principle in protein-ligand interactions. In the case of GPCRs, biophysical techniques, single molecule studies, and the dynophore approach may be helpful to address whether purely orthosteric agonists can form a ligand binding ensemble consisting of active and inactive agonist⅐receptor complexes.
Taken together, the combination of pharmacological experiments and computational simulations based on the inactive and active M 2 AChR crystal structures have led to a molecular description of an agonist binding ensemble. This is a proof of principle that active and inactive agonist⅐receptor complexes can be bound by an agonist in different binding modes. The concept of ligand binding ensembles may be applicable to other GPCRs and receptor classes.

Experimental Procedures
Materials, Chemical Probes, and Buffers-All cell culture media and supplements were purchased form Sigma-Aldrich and Invitrogen. All buffer reagents were from Sigma-Aldrich, Grüssing GmbH Analytika (Filsum, Germany), and Merck Labor und Chemie Vertrieb GmbH (Bruchsal, Germany).
Site-directed Mutagenesis and Generation of Stable Cell Lines-The allosteric triple mutant hM 2 Y177A,W422A,T423A was generated by site-directed mutagenesis (QuikChange kit (Stratagene, La Jolla, CA)) following the manufacturer's instructions. The cDNA of the allosteric double mutant hM 2 Y177A,W422A in pcDNA5/FRT was used as the template. The forward and reverse primers for induction of the T423A mutation by PCR were as follows (the underlined nucleotides indicate the mutated triplet): 5Ј-CCCCAACACTGTGGCGGCAATTG-GTTACTGGC-3Ј (forward) and 5Ј-GCCAGTAACCAATT-GCCGCCACAGTGTT GGGG-3Ј (reverse). The resulting construct was verified by sequencing. The stable Chinese hamster ovary (CHO)-Flp-In TM cell line expressing the triple mutant was generated according to the manufacturer's instructions (Flp-In system, Invitrogen).
Membrane Preparation-When stably transfected CHO cells had reached about 80 -90% confluence, the medium was aspirated and replaced by fresh medium containing 5 mM sodium butyrate, and cells were incubated for another 16 -18 h before membrane preparation. Medium was aspirated, and 2.4 ml of ice-cold harvesting buffer (20 mM HEPES, 10 mM Na 2 EDTA, pH 7.4) was added to the cells. Cells were detached from the culture dish using a cell scraper (Sarstedt, Newton, NC). The cell suspension was shredded twice for 20 s at level 6 using a Polytron homogenizer, and the resulting cell fragments were centrifuged at 40,000 ϫ g for 10 min (4°C). The supernatant was aspirated, and the pellet was resuspended in 15 ml of ice-cold centrifugation buffer (20 mM HEPES, 0.1 mM Na 2 EDTA, pH 7.4) before centrifugation in conditions mentioned above. This step was repeated once, and the remaining pellet was resuspended in ice-cold HEPES buffer (12.5 mM HEPES, 12.5 mM MgCl 2 , 125 mM NaCl, pH 7.4), quickly frozen, and stored in aliquots at Ϫ80°C. The protein concentration was determined by the method of Lowry (71).
[ 3 H]NMS Equilibrium Binding-CHO-hM 2 wild type and triple mutant membrane homogenates (15-20 g/ml) were incubated in binding buffer A or B supplemented with GDP (100 M), the radioligand [ 3 H]NMS (0.2 nM/well), and the test compounds at different concentrations for 18 h at 30°C. The incubation time was calculated as described previously (41). Atropine (10 M) was used to determine unspecific binding. The assay was terminated by rapid vacuum filtration. All experiments were carried out in a 96-well microtiter plate (Thermo Scientific Abgene, Germany) in a final volume of 500 l/well. Membrane-bound radioactivity was separated from non-bound radioligand by vacuum filtration using the TomTec filtration machine. The filter mats (Printed Filtermat A for use with 1450 MicroBeta TM , glass fiber filter, PerkinElmer Life Sciences) were then washed twice with ice-cold distilled water, further dried for 2.5 min in a microwave, and covered with melt-on scintillator sheets (MeltiLex TM melt-on scintillator sheets, PerkinElmer Life Sciences) on a heating block. The filter mats were then transferred into a plastic sample bag (Sample Bag for MicroBeta TM , PerkinElmer Life Sciences) and placed into a reading cassette. Finally, a solid scintillation counter (Beckman Instruments, Palo Alto, CA) was used to quantify the radioactivity on the filter mats.
[ 35 and where "Basal" and E max represent the unstimulated (in the absence of agonist) and maximally ACh-stimulated [ 35 S]GTP␥S accumulation (as a surrogate for the maximum response of the system), respectively. EC 50 is the concentration of the reference ligand ACh that gives the half-maximal response of the system. [X] is the concentration of the reference ligand ACh.
[AB] is the concentration of the dynamic ligand AB. n is the slope of the curve. K obs represents the observed equilibrium dissociation constant of the dynamic ligand. The global fit was set up using Equations 1, 2, and 5 with its subvariables. For fitting, the following parameters were constrained to the indicated values: Basal ϭ 0, log EC 50 to that of ACh determined separately for wild type (log EC 50 ϭ Ϫ7.27) and triple mutant receptors (log EC 50 ϭ Ϫ5.94), B 0 ϭ 100; log[L] ϭ Ϫ9.7, and log K L ϭ Ϫ9.28. ␣Ј was constrained to the cooperativity of the respective allosteric fragment with [ 3 H]NMS. E max , K obs , and n were shared among all data sets. This global fit yields three parameters, i.e. dyn , K obs , and f inactive . From those, the remaining parameters can be calculated consequently as follows. and Error propagation was accounted for by using the following equation.
S.E.ϭ ͱ(S.E. 1 ) 2 ϩ(S.E. 2 ) 2 (Eq. 12) Docking, All-atom Molecular Dynamics Simulations, and Three-dimensional Pharmacophore Analysis-In the presented study, different modeling techniques, such as docking, all-atom MD simulations, and both static and dynamic three-dimensional pharmacophore analyses, were successfully combined (36). All protein-ligand docking experiments reported in this study were carried out with the Cambridge Crystallographic Data Centre's software GOLD version 5.1 (61). Prior to the docking experiments, ligand conformations were generated by CORINA 3.0 (62). Antagonistic ligand poses were obtained by docking into the inactive M 2 AChR crystal structure (Protein Data Bank code 3UON) (15), and agonistic ligand poses were obtained by docking into the active M 2 AChR crystal structure (Protein Data Bank code 4MQT) (16). Due to the superagonistic properties of the co-crystallized ligand iperoxo, the tyrosine lid (Tyr-104 3.33 , Tyr-403 6.51 , and Tyr-426 7.39 ) strongly separates the orthosteric and the allosteric binding sites in the active crystal structure. Therefore, this tyrosine lid was remodeled by a side chain sampling based on MD simulations of the aporeceptor and the active crystal structure. The preparation of protein structures was performed using Molecular Operating Environment (MOE; 2014.09, Chemical Computing Group Inc.). All ligands and water molecules were removed, and correct protonation states were assigned. All residues of the inner core region and the extracellular domains were defined as potential binding sites (10 Å around the co-crystallized ligands; Protein Data Bank code 4MQT). Default settings were used for ligand docking, and GoldScore served as the scoring function. The obtained docking poses and receptor-ligand interactions were analyzed using LigandScout 3.1 (63, 64) using a three-dimensional pharmacophore approach.
All MD simulations described in this study were performed using Desmond 3.2 (65,66). An orthorhombic box was used to build the model systems with periodic boundary conditions in an isothermal-isobaric ensemble with a constant number of particles. The system temperature was kept at 300 K, and the pressure was kept at atmospheric pressure. The definition of transmembrane regions was taken from the OPM database (67). The receptor structures were embedded in a preequilibrated palmitoyloleoylphosphatidylcholine membrane (bilayer) and solvated with simple point charge water and 0.15 M NaCl. All other parameters were set on default values. Each simulation consisted of an equilibration run of 4.52 ns followed by a production run of 50 ns. The simulations were carried out on the Soroban computer cluster (Freie Universität Berlin) by using 24 central processing units. The obtained trajectories were analyzed with the software VMD (68) and LigandScout 3.1 (63,64).
Dynamic Three-dimensional Pharmacophores (Dynophores)-For the comparison of iperoxo and isox as well as their related dualsteric ligands iper-6-naph and isox-6-naph, respectively, a novel approach termed dynophore was developed that combines static three-dimensional pharmacophores with MD-based conformational sampling. Unlike previous applications of gathering pharmacophore information from molecular dynamics (33,69), this new implementation works in a fully automated way: dynophore groups interaction points (such as hydrogen bonds, charges, and lipophilic contacts) of each tra-jectory frame according to their type and ligand atoms involved. All feature groups are graphically represented by three-dimensional volumetric feature density clouds, which are statistically characterized by occurrence frequency and interaction patterns with the protein. The dynophore algorithm was implemented within the LigandScout framework (63,64,70).
Statistical Analysis-Data are shown as mean Ϯ S.E. for n observations. A significant difference between two distinct values was tested using an unpaired t test. Comparisons of groups were performed using one-way analysis of variance with Bonferroni's multiple comparison test.