Skip to main content Accessibility help
×
Home
Hostname: page-component-59b7f5684b-2bkkj Total loading time: 1.337 Render date: 2022-10-04T23:13:25.682Z Has data issue: true Feature Flags: { "shouldUseShareProductTool": true, "shouldUseHypothesis": true, "isUnsiloEnabled": true, "useRatesEcommerce": false, "displayNetworkTab": true, "displayNetworkMapGraph": false, "useSa": true } hasContentIssue true

Accelerated molecular dynamics simulations of ligand binding to a muscarinic G-protein-coupled receptor

Published online by Cambridge University Press:  16 July 2015

Kalli Kappel
Affiliation:
Department of Chemistry and Biochemistry, University of California at San Diego, La Jolla, CA 92093, USA
Yinglong Miao*
Affiliation:
Howard Hughes Medical Institute, University of California at San Diego, La Jolla, CA 92093, USA
J. Andrew McCammon
Affiliation:
Department of Chemistry and Biochemistry, University of California at San Diego, La Jolla, CA 92093, USA Howard Hughes Medical Institute, University of California at San Diego, La Jolla, CA 92093, USA Department of Pharmacology, University of California at San Diego, La Jolla, CA 92093, USA
*
*Author for correspondence: Y. Miao, Howard Hughes Medical Institute, University of California at San Diego, La Jolla, CA 92093, USA. Tel.: 1-850-822-0255; Fax: 1-858-534-4974; Email: yimiao@ucsd.edu
Rights & Permissions[Opens in a new window]

Abstract

Elucidating the detailed process of ligand binding to a receptor is pharmaceutically important for identifying druggable binding sites. With the ability to provide atomistic detail, computational methods are well poised to study these processes. Here, accelerated molecular dynamics (aMD) is proposed to simulate processes of ligand binding to a G-protein-coupled receptor (GPCR), in this case the M3 muscarinic receptor, which is a target for treating many human diseases, including cancer, diabetes and obesity. Long-timescale aMD simulations were performed to observe the binding of three chemically diverse ligand molecules: antagonist tiotropium (TTP), partial agonist arecoline (ARc) and full agonist acetylcholine (ACh). In comparison with earlier microsecond-timescale conventional MD simulations, aMD greatly accelerated the binding of ACh to the receptor orthosteric ligand-binding site and the binding of TTP to an extracellular vestibule. Further aMD simulations also captured binding of ARc to the receptor orthosteric site. Additionally, all three ligands were observed to bind in the extracellular vestibule during their binding pathways, suggesting that it is a metastable binding site. This study demonstrates the applicability of aMD to protein–ligand binding, especially the drug recognition of GPCRs.

Type
Report
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/3.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
Copyright © Cambridge University Press 2015

Introduction

Representing the largest family of membrane proteins, G-protein-coupled receptors (GPCRs) mediate cellular responses to hormones and neurotransmitters and the senses of sight olfaction, and taste. They also play a crucial role in the central and parasympathetic nervous systems. Due to their critical involvement in human diseases, GPCRs are primary targets of about one third of current marketed drugs, including treatments for cancer, heart failure, asthma, schizophrenia, Alzheimer's and Parkinson's diseases (Kow & Nathanson, Reference Kow and Nathanson2012; Lappano & Maggiolini, Reference Lappano and Maggiolini2011).

GPCRs exist in an ensemble of different conformations and are known to bind a wide spectrum of ligands. Binding of agonists and inverse agonists in the orthosteric site biases the receptor conformational equilibrium towards the active and inactive states, respectively. GPCRs also bind neutral antagonists, that have no signalling effects but block the receptors from binding other ligands, as well as partial agonists, which induce only submaximal activity (Spalding & Burstein, Reference Spalding and Burstein2006). Additionally, the dynamics and functions of GPCRs can be further regulated by binding of various allosteric modulators, which can impose cell-signalling effects alone or affect the binding affinity and/or signalling efficacy of the orthosteric ligands (Christopoulos, Reference Christopoulos2002; Jeffrey Conn et al. Reference Jeffrey Conn, Christopoulos and Lindsley2009).

It is of paramount importance to understand how drug molecules bind to protein targets such as GPCRs. Detailed characterization of drug-binding pathways to the proteins would provide useful information for effective design of pharmaceutical therapeutics. Using the specialized supercomputer ‘Anton’, microsecond-timescale conventional molecular dynamics simulations captured the processes of a ligand binding to the Src protein kinase (Shan et al. Reference Shan, Kim, Eastwood, Dror, Seeliger and Shaw2011) and more recently to the β 1- and β 2-adrenergic receptors, which are two prototypical GPCRs (Dror et al. Reference Dror, Pan, Arlow, Borhani, Maragakis, Shan, Xu and Shaw2011b). Anton simulations were also applied to the M2 and M3 muscarinic GPCRs (Dror et al. Reference Dror, Green, Valant, Borhani, Valcourt, Pan, Arlow, Canals, Lane, Rahmani, Baell, Sexton, Christopoulos and Shaw2013; Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012). The binding of the endogenous agonist acetylcholine (ACh) to the orthosteric site in the M3 receptor was observed during a 25 μs simulation. Another three Anton simulations (one 16 μs and two 1 μs) captured the binding of antagonist tiotropium (TTP) to the extracellular vestibule. The extensive conventional molecular dynamics (MD) simulations were not able to capture the binding of TTP to the receptor orthosteric site as observed in the X-ray crystal structure, principally due to the fact that TTP is significantly larger than ACh (Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012).

