Structure and Dynamics of the M3 Muscarinic Acetylcholine Receptor

Acetylcholine (ACh), the first neurotransmitter to be identified1, exerts many of its physiological actions via activation of a family of G protein-coupled receptors (GPCRs) known as muscarinic ACh receptors (mAChRs). Although the five mAChR subtypes (M1-M5) share a high degree of sequence homology, they show pronounced differences in G protein coupling preference and the physiological responses they mediate.2–4 Unfortunately, despite decades of effort, no therapeutic agents endowed with clear mAChR subtype selectivity have been developed to exploit these differences.5–6 We describe here the structure of the Gq/11-coupled M3 mAChR bound to the bronchodilator drug tiotropium and identify the binding mode for this clinically important drug. This structure, together with that of the Gi/o-coupled M2 receptor, offers new possibilities for the design of mAChR subtype-selective ligands. Importantly, the M3 receptor structure allows the first structural comparison between two members of a mammalian GPCR subfamily displaying different G-protein coupling selectivities. Furthermore, molecular dynamics simulations suggest that tiotropium binds transiently to an allosteric site en route to the binding pocket of both receptors. These simulations offer a structural view of an allosteric binding mode for an orthosteric GPCR ligand and raise additional opportunities for the design of ligands with different affinities or binding kinetics for different mAChR subtypes. Our findings not only offer new insights into the structure and function of one of the most important GPCR families, but may also facilitate the design of improved therapeutics targeting these critical receptors.

