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.
Open access