Accelerated molecular dynamics (aMD) is an enhanced sampling technique in which a non-negative boost potential is added to the system's potential energy when it drops below a certain threshold, effectively decreasing the energy barriers and thus accelerating transitions between the low-energy states (Hamelberg et al. Reference Hamelberg, Mongan and McCammon2004, Reference Hamelberg, de Oliveira and McCammon2007; Markwick & McCammon, Reference Markwick and McCammon2011). AMD has been successfully applied to a number of systems (Bucher et al. Reference Bucher, Grant, Markwick and McCammon2011; Gasper et al. Reference Gasper, Fuglestad, Komives, Markwick and McCammon2012; Markwick et al. Reference Markwick, Pierce, Goodin and McCammon2011; Wang et al. Reference Wang, Markwick, de Oliveira and McCammon2011b; Wereszczynski & McCammon, Reference Wereszczynski and McCammon2012) and hundreds-of-nanosecond aMD simulations have been shown to capture millisecond-timescale events (Pierce et al. Reference Pierce, Salomon-Ferrer, de Oliveira, McCammon and Walker2012), including the activation of the M2 and M3 muscarinic receptors (Miao et al. Reference Miao, Nichols, Gasper, Metzger and McCammon2013, Reference Miao, Caliman and McCammon2014a, Reference Miao, Nichols and McCammonb). Based on the funnel-shaped free-energy landscape theory (Frauenfelder et al. Reference Frauenfelder, Sligar and Wolynes1991; Onuchic et al. Reference Onuchic, Luthey-Schulten and Wolynes1997), previous studies suggested that both protein–ligand binding and protein conformational changes (especially folding) involve minimization of the free energy across various energy barriers (Tsai et al. Reference Tsai, Ma and Nussinov1999). It is thus appealing to examine the applicability of aMD to protein–ligand binding.

Here, aMD is used to simulate ligand binding to the M3 muscarinic GPCR, which has been targeted for treating many human diseases, including cancer (Spindel, Reference Spindel2012), diabetes (de Azua et al. Reference de Azua, Scarselli, Rosemond, Gautam, Jou, Gavrilova, Ebert, Levitt and Wess2010; Gregory et al. Reference Gregory, Sexton and Christopoulos2007) and obesity (Weston-Green et al. Reference Weston-Green, Huang, Lian and Deng2012). aMD simulations were performed to observe the binding of three known ligands to the M3 receptor: the antagonist TTP, partial agonist arecoline (ARc) (Kurian et al. Reference Kurian, Hall, Wilkinson, Sullivan, Tobin and Willars2009) and the full agonist ACh (Fig. 1). These simulations elucidate key features of the ligand-binding pathways and highlight metastable binding sites in significantly shorter simulation time than would be required with conventional MD.

Fig. 1. (a) Schematic representation of the M3 muscarinic receptor–ligand-binding simulation system and (b) three known ligands of the M3 receptor that are selected for aMD simulations: antagonist TTP, partial agonist Arc, and full agonist ACh.

Materials and methods

System setup

Simulations of the M3 muscarinic receptor were carried out using the inactive TTP-bound X-ray structure (PDB: 4DAJ) that was determined at 3·40 Å resolution (Haga et al. Reference Haga, Kruse, Asada, Yurugi-Kobayashi, Shiroishi, Zhang, Weis, Okada, Kobilka, Haga and Kobayashi2012). Preparation of the M3 receptor followed the same procedure as previously used and is described briefly here (Miao et al. Reference Miao, Caliman and McCammon2014a). To simulate the ligand binding, TTP was removed from the X-ray structure. The T4 lysozyme that was fused into the protein to replace intracellular loop 3 (ICL3) for crystallizing the receptor was omitted from all simulations, based on previous findings that removal of the bulk of ICL3 does not appear to affect GPCR function and ICL3 is highly flexible (Dror et al. Reference Dror, Arlow, Maragakis, Mildorf, Pan, Xu, Borhani and Shaw2011a). All chain termini were capped with neutral groups (acetyl and methylamide). Two disulphide bonds that were resolved in the crystal structure, i.e., C1403·25–C220ECL2 and C5166·61–C5197·29, were maintained in the simulations. Using the psfgen plugin in visual molecular dynamics (VMD) (Humphrey et al. Reference Humphrey, Dalke and Schulten1996), the Asp1132·50 residue was protonated as in previous microsecond-timescale Anton simulations (Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012). All other protein residues were set to the standard CHARMM protonation states at neutral pH (Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012), including the deprotonated Asp1473·32 residue in the orthosteric site.

The M3 receptor was inserted into a palmitoyl-oleoyl- phosphatidyl-choline (POPC) bilayer with all overlapping lipid molecules removed using the Membrane plugin and solvated in a water box using the Solvate plugin in VMD (Humphrey et al. Reference Humphrey, Dalke and Schulten1996). Four ligand molecules were placed at least 40 Å away from the receptor orthosteric site in the bulk solvent of the starting structures for the antagonist TTP, partial agonist Arc, and full agonist ACh (Fig. 1). The system charges were then neutralized with 18 Cl ions. The simulation systems of the M3 receptor initially measured about 80 × 87 × 97 Å3 with 130 lipid molecules, ~11 200 water molecules and a total of ~55 500 atoms. Periodic boundary conditions were applied to all systems.

MD simulations