differences. [5][6] We describe here the structure of the G q/11 -coupled M 3 mAChR bound to the bronchodilator drug tiotropium and identify the binding mode for this clinically important drug. This structure, together with that of the G i/o -coupled M 2 receptor, offers new possibilities for the design of mAChR subtype-selective ligands. Importantly, the M 3 receptor structure allows the first structural comparison between two members of a mammalian GPCR subfamily displaying different G-protein coupling selectivities. Furthermore, molecular dynamics simulations suggest that tiotropium binds transiently to an allosteric site en route to the binding pocket of both receptors. These simulations offer a structural view of an allosteric binding mode for an orthosteric GPCR ligand and raise additional opportunities for the design of ligands with different affinities or binding kinetics for different mAChR subtypes. Our findings not only offer new insights into the structure and function of one of the most important GPCR families, but may also facilitate the design of improved therapeutics targeting these critical receptors.
The mAChR family consists of five subtypes, M 1 -M 5 , which can be subdivided into two major classes (Fig 1a). The M 1 , M 3 , and M 5 receptors show selectivity for G proteins of the G q/11 family, whereas the M 2 and M 4 receptors preferentially couple to G i/o -type G proteins [2][3][4] . The development of small molecule ligands that can selectively act on specific mAChR subtypes has proven extremely challenging, primarily due to the high degree of sequence similarity in the transmembrane (TM) core of these receptors [2][3][4] . More recently, considerable progress has been made in targeting drugs to non-classical (allosteric) binding sites of certain mAChR subtypes 5 .
Within the mAChR family, the M 3 subtype mediates many important physiological functions, including smooth muscle contraction and glandular secretion 3,4,6-9 . Central M 3 receptors have also been implicated in the regulation of food intake 7 , learning and memory 8 , and the proper development of the anterior pituitary gland 9 . Selective drugs targeted at this receptor subtype may prove clinically useful 4,6-9 , and non-selective muscarinic ligands are already widely used in current practice.
Due to the profound physiological importance of the M 3 receptor and its longstanding role as a model system for understanding GPCR function 3,10 , we used the T4 lysozyme (T4L) fusion protein strategy 11 to obtain crystals of Rattus norvegicus M 3 receptor-T4L fusion protein ( Supplementary Fig. 1) by lipidic cubic phase crystallization. Diffraction data from more than 70 crystals were merged to create a data set to 3.4 A resolution and to solve the structure by molecular replacement. The M 3 receptor structure, together with that of the M 2 receptor 12 , affords for the first time an opportunity to compare two closely related mammalian receptors with divergent G protein coupling selectivities.
The overall structure of the M 3 receptor is similar to that of M 2 (Fig. 1b-d). Surprisingly, structural conservation includes intracellular loops (ICLs) 1 and 2, and extracellular loops (ECLs) 1-3, which share highly similar overall folds despite low sequence conservation (Fig. 1f). Like the M 2 receptor, the M 3 receptor exhibits unique mAChR features, including a large extracellular vestibule as part of an extended hydrophilic channel containing the orthosteric binding site (Fig. 1e). Also like M 2 , the M 3 receptor features a pronounced outward bend at the extracellular end of TM4 ( Fig. 1d; Supplementary Fig. 2b). This bend, not seen in any other GPCR family crystallized to date, is stabilized by a hydrogen bond from the Q207 (Q163 in M 2 ) side chain to the L204 backbone peptide carbonyl ( Supplementary Fig. 2b). This bond is part of a polar interaction network involving four residues absolutely conserved within the mAChR family, suggesting that this unusual feature is important to mAChRs in general. Indeed, mutagenesis of Q207 in M 3 impaired both ligand binding and receptor activation 13 .
The M 3 receptor was crystallized in complex with tiotropium (Spiriva), a potent muscarinic inverse agonist 14,15 used clinically for the treatment of chronic obstructive pulmonary disease. The M 2 receptor was crystallized in complex with R-(-)-3-quinuclidinyl benzilate (QNB) which, like tiotropium, is a non-subtype-selective mAChR blocker 14,16 . The two ligands bind in remarkably similar poses (Fig. 2b), and it is likely that this pose represents a conserved binding mode for structurally similar anticholinergics. In the M 3 receptor, as in M 2 , the ligand is deeply buried within the TM receptor core (Fig. 2a, d) and is covered by a lid comprised of three conserved tyrosines, Y148 3.33 , Y506 6.51 , and Y529 7.39 ( Fig. 2a; superscripts indicate Ballesteros-Weinstein numbers 17 ). The ligand is almost completely occluded from solvent and engages in extensive hydrophobic contacts with the receptor. A pair of hydrogen bonds are formed from N507 6.52 to the ligand carbonyl and hydroxyl, while D147 3.32 interacts with the ligand amine.
Reflecting the difficulty in developing subtype-selective orthosteric ligands, the residues forming the orthosteric binding pocket are absolutely conserved among the five mAChR subtypes (Fig. 1f). However, this conservation at the amino acid level does not preclude the existence of differences in the three-dimensional architecture of the orthosteric site between the different mAChR subtypes. In fact, comparison of the structures of the M 3 and M 2 receptor ligand binding sites reveals for the first time structural divergences that might be exploited in the development of subtype-selective ligands.
One such difference derives from the replacement of Phe181 in ECL2 of M 2 with Leu225 in M 3 (this residue is leucine in all mAChRs except M 2 ). This creates a pocket in M 3 not found in M 2 (Fig. 2c,d). A second difference is a 2.8 A shift of Tyr529 7.39 relative to the position of the corresponding M 2 residue (Tyr426; Fig. 2e). This feature may derive from a difference in the identity of the residue in position 2.61 (Phe124 in M 3 and Tyr80 in M 2 ; Fig. 2f). This residue interacts directly with TM7, influencing the position of this helix and the residues within it, including Tyr529 7.39 . Notably, the residue at position 2.61 is not a part of the orthosteric binding pocket, but is positioned near a probable allosteric binding site 12 . Since tiotropium and QNB are structurally similar but not identical, the observed binding site differences must be interpreted with some degree of caution. However, sitedirected mutagenesis studies with M 1 and M 3 receptors support the concept that the residue at position 2.61 plays a role in receptor activation 18,19 , and ligand binding selectivity 20 . This site does not appear to play a role in determining antagonist dissociation rates, since mutation of M 3 F 2.61 to tyrosine or of M 2 Y 2.61 to phenylalanine had no effect on dissociation rates for [ 3

H]NMS or [ 3 H]QNB.
We used molecular dynamics simulations to characterize the pathway by which tiotropium binds to and dissociates from the M 2 and M 3 receptors. Similar techniques have previously been shown to correctly predict crystallographic ligand binding poses and kinetics in studies of β-adrenergic receptors 21 . In both the M 2 and M 3 receptors, our simulations indicate that as tiotropium binds to or dissociates from the receptor, it pauses at an alternative binding site in the extracellular vestibule (Fig. 3, Supplementary Fig. 3). Intriguingly, this site corresponds to an allosteric site that has been previously identified by mutagenesis 12 , a finding consistent with pharmacological studies showing that orthosteric ligands can act as allosteric modulators at the M 2 receptor 22 .
Tiotropium adopts different preferred allosteric binding poses in M 2 and M 3 (Fig. 3d, Supplementary Fig. 4). These metastable binding poses, which appear independently in both binding and dissociation simulations, may represent the first structural view of a clinically used "orthosteric" GPCR ligand binding to an experimentally validated allosteric site. Conceivably, therapeutic molecules could be rationally engineered to act independently as both allosteric and orthosteric ligands (in contrast to previously described bitopic ligands that bind at both sites simultaneously 23 ). Tiotropium dissociates from M 3 receptors more slowly than from M 2 receptors, a phenomenon thought to provide clinically important "kinetic selectivity" of this drug for M 3 receptors despite similar equilibrium binding affinities for both subtypes 14 . In simulations with tiotropium bound, the portion of ECL2 nearest the binding pocket proved more mobile in M 2 than in M 3 ( Supplementary Fig. 5), likely due to multiple sequence differences between the two receptor subtypes. This increased mobility disrupts a hydrophobic cluster involving a thiophene ring of tiotropium, the ECL2 residue Phe181(M 2 )/Leu225(M 3 ), and Tyr 3.33 , facilitating movement of Phe181/ Leu225 away from the orthosteric site and rotation of Tyr 3.33 toward TM4. In simulations of ligand dissociation, such motions clear a path for tiotropium's egress from the orthosteric site to the extracellular vestibule. The increased mobility of ECL2 in M 2 thus appears to facilitate tiotropium's traversal of the largest energetic barrier on the binding/dissociation pathway (Fig. 3c). Experimental measurements with wild-type and mutant receptors (M 3 L225F and M 2 F181L) suggest that the Leu225/Phe181 sequence difference alone is insufficient to explain the difference in off-rates, (for practical reasons these measurements were performed with QNB rather than tiotropium; see Online Methods).
One of the most interesting features of the M 2 and M 3 receptors is the fact that the two highly similar receptors display pronounced differences in G-protein coupling specificity. For this reason, the M 2 /M 3 receptor pair has long served as an excellent model system to identify features contributing to GPCR-G protein coupling selectivity 3 . Since no simple sequence elements have been identified as general determinants of coupling specificity across GPCR families 24 , it is likely that recognition depends on features such as overall conformation in addition to specific inter-residue contacts.
The M 2 and M 3 receptor structures show a significant difference in the position of the cytoplasmic end of TM5 and of ICL2 (Fig. 4a, b). The highly conserved tyrosine residue at position 5.58 (M 3 Tyr250 5.58 , M 2 Tyr206 5.58 ) shows a clear deviation between the two receptors, pointing toward the core of the protein in M 2 , and away from the receptor toward the surrounding lipid bilayer in M 3 . Interestingly, mutagenesis studies have identified a tetrad of residues ('AALS' in M 3 , 'VTIL' in M 2 ) located on the cytoplasmic end of TM6 that are critical in determining G protein coupling selectivity 25,26 . In both structures, these residues interact directly with TM5 (Fig. 4a), and in the β 2 adrenergic receptor-G s complex 27 two of the four corresponding residues make contact with the carboxy terminal helix of Gα s . M 3 Tyr254 5.62 at the bottom of TM5 also plays a role in activation of G q/ 11 28 .
In the M 2 receptor structure, the corresponding residue (Ser210 5.62 ) is displaced by approximately 4 A relative to Tyr254 5.62 in M 3 (Fig. 4a).
When we compared the position of TM5 in the M 2 and M 3 receptors to that in other GPCR structures, we found that it is M 2 -like in all G i/o -coupled receptors, while the two mammalian G q/11 -coupled receptors solved to date exhibit another conformation (Fig. 4c, d).
An important caveat here is that these structures have been solved using the T4L fusion strategy, and we cannot completely exclude the possibility that this approach perturbs the conformation of TMs 5 and 6. However, in molecular dynamics simulations of M 2 and M 3 receptors without T4L, each of the receptors adopts a set of conformations that includes its own crystallographically observed conformation ( Supplementary Fig. 6, 7). These simulations suggest that the observed conformations are unlikely to be artifacts of the crystallization methodology, though the crystal structures likely represent only one conformation among many adopted by the receptors in a biological context.
The structure of the M 3 receptor, together with that of the M 2 receptor 12 , offers a unique opportunity to directly compare the structural properties of two members of a mammalian GPCR subfamily endowed with different G-protein coupling selectivities. Examination of the M 3 structure has provided the first structural evidence of differences between ligand binding sites of mAChR subtypes that could be exploited for the design of more selective therapeutics. Moreover, computational studies have identified a pathway by which the COPD drug tiotropium may bind to and dissociate from the M 3 receptor, offering the first structural view of an orthosteric GPCR ligand binding to an experimentally validated allosteric site. This information should facilitate the rational design of new muscarinic drugs exhibiting increased receptor subtype selectivity, potentially improving treatment for a wide variety of important clinical disorders.

