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

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.

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, 2012; Lappano & Maggiolini, 2011).

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, 2006). 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, 2002; Jeffrey Conn et al. 2009).

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. 2011) and more recently to the β ₁- and β ₂-adrenergic receptors, which are two prototypical GPCRs (Dror et al. 2011b). Anton simulations were also applied to the M2 and M3 muscarinic GPCRs (Dror et al. 2013; Kruse et al. 2012). 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. 2012).

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. 2004, 2007; Markwick & McCammon, 2011). AMD has been successfully applied to a number of systems (Bucher et al. 2011; Gasper et al. 2012; Markwick et al. 2011; Wang et al. 2011b; Wereszczynski & McCammon, 2012) and hundreds-of-nanosecond aMD simulations have been shown to capture millisecond-timescale events (Pierce et al. 2012), including the activation of the M2 and M3 muscarinic receptors (Miao et al. 2013, 2014a, b). Based on the funnel-shaped free-energy landscape theory (Frauenfelder et al. 1991; Onuchic et al. 1997), 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. 1999). 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, 2012), diabetes (de Azua et al. 2010; Gregory et al. 2007) and obesity (Weston-Green et al. 2012). 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. 2009) 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. 2012). Preparation of the M3 receptor followed the same procedure as previously used and is described briefly here (Miao et al. 2014a). 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. 2011a). All chain termini were capped with neutral groups (acetyl and methylamide). Two disulphide bonds that were resolved in the crystal structure, i.e., C140^3·25–C220^ECL2 and C516^6·61–C519^7·29, were maintained in the simulations. Using the psfgen plugin in visual molecular dynamics (VMD) (Humphrey et al. 1996), the Asp113^2·50 residue was protonated as in previous microsecond-timescale Anton simulations (Kruse et al. 2012). All other protein residues were set to the standard CHARMM protonation states at neutral pH (Kruse et al. 2012), including the deprotonated Asp147^3·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. 1996). 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 Å³ 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. 2005). The CHARMM27 parameter set with CMAP terms was used for the protein (MacKerell et al. 1998; MacKerell et al. 2004), CHARMM36 for the POPC lipids (Klauda et al. 2010), and TIP3P model for the water molecules (Jorgensen et al. 1983). For the ligand molecules, force field parameters of ACh were retrieved from the CHARMM General Force Field (CGenFF) database (Vanommeslaeghe et al. 2012a, 2012b). 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, 2013). With ab initio quantum mechanical calculations, GAAMP (Huang & Roux, 2013) 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. 1995) 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. 2005) 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. 1977) 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 Ų)⁻¹ 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. 2011a), aMD simulations were performed on the M3 receptor–ligand binding using the ‘dual-boost’ version (Hamelberg et al. 2007). 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. 1999; Pronk et al. 2013). 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 (χ ₁ and χ ₂) 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 Phe221^ECL2, Leu225^ECL2 and Lys522^7·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. 2012). 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. 2012).

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 (Tyr148^3·33 and Tyr506^6·51), as well as other residues in the orthosteric site, including Phe239^5·47 and Trp503^6·48 (Fig. 4e ). In Cluster B, ACh was bound above the orthosteric site and interacted with two residues that form the tyrosine lid, Tyr506^6·51 and Tyr529^7·39, as well as Asp147^3·32 (Fig. 4d ). In Cluster C, ACh formed cation–π interactions with residues Phe124^2·60 and Tyr127^2·63. This suggests that Cluster C corresponds to Centre 2 of the extracellular vestibule as defined by Dror et al. (2013) (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, Tyr148^3·33, Tyr506^6·51, and Tyr529^7·39, as well as Asp147^3·32, another key residue in the orthosteric site (Fig. 5e ). In Cluster B, ARc maintained contact with Asp147^3·32, Tyr148^3·33 and Tyr529^7·39 (Fig. 5d ). In Cluster C, ARc formed cation–π interactions with Phe124^2·60 and Tyr127^2·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 χ ₁ and χ ₂ 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, 2011). 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 χ ₁ and χ ₂ 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 χ ₁ and χ ₂ 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. W525^7·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, Y529^7·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. 2008). The extracellular vestibule has also been confirmed as a binding site of allosteric modulators in long-timescale conventional MD simulations (Dror et al. 2013) and a recent active X-ray structure (Kruse et al. 2013) 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. 2013; Kruse et al. 2012) and experimental findings (Kruse et al. 2012; Redka et al. 2008). 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, 2002) and adaptive biasing force calculations (Darve & Pohorille, 2001; Rodriguez-Gomez et al. 2004). Thus, aMD should be of wide applicability to protein–ligand-binding studies.