MD simulations were performed using NAMD2.9 (Phillips et al. Reference Phillips, Braun, Wang, Gumbart, Tajkhorshid, Villa, Chipot, Skeel, Kale and Schulten2005). The CHARMM27 parameter set with CMAP terms was used for the protein (MacKerell et al. Reference MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, Joseph-McCarthy, Kuchnir, Kuczera, Lau, Mattos, Michnick, Ngo, Nguyen, Prodhom, Reiher, Roux, Schlenkrich, Smith, Stote, Straub, Watanabe, Wiorkiewicz-Kuczera, Yin and Karplus1998; MacKerell et al. Reference MacKerell, Feig and Brooks2004), CHARMM36 for the POPC lipids (Klauda et al. Reference Klauda, Venable, Freites, O'Connor, Tobias, Mondragon-Ramirez, Vorobyov, MacKerell and Pastor2010), and TIP3P model for the water molecules (Jorgensen et al. Reference Jorgensen, Chandrasekhar, Madura, Impey and Klein1983). For the ligand molecules, force field parameters of ACh were retrieved from the CHARMM General Force Field (CGenFF) database (Vanommeslaeghe et al. Reference Vanommeslaeghe and Mackerell2012a, Reference Vanommeslaeghe, Raman and Mackerell2012b). However, CGenFF does not include force field parameters of TTP and ARc, so instead they were computed using the General Automated Atomic Model Parameterization (GAAMP) tool (Huang & Roux, Reference Huang and Roux2013). With ab initio quantum mechanical calculations, GAAMP (Huang & Roux, Reference Huang and Roux2013) provides force field parameters that are compatible with CHARMM as used for protein and lipids in the present study. A cut-off distance of 12 Å was used for the van der Waals and short-range electrostatic interactions and the long-range electrostatic interactions were computed with the particle-mesh Ewald summation method (Essmann et al. Reference Essmann, Perera, Berkowitz, Darden, Lee and Pedersen1995) using a grid point density of 1/Å. A 2 fs integration time-step was used for all MD simulations and a multiple-time-stepping algorithm (Phillips et al. Reference Phillips, Braun, Wang, Gumbart, Tajkhorshid, Villa, Chipot, Skeel, Kale and Schulten2005) was employed with bonded and short-range non-bonded interactions computed every time-step and long-range electrostatic interactions every two time-steps. The SHAKE (Ryckaert et al. Reference Ryckaert, Ciccotti and Berendsen1977) algorithm was applied to all hydrogen-containing bonds.

Simulations of the M3 receptor started with equilibration of the lipid tails. With all other atoms fixed, the lipid tails were energy minimized for 1000 steps using the conjugate gradient algorithm and melted with an NVT run for 0·5 ns at 310 K. The systems were further equilibrated using an NPT run at 1 atm and 310 K for 10 ns with 5 kcal (mol Å2)−1 harmonic position restraints applied to the crystallographically identified atoms in the protein and ligand. The system volume was found to decrease with a flexible unit cell applied and level off during the second half of the 10 ns NPT run, suggesting that water molecules, ions, and lipids were well equilibrated surrounding the protein receptor. Final equilibration of the two systems was performed using an NPT run at 1 atm and 310 K for 0·5 ns with all atoms unrestrained. After these minimization and equilibration procedures, the production MD simulations were performed on the systems for 100 ns at 1 atm pressure and 310 K with a constant ratio constraint applied on the lipid bilayer in the X–Y plane.

aMD simulations

With the aMD implemented in NAMD2.9 (Wang et al. Reference Wang, Harrison, Schulten and McCammon2011a), aMD simulations were performed on the M3 receptor–ligand binding using the ‘dual-boost’ version (Hamelberg et al. Reference Hamelberg, de Oliveira and McCammon2007). Boost potential was applied to both dihedral angles and the total energy across all individual atoms with E dihed = V dihed_avg + 0·3 × V dihed_avg, α dihed = 0·3 × V dihed_avg/5; E total = V total_avg + 0·2 × N atoms and α total = 0·2 × N atoms. Three independent 200 ns aMD runs were performed on the binding of ACh, ARc and TTP ligands by restarting from the final structure of the 100 ns conventional MD simulation with random atomic velocity initializations at 310 K. Two of the three simulations of TTP binding were extended to 1000 and 500 ns, respectively. It is important to note that because the free energy and kinetics of the system are modified with aMD, these simulations do not suggest a definitive time-evolution of the binding pathway. Rather, aMD simulations are used to identify metastable ligand-binding sites.

Simulation analysis

To determine how close ligands bound to the orthosteric site of the M3 muscarinic receptor, the root-mean-square deviations (RMSDs) were calculated for the heavy atoms of each diffusing ligand molecule with reference to the X-ray structure (for TTP) or the top-ranked docking pose in the orthosteric site (for ARc and ACh) after aligning simulation frames using the C α atoms of the receptor transmembrane bundle. Structural clustering was performed using the g_cluster tool in GROMACS with the gromos algorithm based on the RMSD of heavy atoms in each ligand molecule after alignment of the C α atoms in the transmembrane helices in each frame to the starting structure (Daura et al. Reference Daura, Gademann, Jaun, Seebach, van Gunsteren and Mark1999; Pronk et al. Reference Pronk, Pall, Schulz, Larsson, Bjelkmar, Apostolov, Shirts, Smith, Kasson, van der Spoel, Hess and Lindahl2013). A 3 Å RMSD cut-off was chosen because it best captured spatially distinct clusters and allowed the top clusters to be representative of the predominant binding sites explored. The clustering was performed on all simulation frames in which the ligand was found within 5 Å of the receptor. The populations of each cluster are given in Table S1. We examined the three most populated clusters and calculated the RMSDs for each ligand found in these clusters (Table S2). Additionally, the receptor residues that stayed in contact with the ligand in at least 90% of the frames belonging to the cluster were identified (Table S3). For each of these contact residues, we also compared their side-chain dihedral angles (χ 1 and χ 2) to their values in the crystal structure (Table S4) and those in the other clusters (Table S5).

Results

The RMSDs were calculated for the heavy atoms in each diffusing ligand relative to the X-ray structure (for TTP) or the top-ranked docking pose in the orthosteric site (for ARc and ACh) (Methods section) and plotted in Fig. 2b and Fig. S1. TTP was observed to bind to the receptor extracellular vestibule on four separate occasions during the 1000 ns simulation, remaining bound for a total of approximately 540 ns, with a minimum RMSD of 7·86 Å relative to the X-ray structure (Table S2). In the other two simulations, TTP bound briefly to the extracellular surface of the protein several times. ARc bound to the receptor orthosteric site in one of the three 200 ns simulations with 2·20 Å minimum RMSD relative to a top-ranked docking pose. ACh bound to the receptor three times in two of the three 200 ns aMD simulations and dissociated from the receptor twice. In one simulation, ACh bound to the extracellular vestibule briefly for approximately 10 ns, then dissociated and bound again, this time reaching the orthosteric site before dissociating again. In the other simulation, ACh bound to the receptor, reached the orthosteric site, and stayed bound for the remainder of the simulation. Both ACh and ARc were very mobile within the receptor ligand-binding pocket.

Fig. 2. (a) Schematic representation of the X-ray crystal structure of the M3 muscarinic receptor bound to the antagonist, TTP. With TTP removed, this structure was used for aMD simulations of the binding of three known ligand molecules: antagonist TTP, partial agonist ARc and full agonist ACh. (b) RMSDs are plotted for the heavy atoms of each ligand relative to the crystal structure (for TTP) or the top-ranked docking pose (for ARc and ACh) in the orthosteric binding site after aligning all simulation frames using the C α atoms of the receptor transmembrane bundle. Here, data are shown only for ligand molecules that bound to the receptor at some point during the simulations (see Fig. S1 for RMSDs of all ligands). Note that in the TTP RMSD plot the blue and black traces represent two different ligand molecules in one simulation, whereas in the ACh RMSD plot the two curves represent ligand molecules in two different simulations. Only one ARc molecule bound to the receptor, thus there is just a single curve in the ARc RMSD plot. In all cases, the ligand is bound to the receptor at RMSD values less than approximately 20 Å.

Next, we performed RMSD-based structural clustering to identify distinct poses of the three ligands bound to the M3 receptor. For antagonist TTP, the three most populated clusters were all located in the extracellular vestibule, but with different orientations of the ligand (Fig. 3a ). In all three clusters, TTP was in contact with residues Phe221ECL2, Leu225ECL2 and Lys5227·32 (Fig. 3b ). These residues are not conserved between the M2 and M3 receptors, and they also contacted TTP during binding in the extracellular vestibule in previous Anton conventional MD simulations of the M3 receptor (Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012). This indicates that the metastable binding site identified in the present aMD simulations is similar to that in the microsecond-timescale Anton conventional MD simulations. Additionally, reorientations of TTP in the extracellular vestibule were observed in both the aMD and Anton conventional MD simulations. Overall, the aMD simulations of TTP binding are consistent with the previous Anton conventional MD simulations. TTP bound in the extracellular vestibule, but did not reach the orthosteric site (Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012).

Fig. 3. (a) The three most populated TTP-binding clusters are located in the extracellular vestibule and are shown in blue, grey and purple, respectively. (b) Key residues interacting with TTP in cluster A (blue) are shown in sticks and the representative receptor structure observed in the aMD simulations is shown in blue ribbons. A full list of contact residues is given in Table S3.

For full agonist ACh, the three most populated clusters identified from the aMD simulations are shown in Fig. 4a : Cluster A represented the orthosteric binding site, Cluster C in the extracellular vestibule, and Cluster B was located between them. In Cluster A, ACh interacted with two of the residues that form the tyrosine lid (Tyr1483·33 and Tyr5066·51), as well as other residues in the orthosteric site, including Phe2395·47 and Trp5036·48 (Fig. 4e ). In Cluster B, ACh was bound above the orthosteric site and interacted with two residues that form the tyrosine lid, Tyr5066·51 and Tyr5297·39, as well as Asp1473·32 (Fig. 4d ). In Cluster C, ACh formed cation–π interactions with residues Phe1242·60 and Tyr1272·63. This suggests that Cluster C corresponds to Centre 2 of the extracellular vestibule as defined by Dror et al. (Reference Dror, Green, Valant, Borhani, Valcourt, Pan, Arlow, Canals, Lane, Rahmani, Baell, Sexton, Christopoulos and Shaw2013) (Fig. 4c ). The trajectory of ACh diffusing between the clusters (Fig. 4b ) suggests that these three metastable binding states interconvert. These clusters are remarkably similar to those identified in the previously published 25 μs Anton conventional MD simulation (Fig. S2), when the same clustering protocol was applied. The top three clusters identified from the Anton simulation also represent the orthosteric binding site (Fig. S2E), Centre 2 in the extracellular vestibule (Fig. S2C) and a site in between (Fig. S2D). Similarly, ACh diffuses between the three clusters in the Anton simulation (Fig. S2B). Overall, the aMD simulations of ACh binding reproduced the key features of the Anton conventional MD simulations.

Fig. 4. (a) The three most populated Ach-binding clusters are shown in purple, blue and grey, respectively. (b) Trajectory of ACh diffusing between the three clusters during a 200 ns aMD simulation. The time evolution for the other 200 ns aMD simulation in which ACh bound to the receptor is plotted in Fig. S3. Key residues in contact with ACh are shown for the (c) cluster C, (d) cluster B and (e) cluster A. A full list of contact residues is given in Table S3.

Like ACh, partial agonist ARc was also very mobile within the receptor–ligand-binding pocket. The three most populated clusters reflected the ligand mobility with one located in the orthosteric site (Cluster A), one in Centre 2 in the extracellular vestibule (Cluster C), and one in between these two sites (Cluster B; Fig. 5a ). In Cluster A, ARc interacted with all three tyrosine residues that form the tyrosine lid, Tyr1483·33, Tyr5066·51, and Tyr5297·39, as well as Asp1473·32, another key residue in the orthosteric site (Fig. 5e ). In Cluster B, ARc maintained contact with Asp1473·32, Tyr1483·33 and Tyr5297·39 (Fig. 5d ). In Cluster C, ARc formed cation–π interactions with Phe1242·60 and Tyr1272·63, the two residues that define the Centre 2 site in the extracellular vestibule (Fig. 5c ). The simulation trajectory of ARc (Fig. 5b ) indicates interconversion between these three clusters, particularly between Clusters A and B.

Fig. 5. (a) The three most populated Arc-binding clusters are shown in blue, purple and grey, respectively. (b) Trajectory of ARc diffusing between the three clusters during a 200 ns aMD simulation in which binding occurred. Key residues interacting with ARc are shown for (c) cluster C, (d) cluster B and (e) cluster A. A full list of contact residues is given in Table S3.

In addition, we examined the receptor residues that interact with the ligand molecules and undergo significant conformational changes during ligand binding. A complete list of residues with significantly different χ 1 and χ 2 angles compared with the crystal structure is provided in Table S4. The average dihedral angles are generally in agreement with the most probable values found in rotamer libraries (Shapovalov & Dunbrack, Reference Shapovalov and Dunbrack2011). For the TTP clusters, many of these residue conformational changes are likely due to the absence of TTP in the orthosteric site during the aMD simulations. For ACh and ARc, several residues have significantly different χ 1 and χ 2 angles when compared with the crystal structure. In Cluster A, this is largely due to the fact that the ACh and ARc ligands are much smaller than TTP. In Clusters B and C, the ligands are not bound in the orthosteric site. Furthermore, the residues that exhibit significantly different χ 1 and χ 2 angles when the bound ligand moves from one cluster to another are listed in Table S5. The few differences observed in residue side-chain dihedrals between the TTP clusters are subtle because TTP stayed in a very similar position in the extracellular vestibule in the three clusters. For ACh, the differences in residue conformations between the clusters is notable. W5257·35 protrudes into the extracellular vestibule in Cluster A, but flipped up parallel to the transmembrane helices to accommodate ACh binding in Cluster B. Additionally, Y5297·39 flipped up from its position in Cluster A to contact ACh in Cluster C. The differences in side-chain dihedrals between the ARc clusters represent reorientations mainly due to large rearrangements of the highly flexible loops ECL1 and ECL2.

Discussion

It is particularly significant that each ligand bound in the extracellular vestibule in the simulations. This suggests that ligand binding in the extracellular vestibule is a common metastable state during the binding of orthosteric ligands. Notably, this result is consistent with the previous experimental finding that orthosteric ligands can bind to the extracellular vestibule of the M2 receptor (Redka et al. Reference Redka, Pisterzi and Wells2008). The extracellular vestibule has also been confirmed as a binding site of allosteric modulators in long-timescale conventional MD simulations (Dror et al. Reference Dror, Green, Valant, Borhani, Valcourt, Pan, Arlow, Canals, Lane, Rahmani, Baell, Sexton, Christopoulos and Shaw2013) and a recent active X-ray structure (Kruse et al. Reference Kruse, Ring, Manglik, Hu, Hu, Eitel, Hubner, Pardon, Valant, Sexton, Christopoulos, Felder, Gmeiner, Steyaert, Weis, Garcia, Wess and Kobilka2013) of the M2 receptor. This allosteric site could be exploited to develop modulators with high muscarinic receptor subtype selectivity.

In summary, aMD captured the binding of ligand molecules to the M3 receptor in significantly shorter simulation time compared with conventional MD (~80 times speedup in the case of ACh). The identified metastable states of the ligands along the binding pathway are in agreement with previous conventional MD simulations (Dror et al. Reference Dror, Green, Valant, Borhani, Valcourt, Pan, Arlow, Canals, Lane, Rahmani, Baell, Sexton, Christopoulos and Shaw2013; Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012) and experimental findings (Kruse et al. Reference Kruse, Hu, Pan, Arlow, Rosenbaum, Rosemond, Green, Liu, Chae, Dror, Shaw, Weis, Wess and Kobilka2012; Redka et al. Reference Redka, Pisterzi and Wells2008). Key residues in the predominant ligand-binding sites have also been identified, which will be highly useful for designing GPCR mutation studies and engineering small molecules for receptor-selective therapeutics.

More generally, this study demonstrates the applicability of aMD to the study of protein–ligand binding. It is important to note that because the free energy and kinetics of the system are modified with aMD, the above simulations do not suggest a definitive time-evolution of the ligand-binding pathway. Rather, aMD simulations are useful for the discovery of metastable ligand-binding sites and to aid the development of effective drugs. In comparison with other methods, aMD has the advantage of significantly shortening the simulation time needed to observe ligand binding without the need to pre-define reaction coordinates as in metadynamics (Laio & Parrinello, Reference Laio and Parrinello2002) and adaptive biasing force calculations (Darve & Pohorille, Reference Darve and Pohorille2001; Rodriguez-Gomez et al. Reference Rodriguez-Gomez, Darve and Pohorille2004). Thus, aMD should be of wide applicability to protein–ligand-binding studies.

Supplementary material

The supplementary material for this article can be found at http://dx.doi.org/10.1017/S0033583515000153

Acknowledgements

We thank Lei Huang at the University of Chicago for assistance with calculating the ligand force field parameters, Yibing Shan for valuable discussions, and Ron Dror, Jodi Hezky and Albert Pan from DE Shaw's Research Group for generously providing the Anton simulation trajectories.

Financial Support

This work was supported by NSF (grant no. MCB1020765), NIH (grant no. GM31749), Howard Hughes Medical Institute, Center for Theoretical Biological Physics (CTBP) and National Biomedical Computation Resource (NBCR). Computing time was provided on the Gordon and Stampede supercomputers through the Extreme Science and Engineering Discovery Environment (XSEDE) award TG-MCB140011 and the Hopper and Edison supercomputers through the National Energy Research Scientific Computing Center (NERSC) project m1395.

Conflict of Interest

None.

References

Bucher, D., Grant, B. J., Markwick, P. R. & McCammon, J. A. (2011). Accessing a hidden conformation of the maltose binding protein using accelerated molecular dynamics. PLoS Computational Biology 7, e1002034.CrossRefGoogle ScholarPubMed
Christopoulos, A. (2002). Allosteric binding sites on cell-surface receptors: novel targets for drug discovery. Nature Reviews Drug Discovery 1, 198210.CrossRefGoogle ScholarPubMed
Darve, E. & Pohorille, A. (2001). Calculating free energies using average force. Journal of Chemical Physics 115, 91699183.CrossRefGoogle Scholar
Daura, X., Gademann, K., Jaun, B., Seebach, D., van Gunsteren, W. F. & Mark, A. E. (1999). Peptide folding: when simulation meets experiment. Angewandte Chemie International Edition 38, 236240.3.0.CO;2-M>CrossRefGoogle Scholar
de Azua, I. R., Scarselli, M., Rosemond, E., Gautam, D., Jou, W., Gavrilova, O., Ebert, P. J., Levitt, P. & Wess, J. (2010). RGS4 is a negative regulator of insulin release from pancreatic beta-cells in vitro and in vivo . Proceedings of the National Academy of Sciences of the United States of America 107, 79998004.CrossRefGoogle Scholar
Dror, R. O., Arlow, D. H., Maragakis, P., Mildorf, T. J., Pan, A. C., Xu, H., Borhani, D. W. & Shaw, D. E. (2011a). Activation mechanism of the β2-adrenergic receptor. Proceedings of the National Academy of Sciences of the United States of America 108, 1868418689.CrossRefGoogle ScholarPubMed
Dror, R. O., Green, H. F., Valant, C., Borhani, D. W., Valcourt, J. R., Pan, A. C., Arlow, D. H., Canals, M., Lane, J. R., Rahmani, R., Baell, J. B., Sexton, P. M., Christopoulos, A. & Shaw, D. E. (2013). Structural basis for modulation of a G-protein-coupled receptor by allosteric drugs. Nature 503, 295299.Google Scholar
Dror, R. O., Pan, A. C., Arlow, D. H., Borhani, D. W., Maragakis, P., Shan, Y., Xu, H. & Shaw, D. E. (2011b). Pathway and mechanism of drug binding to G-protein-coupled receptors. Proceedings of the National Academy of Sciences of the United States of America 108, 1311813123.CrossRefGoogle Scholar
Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H. & Pedersen, L. G. (1995). A smooth particle mesh Ewald method. Journal of Chemical Physics 103, 85778593.CrossRefGoogle Scholar
Frauenfelder, H., Sligar, S. G. & Wolynes, P. G. (1991). The energy landscapes and motions of proteins. Science 254, 15981603.CrossRefGoogle ScholarPubMed
Gasper, P. M., Fuglestad, B., Komives, E. A., Markwick, P. R. L. & McCammon, J. A. (2012). Allosteric networks in thrombin distinguish procoagulant vs. anticoagulant activities. Proceedings of the National Academy of Sciences of the United States of America 109, 2121621222.CrossRefGoogle Scholar
Gregory, K. J., Sexton, P. M. & Christopoulos, A. (2007). Allosteric modulation of muscarinic acetylcholine receptors. Current Neuropharmacology 5, 157167.CrossRefGoogle ScholarPubMed
Haga, K., Kruse, A. C., Asada, H., Yurugi-Kobayashi, T., Shiroishi, M., Zhang, C., Weis, W. I., Okada, T., Kobilka, B. K., Haga, T. & Kobayashi, T. (2012). Structure of the human M2 muscarinic acetylcholine receptor bound to an antagonist. Nature 482, 547551.CrossRefGoogle Scholar
Hamelberg, D., de Oliveira, C. A. F. & McCammon, J. A. (2007). Sampling of slow diffusive conformational transitions with accelerated molecular dynamics. Journal of Chemical Physics 127, 155102.CrossRefGoogle Scholar
Hamelberg, D., Mongan, J. & McCammon, J. A. (2004). Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules. Journal of Chemical Physics 120, 1191911929.CrossRefGoogle ScholarPubMed
Huang, L. & Roux, B. (2013). Automated force field parameterization for nonpolarizable and polarizable atomic models based on ab initio target data. Journal of Chemical Theory and Computation 9, 35433556.CrossRefGoogle Scholar
Humphrey, W., Dalke, A. & Schulten, K. (1996). VMD: visual molecular dynamics. Journal of Molecular Graphics & Modelling 14, 3338.CrossRefGoogle ScholarPubMed
Jeffrey Conn, P., Christopoulos, A. & Lindsley, C. W. (2009). Allosteric modulators of GPCRs: a novel approach for the treatment of CNS disorders. Nature Reviews Drug Discovery 8, 4154.CrossRefGoogle Scholar
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. Journal of Chemical Physics 79, 926935.CrossRefGoogle Scholar
Klauda, J. B., Venable, R. M., Freites, J. A., O'Connor, J. W., Tobias, D. J., Mondragon-Ramirez, C., Vorobyov, I., MacKerell, A. D. & Pastor, R. W. (2010). Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. The Journal of Physical Chemistry B 114, 78307843.CrossRefGoogle ScholarPubMed
Kow, R. L. & Nathanson, N. M. (2012). Structural biology: muscarinic receptors become crystal clear. Nature 482, 480481.CrossRefGoogle ScholarPubMed
Kruse, A. C., Hu, J., Pan, A. C., Arlow, D. H., Rosenbaum, D. M., Rosemond, E., Green, H. F., Liu, T., Chae, P. S., Dror, R. O., Shaw, D. E., Weis, W. I., Wess, J. & Kobilka, B. K. (2012). Structure and dynamics of the M3 muscarinic acetylcholine receptor. Nature 482, 552556.CrossRefGoogle ScholarPubMed
Kruse, A. C., Ring, A. M., Manglik, A., Hu, J., Hu, K., Eitel, K., Hubner, H., Pardon, E., Valant, C., Sexton, P. M., Christopoulos, A., Felder, C. C., Gmeiner, P., Steyaert, J., Weis, W. I., Garcia, K. C., Wess, J. & Kobilka, B. K. (2013). Activation and allosteric modulation of a muscarinic acetylcholine receptor. Nature 504, 101106.CrossRefGoogle Scholar
Kurian, N., Hall, C. J., Wilkinson, G. F., Sullivan, M., Tobin, A. B. & Willars, G. B. (2009). Full and partial agonists of muscarinic M3 receptors reveal single and oscillatory Ca2+ responses by beta 2-adrenoceptors. Journal of Pharmacology and Experimental Therapeutics 330, 502512.CrossRefGoogle Scholar
Laio, A. & Parrinello, M. (2002). Escaping free-energy minima. Proceedings of the National Academy of Sciences of the United States of America 99, 1256212566.CrossRefGoogle ScholarPubMed
Lappano, R. & Maggiolini, M. (2011). G protein-coupled receptors: novel targets for drug discovery in cancer. Nature Reviews Drug Discovery 10, 4760.CrossRefGoogle Scholar
MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. J., Fischer, S., Gao, J., Guo, H., Ha, S., Joseph-McCarthy, D., Kuchnir, L., Kuczera, K., Lau, F. T. K., Mattos, C., Michnick, S., Ngo, T., Nguyen, D. T., Prodhom, B., Reiher, W. E., Roux, B., Schlenkrich, M., Smith, J. C., Stote, R., Straub, J., Watanabe, M., Wiorkiewicz-Kuczera, J., Yin, D. & Karplus, M. (1998). All-atom empirical potential for molecular modeling and dynamics studies of proteins. Journal of Physical Chemistry B 102, 35863616.CrossRefGoogle ScholarPubMed
MacKerell, A. D. Jr., Feig, M. & Brooks, C. L. III (2004). Improved treatment of the protein backbone in empirical force fields. Journal of the American Chemical Society 126, 698699.CrossRefGoogle ScholarPubMed
Markwick, P. R. L. & McCammon, J. A. (2011). Studying functional dynamics in bio-molecules using accelerated molecular dynamics. Physical Chemistry Chemical Physics 13, 2005320065.CrossRefGoogle ScholarPubMed
Markwick, P. R. L., Pierce, L. C. T., Goodin, D. B. & McCammon, J. A. (2011). Adaptive accelerated molecular dynamics (Ad-AMD) revealing the molecular plasticity of P450cam. Journal of Physical Chemistry Letters 2, 158164.CrossRefGoogle ScholarPubMed
Miao, Y., Caliman, A. D. & McCammon, J. A. (2014a). Allosteric effects of sodium ion binding on activation of the M3 muscarinic G-protein coupled receptor. Biophysical Journal 108, 17961806.CrossRefGoogle ScholarPubMed
Miao, Y., Nichols, S. E., Gasper, P. M., Metzger, V. T. & McCammon, J. A. (2013). Activation and dynamic network of the M2 muscarinic receptor. Proceedings of the National Academy of Sciences of the United States of America 110, 1098210987.CrossRefGoogle ScholarPubMed
Miao, Y., Nichols, S. E. & McCammon, J. A. (2014b). Free energy landscape of G-protein coupled receptors, explored by accelerated molecular dynamics. Physical Chemistry Chemical Physics 16, 63986406.CrossRefGoogle ScholarPubMed
Onuchic, J. N., Luthey-Schulten, Z. & Wolynes, P. G. (1997). Theory of protein folding: the energy landscape perspective. Annual Review of Physical Chemistry 48, 545600.CrossRefGoogle ScholarPubMed
Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., Chipot, C., Skeel, R. D., Kale, L. & Schulten, K. (2005). Scalable molecular dynamics with NAMD. Journal of Computational Chemistry 26, 17811802.CrossRefGoogle ScholarPubMed
Pierce, L. C. T., Salomon-Ferrer, R., de Oliveira, C. A. F., McCammon, J. A. & Walker, R. C. (2012). Routine access to millisecond time scale events with accelerated molecular dynamics. Journal of Chemical Theory and Computation 8, 29973002.CrossRefGoogle ScholarPubMed
Pronk, S., Pall, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., Shirts, M. R., Smith, J. C., Kasson, P. M., van der Spoel, D., Hess, B. & Lindahl, E. (2013). GROMACS 4·5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 29, 845854.CrossRefGoogle Scholar
Redka, D. S., Pisterzi, L. F. & Wells, J. W. (2008). Binding of orthosteric ligands to the allosteric site of the M(2) muscarinic cholinergic receptor. Molecular Pharmacology 74, 834843.CrossRefGoogle Scholar
Rodriguez-Gomez, D., Darve, E. & Pohorille, A. (2004). Assessing the efficiency of free energy calculation methods. Journal of Chemical Physics 120, 35633578.CrossRefGoogle ScholarPubMed
Ryckaert, J.-P., Ciccotti, G. & Berendsen, H. J. C. (1977). Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. Journal of Computational Physics 23, 327341.CrossRefGoogle Scholar
Shan, Y., Kim, E. T., Eastwood, M. P., Dror, R. O., Seeliger, M. A. & Shaw, D. E. (2011). How does a drug molecule find its target binding site? Journal of the American Chemical Society 133, 91819183.CrossRefGoogle ScholarPubMed
Shapovalov, M. V. & Dunbrack, R. L. (2011). A smoothed backbone-dependent rotamer library for proteins derived from adaptive kernel density estimates and regressions. Structure 19, 844858.CrossRefGoogle ScholarPubMed
Spalding, T. A. & Burstein, E. S. (2006). Constitutive activity of muscarinic acetylcholine receptors. Journal of Receptors and Signal Transduction 26, 6185.CrossRefGoogle ScholarPubMed
Spindel, E. R. (2012). Muscarinic receptor agonists and antagonists: effects on cancer. Handbook of Experimental Pharmacology 208, 451468.CrossRefGoogle Scholar
Tsai, C. J., Ma, B. Y. & Nussinov, R. (1999). Folding and binding cascades: shifts in energy landscapes. Proceedings of the National Academy of Sciences of the United States of America 96, 99709972.CrossRefGoogle ScholarPubMed
Vanommeslaeghe, K. & Mackerell, A. D. (2012a). Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing. Journal of Chemical Information and Modeling 52, 31443154.CrossRefGoogle Scholar
Vanommeslaeghe, K., Raman, E. P. & Mackerell, A. D. (2012b). Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. Journal of Chemical Information and Modeling 52, 31553168.CrossRefGoogle ScholarPubMed
Wang, Y., Harrison, C. B., Schulten, K. & McCammon, J. A. (2011a). Implementation of accelerated molecular dynamics in NAMD. Computational Science and Discovery 4, 015002.CrossRefGoogle ScholarPubMed
Wang, Y., Markwick, P. R. L., de Oliveira, C. A. F. & McCammon, J. A. (2011b). Enhanced lipid diffusion and mixing in accelerated molecular dynamics. Journal of Chemical Theory and Computation 7, 31993207.CrossRefGoogle ScholarPubMed
Wereszczynski, J. & McCammon, J. A. (2012). Nucleotide- dependent mechanism of Get3 as elucidated from free energy calculations. Proceedings of the National Academy of Sciences of the United States of America 109, 77597764.CrossRefGoogle ScholarPubMed
Weston-Green, K., Huang, X. F., Lian, J. M. & Deng, C. (2012). Effects of olanzapine on muscarinic M3 receptor binding density in the brain relates to weight gain, plasma insulin and metabolic hormone levels. European Neuropsychopharmacology 22, 364373.CrossRefGoogle ScholarPubMed
Figure 0

Fig. 1. (a) Schematic representation of the M3 muscarinic receptor–ligand-binding simulation system and (b) three known ligands of the M3 receptor that are selected for aMD simulations: antagonist TTP, partial agonist Arc, and full agonist ACh.

Figure 1

Fig. 2. (a) Schematic representation of the X-ray crystal structure of the M3 muscarinic receptor bound to the antagonist, TTP. With TTP removed, this structure was used for aMD simulations of the binding of three known ligand molecules: antagonist TTP, partial agonist ARc and full agonist ACh. (b) RMSDs are plotted for the heavy atoms of each ligand relative to the crystal structure (for TTP) or the top-ranked docking pose (for ARc and ACh) in the orthosteric binding site after aligning all simulation frames using the Cα atoms of the receptor transmembrane bundle. Here, data are shown only for ligand molecules that bound to the receptor at some point during the simulations (see Fig. S1 for RMSDs of all ligands). Note that in the TTP RMSD plot the blue and black traces represent two different ligand molecules in one simulation, whereas in the ACh RMSD plot the two curves represent ligand molecules in two different simulations. Only one ARc molecule bound to the receptor, thus there is just a single curve in the ARc RMSD plot. In all cases, the ligand is bound to the receptor at RMSD values less than approximately 20 Å.

Figure 2

Fig. 3. (a) The three most populated TTP-binding clusters are located in the extracellular vestibule and are shown in blue, grey and purple, respectively. (b) Key residues interacting with TTP in cluster A (blue) are shown in sticks and the representative receptor structure observed in the aMD simulations is shown in blue ribbons. A full list of contact residues is given in Table S3.

Figure 3

Fig. 4. (a) The three most populated Ach-binding clusters are shown in purple, blue and grey, respectively. (b) Trajectory of ACh diffusing between the three clusters during a 200 ns aMD simulation. The time evolution for the other 200 ns aMD simulation in which ACh bound to the receptor is plotted in Fig. S3. Key residues in contact with ACh are shown for the (c) cluster C, (d) cluster B and (e) cluster A. A full list of contact residues is given in Table S3.

Figure 4

Fig. 5. (a) The three most populated Arc-binding clusters are shown in blue, purple and grey, respectively. (b) Trajectory of ARc diffusing between the three clusters during a 200 ns aMD simulation in which binding occurred. Key residues interacting with ARc are shown for (c) cluster C, (d) cluster B and (e) cluster A. A full list of contact residues is given in Table S3.

Supplementary material: File

Kappel supplementary material

Figures S1-S3 and Tables S1-S5

Download Kappel supplementary material(File)
File 2 MB
You have Access Open access
99
Cited by

Save article to Kindle

To save this article to your Kindle, first ensure coreplatform@cambridge.org is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about saving to your Kindle.

Note you can select to save to either the @free.kindle.com or @kindle.com variations. ‘@free.kindle.com’ emails are free but can only be saved to your device when it is connected to wi-fi. ‘@kindle.com’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

Find out more about the Kindle Personal Document Service.

Accelerated molecular dynamics simulations of ligand binding to a muscarinic G-protein-coupled receptor
Available formats
×

Save article to Dropbox

To save this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you used this feature, you will be asked to authorise Cambridge Core to connect with your Dropbox account. Find out more about saving content to Dropbox.

Accelerated molecular dynamics simulations of ligand binding to a muscarinic G-protein-coupled receptor
Available formats
×

Save article to Google Drive

To save this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you used this feature, you will be asked to authorise Cambridge Core to connect with your Google Drive account. Find out more about saving content to Google Drive.

Accelerated molecular dynamics simulations of ligand binding to a muscarinic G-protein-coupled receptor
Available formats
×
×

Reply to: Submit a response

Please enter your response.

Your details

Please enter a valid email address.

Conflicting interests

Do you have any conflicting interests? *