Methods summary
The M 3 muscarinic receptor-T4 lysozyme fusion protein was expressed in Sf9 insect cells and purified by nickel affinity chromatography followed by FLAG antibody affinity chromatography and then size exclusion chromatography. It was crystallized using the lipidic cubic phase technique, and diffraction data were collected at the GM/CA-CAT beamline at the Advanced Photon Source at Argonne National Lab. The structure was solved by molecular replacement using merged data from 76 crystals. All-atom classical molecular dynamics simulations with explicitly represented lipids and water were performed using the CHARMM force field 29 on Anton 30 . Ligand-binding simulations included no artificial forces. Dissociation studies included a time-varying biasing term that gradually forces the ligand away from its crystallographic position, but not along any prespecified pathway or direction. Full details are provided in the online methods.

Expression and purification of M 3 muscarinic receptor
The wild-type M 3 mAChR contains several long, likely poorly ordered regions, including the extracellular amino terminal domain and the third intracellular loop, making it a challenging candidate for crystallographic studies. To alleviate this problem, the M 3 receptor from Rattus Norvegicus was modified to include a TEV protease recognition site in the amino terminus and a hexahistidine tag at the carboxy terminus. Moreover, the third intracellular loop (residues 260-481) was replaced with T4 lysozyme residues 1-161 in a manner described previously 11 , with two different fusions tested. These modifications are diagrammed in figure S1, which also shows the final crystallization construct.
The pharmacological properties of the construct were tested and compared to those of the wild-type receptor (Fig. S8, supplementary table 1; see below for methods details). Both constructs showed almost identical affinity for antagonists, while the crystallization construct (M 3 -crys) showed somewhat higher affinity for the agonist acetylcholine than the wild-type construct. A similar observation has been noted previously in the β 2 adrenergic receptor 11 . Studies with membranes prepared from transfected COS-7 cells showed that TEV cleavage of M 3 -crys (to remove most of the N-terminal tail) had no significant effect on ligand binding affinities (Fig. S9). Moreover, the wild-type receptor and M 3 -crys, either cleaved with TEV or left uncleaved, showed very similar [ 3 H]-QNB dissociation rate kinetics (Fig. S10). As expected, the crystallization construct failed to stimulate agonistdependent phosphoinositide hydrolysis in transfected COS cells (data not shown), likely because essential G protein interacting regions in intracellular loop 3 were omitted from the construct and also because the T4 lysozyme fusion protein sterically blocks G protein association.
The crystallization construct was expressed in Sf9 cells using the baculovirus system in the presence of 1 μM atropine. M 3 receptors expressed in Sf9 cells are known to exhibit functional and pharmacological properties similar to M 3 receptors expressed in mammalian cells 31 . Infection was performed at 4 × 10 6 cells per mL and flasks were shaken at 27 °C for 60 hr following infection.
Cells were harvested by centrifugation, then lysed by osmotic shock in the presence of 1 μM tiotropium bromide (obtained from W & J PharmaChem, Inc., Silver Spring, MD, USA), which was present in all subsequent buffers. Receptor was extracted from cells using a dounce homogenizer with a buffer of 0.75 M NaCl, 1% dodecyl maltoside (DDM), 0.03% cholesterol hemisuccinate (CHS), 30 mM HEPES pH 7.5, and 30% glycerol. Iodoacetamide (2 mg/mL) was added to block reactive cysteines at this stage. Nickel-NTA agarose was added to the solubilized receptor without prior centrifugation, stirred for two hr, and then washed in batch with 100 × g spins for 5 min each. Washed resin was poured into a glass column, and receptor was eluted in 0.1% DDM, 0.03% CHS, 20 mM HEPES pH 7.5, 0.75 M NaCl, and 250 mM imidazole.
Nickel resin-purified receptor was then loaded by gravity flow over anti-FLAG M1 affinity resin. Following extensive washing, detergent was gradually exchanged over 1.5 hr into a buffer in which DDM was replaced with 0.01% lauryl maltose neopentyl glycol (MNG), and the NaCl concentration was lowered to 100 mM. MNG has been shown to be more effective at stabilizing muscarinic receptors than DDM 32 . Receptor was eluted with 0.2 mg/mL FLAG peptide and 5 mM EDTA. TEV protease (1:10 w/w) was added and incubated with receptor for 1.5 hr at room temperature to remove the flexible amino terminal tail. Receptor was then separated from TEV by size exclusion chromatography (SEC) on a Sephadex S200 column (GE Healthcare) in a buffer of 0.01% MNG, 0.001% CHS, 100 mM NaCl, and 20 mM HEPES pH 7.5. Tiotropium was added to a final concentration of 10 μM following SEC. The resulting receptor preparation was pure and monomeric (Fig. S11). Purification of unliganded M 3 receptor was also possible by this procedure, but the resulting preparation was polydisperse and unsuitable for crystallographic study.

Crystallization and data collection
Purified M 3 receptor was concentrated to 60 mg/mL, then mixed with 1.5 parts by weight of a 10:1 mix of monoolein with cholesterol (Sigma) using the two syringe reconstitution method 33 . The resulting lipidic cubic phase mix was dispensed in 15 nL drops onto glass plates and overlaid with 600 nL precipitant solution using a Gryphon LCP robot (Art Robbins Instruments). Crystals grew after 2 -3 days in precipitant solution consisting of 27 -38 % PEG 300, 100 mM HEPES pH 7.5, 1 % (w/w) 1,2,3-heptanetriol, and 100 mM ammonium phosphate. Typical crystals are shown in figure S12.
Data collection was performed at Advanced Photon Source GM/CA-CAT beamlines 23ID-B and 23ID-D using a beam size of 5 or 10 microns for most crystals. Diffraction quality rapidly decayed following exposure, and wedges of typically 5 degrees were collected and merged from 76 crystals using HKL2000 34 . Diffraction quality ranged from 3-4 A in most cases, with strong anisotropy evident in many frames. Most crystals tested showed evidence of epitaxial twinning, though in most cases one of the two twins dominated the observed diffraction pattern, allowing processing as a single crystal. A more extensive discussion of the twinning is given below. Some contamination of diffraction measurements due to the twin-related reflections was unavoidable, leading to slightly poorer merging statistics than is typical for datasets collected from many small crystals (Supplementary table 2). Despite this, maps were generally of high quality and electron density was easily interpretable (Figs. S13, S14), in part due to the availability of non-crystallographic symmetry.
Analysis of <F>/<σF> along each of the three reciprocal space axes indicated that the diffraction was strong in two directions, and weak in the third direction, along the reciprocal space axis c* (Fig. S15). Using <F>/<σF> greater than 3 as a guideline suggested a resolution cutoff of better than 3.2 A along a* and b*, and of 4.0 A along c*. We therefore applied an ellipsoidal truncation along these limits, and then applied an overall spherical truncation at 3.4 A due to low completeness in higher resolution shells. Fortunately, 4-fold non-crystallographic symmetry (NCS) allowed for improved map quality with map sharpening followed by NCS averaging, largely alleviating the effects of anisotropic diffraction and epitaxial twinning to give highly interpretable maps (Fig. S13, S14) and allowing details of ligand recognition to be clearly resolved (Supplementary table 3).
The structure of the M 3 receptor was solved using the structure of the M 2 muscarinic receptor (companion manuscript) as the search model in Phaser 35 . The model was improved through iterative refinement in Phenix 36 and manual rebuilding in Coot guided by both NCS averaged and unaveraged maps. NCS restraints were applied in initial refinement stages, and omitted in final refinement cycles to account for differences between NCS-related copies. The quality of the resulting structure was assessed using MolProbity 37 , and figures were prepared using PyMOL 38 .

Epitaxial twinning
Crystals of the M 3 receptor showed hallmarks of epitaxial twinning, such as mixed sharp and split spots, poor indexing, and many unpredicted reflections in some frames. In some cases diffraction from two distinct lattices was clearly visible, with a small fraction of reflections exactly superimposed from both lattices (Fig. S16). In most cases one lattice dominated the diffraction pattern to such an extent that it could be easily processed as a single crystal. Intriguingly, the two indexing solutions were not equivalent cells but rather were two enantiomorphic P1 cells (supplementary table 2).
As one of these two cells gave significantly better diffraction data than the other, data processing and refinement were only pursued in this case. Within the asymmetric unit, two layers of receptors and two layers of T4 lysozyme are present, but each of these four layers exhibits a different lattice packing (Fig. S17, S18). The order in which these layers are stacked in the crystal defines a unique direction along c, the axis normal to the membrane plane. Since P1 is a polar space group, the positive direction along c is uniquely defined, and the two possible orientations of the stacked layers of membrane relative to the positive direction along c distinguish the two twin crystal forms.

Transient expression of wild-type and modified M 3 receptors in COS-7 cells, membrane preparation, and TEV treatment
COS-7 cells were cultured as described previously 39 . About 24 h prior to transfections, ~1 × 10 6 cells were seeded into 100 mm dishes. Cells were transfected with 4 μg/dish of receptor plasmid DNA using the Lipofectamine Plus kit (Invitrogen, Carlsbad, CA), according to the manufacturer's instructions. The mammalian expression plasmid coding for the wild type rat M 3 receptor has been described previously 40

[ 3 H]-QNB dissociation rate assays
[ 3 H]-QNB (PerkinElmer; specific activity: 50.5 Ci/mmol) dissociation rate assays were carried out as described previously 43 . Measurements were carried out at 37 o C in a total volume of 620 μl using a buffer consisting of 4 mM Na 2 HPO 4 and 1 mM KH 2 PO 4 (pH 7.4). Membranes prepared from transfected COS-7 cells (final protein concentration: 10 μg protein/ml) were prelabeled with 1 nM [ 3 H]-QNB for 30 min. Dissociation of the labeled ligand was initiated by the addition of atropine (final concentration: 3 μM). Incubations were terminated by filtration through GF/C Brandel fiber filters that had been pretreated with 0.1% polyethyleneimine, followed by two rinses with ice-cold distilled water. The amount of radioactivity that remained bound to the filters was determined by liquid scintillation spectrometry.

Molecular dynamics
In all simulations, the receptor was embedded in a hydrated lipid bilayer with all atoms, including those in the lipids and water, represented explicitly. Simulations were performed on Anton 30 , a special-purpose computer designed to accelerate standard molecular dynamics simulations by orders of magnitude.

System setup and simulation protocol
Simulations of the M 2 receptor were based on the crystal structure of the QNB-M 2 complex, and simulations of M 3 were based on the structure of the tiotropium-M 3 complex (chain A). These crystal structures were determined using a T4 lysozyme (T4L) fusion strategy, in which intracellular loop 3 (ICL3) of each receptor was replaced by T4L; the T4L sequence was omitted in our simulations. All chain termini were capped with neutral groups (acetyl and methylamide). Residues 6.31-6.33 near the intracellular end of TM6 were unresolved in the M 3 crystal structure, and residues 6.27-6.30 were resolved in an unstructured conformation packed against T4L. Residues 6.27-6.36 were modeled manually as a helical extension of TM6, with side chains then placed using Prime. Hydrogens were added to the crystal structures using Maestro (Schrodinger LLC, New York, NY), as described in previous work 44 . All titratable residues were left in the dominant protonation state at pH 7.0, except for Asp69 2.50 in M 2 and Asp114 2.50 in M 3 , which were protonated. Asp69 2.50 and Asp114 2.50 correspond to rhodopsin Asp83 2.50 , which is protonated during the entire photocycle 45 .
Prepared protein structures were inserted into an equilibrated POPC bilayer as described in previous work 46 . Sodium and chloride ions were added to neutralize the net charge of the system and to create a 150 mM solution.
Simulations of the M 3 receptor initially measured 80 × 80 × 87 A 3 and contained 163 lipid molecules, 26 sodium ions, 41 chloride ions, and approximately 9,897 water molecules, for a total of ~56,000 atoms. Simulations of the M 2 receptor initially measured 79 × 79 × 85 A 3 and contained 156 lipid molecules, 24 sodium ions, 35 chloride ions, and approximately 9,165 water molecules, for a total of ~53,000 atoms. To simulate M 2 with tiotropium bound, we removed the co-crystallized ligand, QNB, and docked in tiotropium using Glide (Schrodinger LLC, New York, NY).
All simulations were equilibrated using Anton in the NPT ensemble at 310 K (37 °C) and 1 bar with 5 kcal mol 1 A 2 harmonic position restraints applied to all non-hydrogen atoms of the protein and the ligand (except for the tiotropium-M 2 complex, where the ligand was unrestrained); these restraints were tapered off linearly over 50 ns. All bond lengths to hydrogen atoms were constrained using M-SHAKE 47 . A RESPA integrator 48 was used with a time step of 2 fs, and long-range electrostatics were computed every 6 fs. Production simulations were initiated from the final snapshot of the corresponding equilibration runs, with velocities sampled from the Boltzmann distribution at 310 K, using the same integration scheme, long-range electrostatics method, temperature, and pressure. Van der Waals and short-range electrostatic interactions were cut off at 13.5 A and long-range electrostatic interactions were computed using the k-space Gaussian Split Ewald method 49 with a 32 × 32 × 32 grid, σ = 3.33 A, and σ s = 2.33 A.

Spontaneous binding of tiotropium and acetylcholine
We performed simulations where tiotropium was placed arbitrarily in the bulk solvent (at least 40 A from the entrance to the extracellular vestibule) and allowed to diffuse freely until it associated spontaneously with the M 2 or M 3 receptor, following methodology described in previous work 21 . In these simulations (Supplementary table 4, conditions D and E), the cocrystallized ligand was removed and four tiotropium molecules were placed in the bulk solvent. A tiotropium molecule bound to the extracellular vestibule at least once in each simulation. In the longer simulations, tiotropium bound to and dissociated from the extracellular vestibule multiple times. Tiotropium assumed several different poses when bound to the extracellular vestibule of either M 2 or M 3 (Fig. S4). Tiotropium never entered the orthosteric binding pocket, presumably because the simulations were not of sufficient length.
The fact that tiotropium associated with and dissociated from the vestibule multiple times, but did not enter the binding pocket, suggests that tiotropium must traverse a larger energetic barrier to enter the binding pocket of M 2 or M 3 from the extracellular vestibule than to enter the vestibule from bulk solvent. This contrasts with earlier simulations on alprenolol binding to the β 2 -adrenergic receptor, in which the largest energetic barrier (by a small margin) was between the bulk solvent and the extracellular vestibule 21 . This difference likely reflects the fact that ligands must pass through a much tighter passageway to enter the binding pocket of the M 2 and M 3 receptors from the vestibule than is the case for the β 2 -adrenergic receptor. Tiotropium lost the majority of its hydration shell as it entered the vestibule (Fig. S19), as observed previously for ligands binding to β-adrenergic receptors 21 .
We followed a similar protocol in a simulation of the M 3 receptor in the presence of the agonist acetylcholine, a smaller molecule which might be expected to bind faster (Supplementary table 4, condition F). Indeed, an acetylcholine molecule bound in the orthosteric binding pocket after 9.5 μs and remained there for the remainder of the 25-μs simulation. Although acetylcholine quickly passed through the extracellular vestibule en route to the binding pocket, it did not exhibit metastable binding in the vestibule. Acetylcholine exhibited significant mobility in the binding pocket, likely reflecting the low affinity of the crystallized inactive state for agonists.

Forced dissociation of tiotropium
To identify the entire binding/dissociation pathway, we "pushed" tiotropium out of the binding pocket of both the M 2 and M 3 receptors 50,51 . Production simulations were initiated from configurations of the corresponding unbiased trajectory. These simulations employed a time-dependent harmonic biasing potential, U(t): where t is time, k is a force constant in units of kcal/mol/A 2 , d is the distance between the center-of-mass of the heavy atoms of tiotropium and the center-of-mass of the protein Cα atoms, and d 0 (t) varied linearly over 1.0 μs from 9.6 A to 33 A for M 2 and from 8.6 A to 32 A for M 3 . This biasing term does not impose any preferred direction of ligand exit. We performed seven such simulations for each of M 2 and M 3 ,with k = 5, starting from configurations extracted from the tiotropium bound simulations of M 2 and M 3. Each initial configuration was separated in time by 36 ns. Results were similar across all simulations.

Force field parameters
The CHARMM27 parameter set 29 with CMAP terms 52 and a recently introduced correction to charged side-chain electrostatics 53 was used for all protein molecules and salt ions in conjunction with the CHARMM TIP3P 54 water model and a modified CHARMM lipid force field 55 . Force field parameters for tiotropium and QNB were obtained from the CHARMM ParamChem web server 56 , version 0.9.1 beta. QNB was simulated in its protonated (ammonium) state. To evaluate the assigned partial charges assigned by ParamChem, we performed a quantum mechanical computation of the electrostatic potential at a collection of points surrounding each ligand (in vacuo at the HF/6-31G* level of theory using MOLPRO 57 ), and compared it to the potential generated by the assigned charges.

Analysis protocols
Trajectory snapshots, each containing a record of all atom positions at a particular instant in time, were saved every 180 ps during production simulations. Distance measurements were computed using the HiMach parallel analysis framework 58 . VMD 59 was used to visualize trajectories and to produce Fig. 3a and Fig. 3b.
To determine the most common vestibule-bound poses of tiotropium shown in Fig. S4, we performed a clustering analysis on the 14.2-μs spontaneous binding simulation of M 2 (Table  S4, condition D) and the 16.0-μs spontaneous binding simulation of M 3 (Table S4, condition E). We performed k-means clustering on the set of trajectory snapshots in which a tiotropium molecule was in the extracellular vestibule, using the positions of atoms indicated in Fig. S4c. Clusters representing highly similar poses were merged.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.    have been shown to be important coupling specificity determinants. b, ICL2 is also divergent between the two structures. Four residues previously implicated as specificity determinants 25 are shown, with residue numbers for M 2 followed by M 3 c, Plot of interhelical distances for crystallographically unique inactive GPCR structures published to date. Distances were measured between C α atoms of TM5 residue 5.62 and TM3 residue 3.54 (x-axis), and TM5 residue 5.62 and TM6 residue 6.37 (y-axis). GPCRs cluster by coupling specificity, although squid rhodopsin is an exception. GPCRs coupling preferentially to G i/o and those coupling to the homologous G protein G t cluster together. d, Structural alignment of mammalian G i/o -coupled and G q/11coupled receptor structures.