Insights from Molecular Dynamics Simulations - American Chemical

Jan 21, 2015 - stress, anxiety, depression, and irritable bowel syndrome.38 The two available crystal structures of this class (PDB IDs 4K5Y and. 4L6R...
0 downloads 0 Views 6MB Size
Article pubs.acs.org/JPCB

Activation of Corticotropin-Releasing Factor 1 Receptor: Insights from Molecular Dynamics Simulations Rajesh Singh, Navjeet Ahalawat, and Rajesh K. Murarka* Department of Chemistry, Indian Institute of Science Education and Research Bhopal, Indore By-pass Road, Bhauri, Bhopal 462066, MP, India S Supporting Information *

ABSTRACT: G-protein-coupled receptors (GPCRs) constitute the largest family of membrane-bound proteins involved in translation of extracellular signals into intracellular responses. They regulate diverse physiological and pathophysiological processes, and hence, they are prime drug targets for therapeutic intervention. In spite of the recent advancements in membrane protein crystallography, limited information is available on the molecular signatures of activation of GPCRs. Although few studies have been reported for class A GPCRs, the activation mechanism of class B GPCRs remains unexplored. Corticotropin-releasing factor 1 receptor (CRF1R), a class B GPCR, is associated with various disease conditions including stress, anxiety, and irritable bowel syndrome. Here, we report the activation of CRF1R using accelerated molecular dynamics simulations of the apo receptor. The breakage of His1552.50−Glu2093.50 and Glu2093.50−Thr3166.42 interactions is found to be crucial in transition of the receptor to its active conformation. Compared to the inactive crystal structure, major structural rearrangements occurred in the intracellular region of the transmembrane (TM) domain upon activation: TM3 twisted away from TM2, and an opening of the G-protein binding site occurred as a result of the outward movements of TM5 and TM6 from the helical bundle. Further, an inward tilt of TM7 toward the helical core is observed at the extracellular side, in agreement with recent findings (Coin et al. Cell 2013, 155, 1258−1269), where it is proposed that this movement helps in establishing favorable interactions with peptide agonist. Moreover, different allosteric pathways in the inactive and active states are identified using the correlations in torsion angle space. The inactive state is found to be less dynamic as compared to the putative active state of the receptor. Results from the current study could present a model for class B GPCRs activation and aid in the design of CRF1R modulators against brain and metabolic disorders.



INTRODUCTION One of the important superfamily of transmembrane proteins from eukaryotes consists of the G-protein-coupled receptors (GPCRs). These receptors are responsible for transducing the extracellular signal into a physiological response and are involved in vital physiological processes that includes cell growth, differentiation, metabolism, sense, neurotransmission, secretion, and defense.1 Due to their essential physiological roles and occurrence at the cell surface, these receptors act as excellent drug targets. More than 30% of the marketed drugs modulate the activity of these receptors either directly or indirectly.2,3 GPCRs exist in multiple conformational states4−6 and show a considerable amount of basal activity7,8 by activating the G-proteins even in the absence of an agonist. The advent of high-resolution crystal structures of a number of © XXXX American Chemical Society

GPCRs bound to variety of ligands provided a wealth of information for the design of small-molecule modulators of this important superfamily of receptors.9−11 However, the crystal structures only deduce the static pictures of these highly dynamic proteins and are unable to elucidate the detailed mechanism of their activation.11−15 A thorough understanding of the activation process at the atomic level is absolutely essential to rationalize the structure based drug design process. Computational simulations and/or state-of-the-art experimental techniques are the current rational approaches that can provide insights into the structural dynamics of these receptors.13,16−23 Received: September 29, 2014 Revised: January 13, 2015

A

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article



METHODS System Preparation. X-ray crystal structure of CRF1R37 in complex with the antagonist CP-376395 (PDB ID: 4K5Y, resolution: 2.98 Å) and without antagonist (apo) were taken for this study after removing the T4 lysozyme existing in the crystal structure. Missing residues in the loops IL2 and EL3 (i.e., 221 to 223 and 334 to 337, respectively) were modeled using Modeller9v10.43 The disulfide bond between Cys188− Cys258, present in the crystal structure, was retained. Chain termini were capped with neutral groups (acetyl and methylamide). CRF1R residues were set to the standard CHARMM protonation states at neutral pH. The receptor was inserted into a palmitoyl−oleoyl−phosphatidyl−choline (POPC) bilayer with all overlapping lipid molecules removed. The system was then neutralized using 0.15 M NaCl. The resulting system for the apo receptor consist of 178 POPC lipid molecules, 22 Na+, 30 Cl− and 4536 TIP3P44 water molecules, with a total of 50 663 atoms, and measured as 78 × 78 × 89 Å3. Molecular Dynamics Simulations. All-atom molecular dynamics (MD) simulations of the CRF1 receptor were carried out using NAMD2.9 for both normal MD (nMD) and accelerated MD (aMD) simulations.45,46 The CHARMM27 force-field parameter set was used for the receptor (with CMAP terms included),47,48 CHARMM36 force-field for the POPC lipids,49 and the CHARMM general force field parameter set was used for CP-376395.50 Periodic boundary conditions were applied, and the cutoffs for the short-range electrostatic and van der Waals interactions were set to 12 Å, with a switching function applied at 10 Å. Long-range electrostatic interaction was calculated using the particle-mesh Ewald method.51 The SHAKE algorithm52 was applied to constrain all the hydrogencontaining bonds, and a time step of 2 fs was used. For performing nMD simulations of the apo and CP-376395bound structures, the two systems were equilibrated at constant pressure (1 atm) and constant temperature (310 K) for 10 ns, and for each system, five such independent simulation runs were performed starting from different initial atomic velocities. After minimization and equilibration procedures, the production runs of the nMD simulations were performed on the two systems for 150 ns at 1 atm pressure and 310 K with a constant ratio constraint applied on the lipid bilayer in the X−Y plane. To enhance the sampling of the conformational space, the accelerated molecular dynamics technique46 was applied using the dual-boost procedure in NAMD2.9.45 Dual-boost aMD was applied to the final structures acquired from the nMD simulations of both the apo and the antagonist-bound structures. A dihedral and a total boost potentials (Edih, αdih, Etotal, αtotal) were added to all the atoms in the system.

Class A (rhodopsin-like) constitute the largest and most diverse class of GPCRs. Extensive experimental and computational studies carried out for this class of receptors provided atomistic details of their structures,2,15,24,25 (including in complex with G-proteins15,26), allosteric modulation by ligands,17,27 and activation.14,19,28,29 The second largest family of GPCRs is class B (secretin receptor family), which is an important class for targeting diseases including diabetes, osteoporosis, cardiovascular, anxiety, depression, and other central nervous system disorders.30 In spite of much effort by the pharmaceutical industry, only a few compounds were able to successfully modulate the action of endogenous peptide agonists of class B receptors.31−33 Most of the compounds under investigation have poor pharmacokinetic profiles due to the limited understanding of the druggable site. Structural information for class B was previously limited to the structures of the N-terminal extracellular domains bound to the peptide ligands.34,35 Crystal structures of the transmembrane domains of human glucagon receptor (GCGR) and human corticotropin-releasing factor 1 receptor (CRF1R) have only been solved very recently.36,37 GCGR, a potential drug target for diabetes, binds to glucagon peptide and plays an important role in glucose homeostasis. CRF1R binds to corticotropin-releasing hormone or urocortin peptide UCN1 and is an important target for the diseases including chronic stress, anxiety, depression, and irritable bowel syndrome.38 The two available crystal structures of this class (PDB IDs 4K5Y and 4L6R) show similar overall topology.30 GCGR was crystallized in the presence of an antagonist; however, due to insufficient electron density, the small-molecule-binding region was not ascertained. The crystal structure of CRF1R revealed an antagonist-binding site that reside deep into the TM domain, ∼18 Å away from the small-molecule binding sites, as reported for class A GPCRs.39 Residues from transmembranes TM3, TM5, and TM6 are involved in lining the pocket characterized for the binding of antagonist CP-376395. The crystal structures of both GCGR and CRF1R represent inactive states, and currently, no active-state structure of class B receptor is available. The wealth of information available for class A GPCRs revealed the presence of molecular switches/conserved motifs responsible for receptor activation through a series of conformational changes. These include the conserved DRY motif that forms an “ionic lock” (between TM3 and TM6), tryptophan rotamer toggle switch (Trp6.48: Ballesteros− Weinstein numbering) and the conserved NPXXY motif.40 On the other hand, class B receptors lack these conserved motifs, and the mechanism by which these receptors transition between inactive and active conformations remains unclear. In the present study, we aim to identify the active state of CRF1R with the help of accelerated molecular dynamics (aMD), capable of sampling millisecond time scale events in nanoseconds of simulations. This technique modifies the complex potential energy landscape by raising energy minima that are below a certain threshold level while the areas above this threshold are unaltered.41,42 Comparison of simulation results obtained for antagonist-bound and unliganded (apo) states of CRF1R allowed us to characterize the structural rearrangements in the transmembranes while transitioning to its active state. The activation mechanism proposed here provides an insight into the functioning of this class B GPCR.

V *(r ) = V (r ),

V (r ) ≥ E

V *(r ) = V (r ) + ΔV(r ),

V (r ) < E

where V(r) is the original potential, E is the reference energy, and V*(r) is the modified potential. The boost potential, ΔV(r), is expressed as ΔV (r ) =

(E − V (r ))2 α + E − V (r )

where α is the acceleration factor. Decrease in α results into a more flattened potential energy surface and the transitions between the low-energy states are accelerated due to low B

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the 12 residues toward the respective side. Active and inactive clusters are identified taking the low energy conformations from the regions defined by the order parameters, cytoplasmic TM2−TM6 distance (inactive = 15−17 Å; active = 21.5−23.5 Å) and His155−Glu209 distance (inactive = 3.5−6.0 Å, active = 10−13 Å), respectively (Figure 1).

energy barriers. The dihedral and total boost acceleration parameters are obtained as Edih = Vdih_avg + λ*Vdih_avg

αdih = λ*Vdih_avg /5 Etotal = Vtotal_avg + 0.2*Natoms αtotal = 0.2*Natoms

where Natoms is the total number of atoms, Vdih_avg and Vtotal_avg are the average dihedral and total potential energies calculated from 150 ns nMD simulations, respectively, and λ is an adjustable parameter for acceleration. Applying the boost potential results in enhanced sampling of the conformational space and thus increases the probability of observing the active state. It should be noted that the application of a large boost to the potential could result into a considerable loss (>5%) of the secondary structure elements, and hence, it is not acceptable. We have applied an optimum boost to the system while retaining the acceptable secondary structure elements. We performed three test aMD simulations (50 ns each) in order to choose optimal acceleration with λ values of 0.2, 0.3, and 0.4. We achieved the proper acceleration at λ = 0.3 with MIavg) and were >10 Å apart in the receptor structure. The top 500 allosteric pathways were used to show pipelines (Table S1, SI). For more detail about clustering of the pathways to generate allosteric pipelines and calculation of MI, see the SI. Allosteric network in the CRF1R is examined through dynamical network analysis using the NetworkView56,57 plugin in VMD.58 In dynamical network, each Cα atom of residues represents a node and two residues are connected by an edge if in any two of their heavy atoms are within 5 Å for more than 75% of the simulation time. Edge weight, wij, that connects node i and j is calculated by the formula: wij = −log(rij), where rij represents the corresponding normalized mutual information (between two residues) derived from torsional degrees of freedom. The Girvan−Newman algorithm is used to divide the network into communities. The communities are groups of nodes within which the nodes are densely intraconnected, and sparsely interconnected.57 The number of shortest paths that cross a given edge is calculated for all edges and the edge with greatest betweenness is removed. This process is repeated and a modularity score is tracked to identify the division that results in an optimal community network. Critical nodes located at the interface of neighboring communities and the connecting edges that describe the probability of information transfer (betweenness) are then obtained.



RESULTS Schematic representation of the system used in simulations, CRF1R inserted into a POPC bilayer and solvated in an aqueous medium neutralized by 0.15 M NaCl, is shown in Figure S1, SI. Five normal molecular dynamics (nMD) simulations of 150 ns each (cumulative total of 750 ns; Table 1) of the inactive X-ray structure of CRF1 receptor (PDB ID: 4K5Y) bound to the antagonist CP-37639537 do not show significant changes as compared to the crystal structure. Higher fluctuations were observed in the loop regions as compared to the seven transmembrane core of the receptor. Capturing GPCR activation using nMD simulations is a very challenging task,59 though the reverse process (i.e., active to inactive transition) has been studied in the case of β2-adrenergic receptor.16 In the absence of an active class B GPCR structure, D

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. Interactions of (A) His1552.50 and Glu2093.50 polar lock; and (B) Arg1652.60 and Glu3527.46 salt bridge observed in inactive (green; PDB: 4K5Y) and active (red) states, respectively.

which may lead to receptor activation.37,60−62 Note that the numbering in the superscripts follow the scheme suggested by Wootten et al. for class B GPCRs.62 In all MD simulations, TM2 showed the least fluctuations and constitutes the most invariant region. PMF profiles have been calculated for various collective variables to get insights into the activation of CRF1R. The PMF profile for the TM2− TM6 distance in the cytoplasmic side versus the distance observed for the His1552.50 and Glu2093.50 polar interaction is plotted in Figure 1 (the corresponding profile for the antagonist-bound receptor is shown in Figure S4, SI). The distance between TM2 and TM6 is plotted as a measure of the cytoplasmic displacement of TM6. The intracellular TM6 is observed to be highly flexible, displaying large displacements mostly away from the helical bundle. In fact, the outward shifting of TM6 is considered as one of the important requisites for the activation of GPCRs.12,40,63,64 The activation of apoCRF1R from the inactive crystal structure is observed as the receptor transition to a conformational state, representing the shallow free energy basin (Figure 1). The distance of His1552.50−Glu2093.50 polar lock in this putative active state increased from the crystal structure value of 3.2 Å to ∼12 Å reflecting breakage of this lock. In addition, intracellular sides of TM5, TM6, and TM7 in this well are moved outward to ∼2, ∼7, and ∼3 Å, respectively, from the crystal structure position (Figure 2A). The observed intracellular rearrangements in TM5 and TM6 may facilitate the binding of G-proteins due to the widening of the binding pocket. Besides, intracellular TM3 is bent toward TM4 at Gly2083.49 with an angle of ∼31°, and a twisting in this region is observed with an angle of ∼35° at Glu2093.50 (Figure 2A). These rearrangements facilitated the weakening of the polar contact between His1552.50 and Glu2093.50 (Figure 3A). The flexibility toward the cytoplasmic side of TM3 can be attributed to the presence of two glycine residues (Gly2083.49 and Gly2103.51) around Glu2093.50. Further, interactions that connect intracellular TM3, TM5, TM6, and TM7 are lost in transition of the receptor to its active state. These include hydrophobic contacts between Ile2173.58− Ile2905.57 and Leu2945.61−Val3136.39 along with cation−pi interaction between Lys3116.37−Phe3657.59 (Figure S6, SI). We note that Leu2945.61 and Lys3116.37 are part of LxxxL and KxxK motifs, respectively; these motifs are reported to be essential for G-protein coupling and activation.65−67 Moreover, extracellular halves of TM6 and TM7 moved ∼2 and ∼3 Å, respectively, inward toward the helical bundle (Figure 2B). In these transmembranes, the conserved glycine residues, Gly3246.50 and Gly3567.50, provide structural flexibility. Indeed, introduction of the backbone flexibility in TM7 is recently been

reported to produce an inward shift of the extracellular tip that facilitates the binding of the peptide-agonist.68 We found that the observed inward shift of extracellular TM7 (Figure S7, SI) caused electrostatic interaction between Arg165 2.60 and Glu3527.46 in the apo-active state of the receptor (Figure 3B). Taken together, these rearrangements in TM3, TM5, TM6, and TM7 combined with the breakage of the polar lock describes the putative active state of the receptor. It should be noted that the active state of the receptor identified in this work is based on the two distances, viz., intracellular TM2−TM6 and His1552.50−Glu2093.50. Allosteric Communication Pathways and Pipelines in CRF1 Receptor. To find the allosteric communication pathways and associated pipelines for inactive and active states, we have used mutual information (MI) in torsional degrees of freedom obtained from the extensive molecular dynamics simulations. We have constructed a network of allosteric pathways considering only the top 500 MI pairs having >10 Å distance between Cα atoms. The shortest paths calculated using a graph theory algorithm are then clustered to generate the allosteric pipelines, as described in Methods section. The significant overlap of the multiple allosteric communication pathways designate the allosteric pipelines, and allosteric hubs are the residues in the allosteric pipelines that mediates multiple pathways, as described by Bhattacharya and coworkers.55 The active state shows increase in the number of allosteric pipelines as compared to the inactive state, but most pipelines showed a lower population of allosteric pathways than that of inactive state pipelines (Table S1, SI). The thick red pipeline passing through TM2−TM3−TM5 in the inactive state is otherwise absent in the active state. Further, the intracellular TM region of TM2, TM3, and TM6 is connected via a thick green pipeline (Figure 4A,C). These intracellular pipelines stabilize the receptor in its inactive conformations. Increased allosteric communication in the extracellular region and extracellular part of the TM core was noted in the active state. The allosteric communication among the extracellular loops and extracellular TM7 was observed via blue pipelines which could help in the agonist binding (Figure 4B). This pipeline is further connected to the green pipeline, which is passing through the central region of TM6 and the intracellular parts of TM6 and TM7 in the active state. Moreover, a thin marine pipeline is observed between intracellular TM5 and TM6. The green pipeline together with this thin marine pipeline could contribute in stabilizing the receptors’ active conformation (Figure 4B,D). The important residues found in major allosteric pipelines, which preferably stabilize the inactive state, are Tyr3467.40, E

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

interactions with nearby residues from TM6 and the hydrophobic interaction of this residue with the cocrystallized ligand is also reported in the crystal structure of CRF1R.37 The presence of an aromatic or bulky aliphatic side chain at Phe2033.44 is found to be important in ligand binding and receptor function.69 We observed that this residue resides in the allosteric pipelines of both the active and inactive states of the receptor, though the number of pathways passing through this residue is significantly large in the active state. Asn2835.50, a highly conserved residue from TM5, forms an essential hydrogen bond with the pyridine nitrogen (Figure S8, SI) of the antagonist, and its mutation to alanine results in a complete loss of antagonist binding.37 In the inactive state, Asn2835.50 is present in the thick red pipeline passing through TM2−TM3− TM5 (Figure 4A), whereas in the active state, the residue is present in the thin slate pipeline (Figure 4B). Arg1652.60, Asn2023.43, and Gln3557.49 were suggested as important in the stabilization of receptor in the absence of an agonist and in early steps of receptor activation.70 A study performed on VPAC1R, a class B GPCR, suggest that glutamine of TM7 (Gln3557.49 in CRF1R) interacts with asparagine (Asn2023.43 in CRF1R) located in the central region of the receptor during activation.70 However, mutations of asparagine and glutamine do not affect the binding affinity of agonist.62,70 Asn2023.43 and Gln3557.49 are present in the thick red and green pipelines (Figure 4B), respectively. Our results suggest that Arg1652.60−Glu3527.46 electrostatic interaction helps in the stabilization of the apo-active state of CRF1R. We propose that the polar residues Arg1652.60 and Glu3527.46 will contribute toward efficient binding of the peptide agonist to the receptor. Arg1652.60 mutation supports this fact; however, mutational data for Glu3527.46 is still warranted. In the active state, the number of allosteric communication pathways mediated by residues Arg1652.60 , Asn202 3.43 , Phe2033.44, Glu3527.46, and Gln3557.49 is significantly higher, showing that these residues stabilize the active state by mediating allosteric communication from extracellular side to the intracellular G-protein binding site (Figure S9B and S9D, SI). The number of allosteric communication pathways mediated by residues His1552.50, Glu2093.50, Thr3166.42, and Asn2835.50 is higher in the inactive state of CRF1R, showing that these residues stabilize the inactive state by mediating allosteric communication throughout the receptor even in the absence of any bound ligand (Figure S9A and S9C, SI). Rearrangement of the Residue Interaction Network in Receptor Activation. The dynamic network for inactive and active states are examined through community network analysis (see Methods) of apo-CRF1R. Significant differences are detected between the networks of inactive and active states. In the inactive state, highly connected networks of residues exist in the central and the intracellular parts of TM2, TM3, TM5, TM6, and TM7. These networks are important to keep the receptor in its inactive conformation. Communities formed by the residues in this network are connected via Glu2093.50− Thr3166.42, Asn2835.50−Phe2073.48, Phe1622.57−Asn2023.43, Phe2033.44−Val2795.56, Phe2033.44−Tyr3276.53, Phe1622.57− Gln3557.49, and Thr3266.52−Leu3517.45. Communities near the vicinity of the polar lock (between His1552.50 and Glu2093.50) are modulated in the apo receptor as the conformation changes from inactive to active state. His1552.50 and Glu2093.50 share the same community in the inactive state, whereas they are members of different communities in the active state (Figure 5A,B). This is particularly so because of the loss of polar

Figure 4. Comparison of allosteric pipelines in the inactive (A and C) and active (B and D) states: side (top) and intracellular (bottom) views. Allosteric pipelines are represented by different colors, and the thickness of a pipeline is proportional to the number of neighboring allosteric pathways belongs to that pipeline.

Arg341 (EL3), Tyr267 (EL2), Glu209 3.50 , Phe358 7.52 , Thr3166.42, His1993.40, His1552.50 and Met2765.43. These residues are considered important due to the fact that the number of allosteric pathways passing through them in the inactive state is more (>25) as compared to the active state. Similarly, key residues identified by the allosteric model, stabilizing the active state are Phe3657.59, Phe2033.44, Ala3126.38, Tyr3276.53, Trp2464.60, Pro2394.53, and Cys2504.64. Mutational studies37,60−62,71,72,76 performed on class B GPCRs show the importance of some residues in the activation/inactivation of the receptor. However, these studies could not provide a detailed mechanistic role of these residues. Our results suggest that these residues are located in the allosteric pipelines and thus are responsible to stabilize one state over the other. The pathways passing through His1552.50, Glu2093.50 and Thr3166.42 are stabilizing the inactive state of the receptor. Also, mutational data suggests that these residues stabilize inactive state hence their mutation leads to increased constitutive activity.37,60−62,71,72,76 His1552.50 and Glu2093.50 are present in the thick red pipeline passing through TM2−TM3− TM5 in the inactive state, whereas Thr3166.42 is present in the red subpipeline connecting TM3−TM6 (Figure 4A). Phe2033.44 acts as a central hub that mediates the signals from extracellular side to the intracellular side. Experimental mutational data suggest that mutating this residue to alanine drastically reduces the binding of both the antagonist as well as agonist.69 Phe2033.44 forms crucial hydrophobic−aromatic F

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. Dynamic networks identified from aMD simulations of apo-CRF1R through community network analysis: (A) inactive-state, (B) activestate portrayed on the crystal structure (PDB: 4K5Y). Network communities are colored separately by their identity number. Critical nodes connecting neighboring communities are represented as spheres and the connecting edges are depicted as red lines with their width weighted by betweenness centrality.

Figure 6. Motions (porcupine plots) along the first principal component (PC1) calculated from aMD trajectories of antagonist-bound (left) and apo (right) structures (intracellular view).

interaction between His1552.50 and Glu2093.50. Further, intracellular communities of TM3 and TM6 are strongly communicating via Glu2093.50−Thr3166.42 and thus restricting the rearrangements in these regions that favor the inactive conformation of the receptor (Figure 5A). However, this strong communication at the intracellular side is absent in the active state (Figure 5B) as a result of the outward movement of cytoplasmic TM6 during activation of the receptor. Interestingly, we observed that Phe2033.44 is present at an important central junction in both the inactive and active states. In the inactive state, this TM3 residue forms strong communication with residues from TM5 (Val2795.56) and TM6 (Tyr3276.53). In the active state, Phe2033.44 is connected to TM6 residues Tyr3276.53 and Met3286.54; Tyr3276.53 from this transmembrane is further connected to TM7 via Gln3557.49. Mutation of Phe2033.44 to alanine is reported to reduce the binding affinities of antagonist as well as agonist and the ability of the agonist to activate CRF1R.69 In fact, Phe2033.44 is found as a central hub in the allosteric pipelines of both the inactive and active states, as discussed in the previous section. Overall, the intracellular region of TM2, TM3, TM5, and TM6 shows a very weak communication in the active state as compared to the inactive state.

The community-based analysis therefore reveals that the network toward the intracellular side is significantly weakened and dynamically modulated during activation of the receptor. This is realized due to the changes in polar and hydrophobic interactions in TM2, TM3, TM5, and TM6, thereby making a shift of the conformational population toward the active state of the receptor. Global Conformational Dynamics of the CRF1 Receptor. To reveal functionally important collective motions with large amplitudes, principal component analysis (PCA) was performed on aMD trajectories of the antagonist-bound and the apo receptor. The first two modes, PC1 and PC2, represent a significant contribution to the total fluctuations; PC1 contributes most significantly to global motions of the receptor (Table 2). The free energy landscape obtained considering PC1 and PC2 as reaction coordinates shows that the inactive to active transition requires contribution from both these global modes (Figure S10, SI). The energy landscape corresponding to the apo state of the receptor consists of four types of clusters or minima, viz., “a”, “b”, “c”, and “d” (Figure S10, SI). In cluster a, there are moderate conformational motions for both the collective variables (i.e., His1552.50−Glu2093.50) and TM2− TM6 distances; in cluster b, the major movements are observed in the intracellular TM6 region of the receptor; receptor G

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

The two highly conserved glycine residues, Gly3246.50 in TM6 and Gly3567.50 in TM7, provide structural flexibility to these helices. The antagonist CP-376395 interacts with residues on both sides of the kink induced by Pro3216.47 and Gly3246.50 in TM6.37 Mutagenesis data supports the involvement of prolines (from TM6) in agonist binding and receptor activation.73 We note that the removal of antagonist in general facilitates higher movements in the intracellular and extracellular halves of TM6. In our apo trajectory of the receptor, the cytoplasmic half of TM6 mostly adapts conformations that are away from the helical bundle. However, inward movements of this TM region toward the helical bundle are also observed in initial conformations representing the apo-inactive state compared to the crystal structure (Figure 1). This suggests that the “breathing motion” of the intracellular TM6 is required for the receptor to adapt different inactive conformations − “bound” and “unbound”. Our simulation results also showed that the extracellular end of TM7 is inward tilted (∼3 Å) in the active state. Indeed, this tilt assists TM7 to come closer to the peptide agonist, aided by the flexibility provided by Gly3567.50, and thus improve hydrophobic and polar interactions between their side chains.68 Notably, a two-domain model proposed for class B GPCRs suggest that the binding of the C-terminus of the peptide agonist to N-terminal extracellular domain of the receptor provides an affinity trap. This is subsequently followed by the binding of N-terminus of the ligand to the juxtamembrane domain (J-domain consists of seven TM helices with intervening loops and C-terminal), and thereby activates the receptor.33 Analysis of the residue interaction network suggests that communities strongly connecting TM3, TM5, and TM6 stabilize the inactive conformations. In the inactive state the community shared by the intracellular domain of TM2 and TM3 is strongly connected to TM6 through Glu2093.50− Thr3166.42 interaction (Figure 5A). The strong network of residues in TM2, TM3, and TM6 is dynamically modulated and weakened toward the intracellular side when the receptor shifts from inactive to active state. Twisting of intracellular TM3 and outward movements of intracellular TM5 and TM6 are few important characteristic features of the active state (Figure 2). These rearrangements in TM5 and TM6 create a larger cavity that facilitates coupling of G-proteins (Figure 2). In fact, opening of the G-protein coupling site is suggested as fundamental for GPCR activation.12,63,64 No crystal structure of any class B GPCR coupled to G-protein is available and hence we have compared the volume at the G-protein binding site in the crystal structure of the β2 adrenergic receptor-Gs protein complex (PDB: 3SN6)15 with the binding-pocket in the putative active state. The volume at the G-protein binding site in β2AR is 1126 Å3 and in the present study it is 1070 Å3. However, the binding-pocket volume in the inactive crystal structure of CRF1R is found to be relatively small (307 Å3). The calculated volumes thus suggest that the binding site is large enough to accommodate the G-proteins. Here, the POVME program74,75 was used to calculate the binding-pocket volumes of the G-protein binding site in the receptors (see SI text). Results from this study suggest that initially the G-protein binding site adopts a relatively open conformation owing to the outward movements in the intracellular TM6. It should be noted that the outward movements of intracellular TM6 have also been observed in the other two aMD trajectories (Sim-II and Sim-IV; Table 1), which provide support to our belief that the rearrangements in intracellular TM6 are required as the

conformations in the cluster c are very close to the inactive crystal structure (PDB: 4K5Y); and in cluster d, there are large movements in both the collective variables including the breaking of the polar lock and the outward movements of the intracellular TM6 (Figure S11, SI). Displacement vectors of the collective motions along PC1 for the antagonist-bound and the apo form show notable differences in directions and amplitudes of motions of Cα atoms in TM5, TM6, and TM7. PC1 of the apo trajectory captures the movement of intracellular TM5 and TM6 away from the helical core (Figure 6). Indeed, outward displacements of these TMs are identified as crucial for the receptor activation. Moreover, extracellular halves of TM6 and TM7 show inward movements toward the helical bundle along PC1 (Figure S12, SI). These movements in PC1 toward the extracellular side facilitate the binding of peptide agonist.68 PC2 is the second largest contributing mode to global fluctuations of the receptor that cause bending of cytoplasmic TM3 (Figure S13, SI) and facilitates the breakage of the important “polar lock” between His1552.50 and Glu2093.50 during activation.



DISCUSSION The results of this study signify successful application of an advanced sampling method, viz., aMD for capturing long time scale (of the order of milliseconds) events of functional importance in membrane-bound proteins. Our simulations provided detailed insights into the dynamics of CRF1 receptor, a class B GPCR. Specific modulations in various regions of the receptor are observed in transition from inactive to active state and thereby helped in the characterization of the active state of the receptor. Residues that are highly conserved in the 7TM core of class B receptors are found to be crucial to maintain the structural integrity and/or in the activation of CRF1R. Hydrogenbonding interaction network that exists between TM1 and TM7 in the crystal structure, due to the conserved residues Ser1301.50, Ser3537.47, Gly3567.50, and Phe3577.51, is reported as important for the stabilization of the kink in TM7.30,37 Apart from the interaction of Ser1301.50 with the backbone carbonyl of TM7 at Ser3537.47, which is found to be highly dynamic in nature, we observed loss of other interactions in the active state of the receptor with introduction of added flexibility in the kinked region. Mutation of the conserved histidine residue in TM2 (His1552.50 of CRF1R) is reported to induce the ligandindependent, constitutive activation of class B receptors.62,71,72 Indeed, our simulation results showed that the breakage of the polar interaction between His1552.50−Glu2093.50 is essential for the receptor to adapt the active state (Figure 3). The side chain and the backbone of Trp2364.50, an important residue of the conserved motif GWGxP in TM4, form hydrogen bonds with Asn1572.52 and Tyr1973.38, respectively. This interaction network in fact provides stability to the configuration of TM2, TM3, and TM4.37,39 We found that these hydrogenbonding interactions remain unaltered in the active state of the receptor. This indicates that they are not involved in the activation of the receptor; rather, they are essential for maintaining the structural integrity of the receptor. Further, the free energy surface along χ1 and χ2 dihedral angles of Trp2364.50, show that the inactive and active conformations overlap in the same region centered at χ1 and χ2 dihedral angles of −72° and 100°, respectively (Figure S14, SI). Clearly, this rules out the possibility of the conserved tryptophan in TM4 to act as a rotamer toggle switch in receptor activation. H

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 7. Direct (left) and water-mediated (right) hydrogen-bonding interactions between Glu2093.50 and Thr3166.42 observed in the apo-inactive state of the receptor.

between the communities of TM3 and TM6 via Glu2093.50− Thr3166.42 in the apo-inactive state (Figure 5A). In this context, it is noteworthy to mention that in class A GPCRs, mutations of residues involved in “ionic lock”, viz., the conserved D(E)RY motif in TM3 and Glu6.30/Thr6.34 in TM6, often results in constitutive activation of the receptors.40,77 Further, the water density observed around the conserved polar residues Arg1652.60, Asn2023.43 and Gln3557.49 is reduced upon activation as calculated using Volmap tool in VMD58 (Figure S16, SI). The water molecules observed in the 5 Å region around these three conserved residues in inactive state are 12− 14, which is reduced to 6−8 in the active state. In summary, extensive molecular dynamics simulations performed in this study provided fundamental insights into the activation of CRF1 class B GPCR. Results of this study are crucial for a better understanding of the structural dynamics of CRF1R at the atomistic level, and the information may be judiciously utilized for rational drug design. The structural features that characterize the active state of the receptor are found well correlated with the available experimental data for class B GPCRs, and furthermore, our results may assist in understanding the activation of other class B GPCRs. Nonetheless, the results may be affected by the limitations of the force-field-based approaches and complexity of these biological systems. Other studies using ternary complexes of receptor (with the complete N-domain), agonist, and Gproteins in a hydrated lipid bilayer can be performed to gather additional insights into the activation mechanism of class B GPCRs.

early events of receptor activation (Figure S15, SI). This is followed by the breakage of His1552.50−Glu2093.50 interaction during activation; the putative active state of CRF1R is found stabilized by the electrostatic interaction between polar residues Arg1652.60 and Glu3527.46 (Figure 3B). The water density observed around the polar residues Arg1652.60, Asn2023.43, Glu3527.46, and Gln3557.49 is significant in the inactive cluster than in active cluster. Upon activation, the water present in the cavity is expelled out, moving the side chains of Arg1652.60 and Glu3527.46 within the direct hydrogen bonding distances (Figures 3B and S16, SI). Due to their location in the extracellular TM domain, it is highly probable that these charged residues would prefer to interact with the residues of the peptide agonist in the agonist-bound CRF1R. However, it should be noted that molecular dynamics simulations are prone to errors from different sources and hence the actual activation mechanism may differ from that proposed in this work. The analysis of the allosteric communication pipelines in CRF1R suggest that the active state has an increase in the number of allosteric pipelines than the inactive state, but in the active state, the pipelines are thin due to a lower population of allosteric pathways (Figure 4 and Table S1, SI). His1552.50, Glu2093.50, Thr3166.42 , and Asn2835.50 stabilize the inactive state by mediating more number of allosteric communication pathways through these residues even in the absence of any bound ligand (Figures S9C, SI). In the active state, the number of allosteric communication pathways preferably is mediated by residues Arg1652.60, Asn2023.43, Phe2033.44, Glu3527.46, and Gln3557.49 and thus facilitate the allosteric communication between the extracellular sides to the intracellular G-protein binding site required for receptor’s activation (Figures S9B, SI). In the apo-inactive state, direct as well as water-mediated hydrogen bonding interactions are observed between Glu2093.50 and Thr3166.42 with a population of ∼45% and ∼40%, respectively (Figures 7); however, neither of them is observed in the active state of the receptor. In fact, a recent study on GLP1R, another class B GPCR, proposed that the water-mediated hydrogen-bonding interaction between Glu3.50 and Thr6.42 stabilize the receptor in its inactive conformation.62 Furthermore, the mutation of Thr6.42 in PTH and SecR class B receptors was observed to lead to constitutive activation.71,72,76 Our results clearly demonstrate that the hydrogen-bonding interaction identified here between Glu2093.50 and Thr3166.42 residues restrict the movements of TM3 and TM6, and thereby helps apo-CRF1R to adapt inactive conformations. This is also consistent with the observation of a strong flow of information



ASSOCIATED CONTENT

S Supporting Information *

SI text, Figures S1−S16 and Table S1, providing additional details concerning the results obtained in this study. This material is available free of charge via the Internet at http:// pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We gratefully acknowledge the financial support from IISER Bhopal. This work was performed using the high performance I

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(18) Kohlhoff, K. J.; Shukla, D.; Lawrenz, M.; Bowman, G. R.; Konerding, D. E.; Belov, D.; Altman, R. B.; Pande, V. S. Cloud-based simulations on Google Exacycle reveal ligand modulation of GPCR activation pathways. Nat. Chem. 2014, 6 (1), 15−21. (19) Manglik, A.; Kobilka, B. The role of protein dynamics in GPCR function: insights from the β2AR and rhodopsin. Curr. Opin. Cell Biol. 2014, 27, 136−143. (20) Sobhia, M. E.; Singh, R.; Kare, P.; Chavan, S. Rational design of CCR2 antagonists: a survey of computational studies. Expert Opin. Drug Discovery 2010, 5 (6), 543−557. (21) Singh, R.; Sobhia, M. E. Structure prediction and molecular dynamics simulations of a G-protein coupled receptor: Human CCR2 receptor. J. Biomol. Struct. Dyn. 2013, 31 (7), 694−715. (22) Niesen, M. J.; Bhattacharya, S.; Vaidehi, N. The role of conformational ensembles in ligand recognition in G-protein coupled receptors. J. Am. Chem. Soc. 2011, 133 (33), 13197−13204. (23) Kim, T. H.; Chung, K. Y.; Manglik, A.; Hansen, A. L.; Dror, R. O.; Mildorf, T. J.; Shaw, D. E.; Kobilka, B. K.; Prosser, R. S. The role of ligands on the equilibria between functional states of a G proteincoupled receptor. J. Am. Chem. Soc. 2013, 135 (25), 9465−9474. (24) Costanzi, S. Modeling G protein-coupled receptors and their interactions with ligands. Curr. Opin. Struct. Biol. 2013, 23 (2), 185− 190. (25) Venkatakrishnan, A.; Deupi, X.; Lebon, G.; Tate, C. G.; Schertler, G. F.; Babu, M. M. Molecular signatures of G-proteincoupled receptors. Nature 2013, 494 (7436), 185−194. (26) Chung, K. Y.; Rasmussen, S. G.; Liu, T.; Li, S.; DeVree, B. T.; Chae, P. S.; Calinski, D.; Kobilka, B. K.; Woods, V. L., Jr; Sunahara, R. K. Conformational changes in the G protein Gs induced by the β2 adrenergic receptor. Nature 2011, 477 (7366), 611−615. (27) Wootten, D.; Christopoulos, A.; Sexton, P. M. Emerging paradigms in GPCR allostery: implications for drug discovery. Nat. Rev. Drug Discovery 2013, 12 (8), 630−644. (28) Rosenbaum, D. M.; Zhang, C.; Lyons, J. A.; Holl, R.; Aragao, D.; Arlow, D. H.; Rasmussen, S. G.; Choi, H.-J.; DeVree, B. T.; Sunahara, R. K. Structure and function of an irreversible agonist-β2 adrenoceptor complex. Nature 2011, 469 (7329), 236−240. (29) Elgeti, M.; Rose, A. S.; Bartl, F. J.; Hildebrand, P. W.; Hofmann, K.-P.; Heck, M. Precision vs Flexibility in GPCR signaling. J. Am. Chem. Soc. 2013, 135 (33), 12305−12312. (30) Hollenstein, K.; de Graaf, C.; Bortolato, A.; Wang, M.-W.; Marshall, F. H.; Stevens, R. C. Insights into the structure of class B GPCRs. Trends Pharmacol. Sci. 2014, 35 (1), 12−22. (31) Poyner, D. R.; Hay, D. L. Secretin family (Class B) G proteincoupled receptors−from molecular to clinical perspectives. Br. J. Pharmacol. 2012, 166 (1), 1−3. (32) Sloop, K. W.; Willard, F. S.; Brenner, M. B.; Ficorilli, J.; Valasek, K.; Showalter, A. D.; Farb, T. B.; Cao, J. X.; Cox, A. L.; Michael, M. D. Novel small molecule glucagon-like peptide-1 receptor agonist stimulates insulin secretion in rodents and from human islets. Diabetes 2010, 59 (12), 3099−3107. (33) Hoare, S. R. Mechanisms of peptide and nonpeptide ligand binding to Class B G-protein-coupled receptors. Drug Discovery Today 2005, 10 (6), 417−427. (34) Siu, F. Y.; Stevens, R. C. RAMP-ing up class-B GPCR ECD structural coverage. Structure 2010, 18 (9), 1067−1068. (35) Dong, M.; Koole, C.; Wootten, D.; Sexton, P.; Miller, L. Structural and functional insights into the juxtamembranous aminoterminal tail and extracellular loop regions of class B GPCRs. Br. J. Pharmacol. 2014, 171 (5), 1085−1101. (36) Siu, F. Y.; He, M.; de Graaf, C.; Han, G. W.; Yang, D.; Zhang, Z.; Zhou, C.; Xu, Q.; Wacker, D.; Joseph, J. S. Structure of the human glucagon class B G-protein-coupled receptor. Nature 2013, 499, 444− 449. (37) Hollenstein, K.; Kean, J.; Bortolato, A.; Cheng, R. K.; Doré, A. S.; Jazayeri, A.; Cooke, R. M.; Weir, M.; Marshall, F. H. Structure of class B GPCR corticotropin-releasing factor receptor 1. Nature 2013, 499, 438−443.

computing resources of IISER Bhopal. R.S. and N.A. are supported by the research fellowships provided by IISER Bhopal.



REFERENCES

(1) Pierce, K. L.; Premont, R. T.; Lefkowitz, R. J. Seventransmembrane receptors. Nat. Rev. Mol. Cell Biol. 2002, 3 (9), 639−650. (2) Rosenbaum, D. M.; Rasmussen, S. G.; Kobilka, B. K. The structure and function of G-protein-coupled receptors. Nature 2009, 459 (7245), 356−363. (3) Vaidehi, N.; Pease, J. E.; Horuk, R. Modeling Small MoleculeCompound Binding to G-Protein-Coupled Receptors. Methods Enzymol. 2009, 460, 263−288. (4) Altenbach, C.; Kusnetzow, A. K.; Ernst, O. P.; Hofmann, K. P.; Hubbell, W. L. High-resolution distance mapping in rhodopsin reveals the pattern of helix movement due to activation. Proc. Natl. Acad. Sci. U.S.A. 2008, 105 (21), 7439−7444. (5) Liu, J. J.; Horst, R.; Katritch, V.; Stevens, R. C.; Wüthrich, K. Biased signaling pathways in β2-adrenergic receptor characterized by 19F-NMR. Science 2012, 335 (6072), 1106−1110. (6) Yao, X.; Parnot, C.; Deupi, X.; Ratnala, V. R.; Swaminath, G.; Farrens, D.; Kobilka, B. Coupling ligand structure to specific conformational switches in the β2-adrenoceptor. Nat. Chem. Biol. 2006, 2 (8), 417−422. (7) Hoare, S. R.; Fleck, B. A.; Gross, R. S.; Crowe, P. D.; Williams, J. P.; Grigoriadis, D. E. Allosteric ligands for the corticotropin releasing factor type 1 receptor modulate conformational states involved in receptor activation. Mol. Pharmacol. 2008, 73 (5), 1371−1380. (8) Kobilka, B. K.; Deupi, X. Conformational complexity of Gprotein-coupled receptors. Trends Pharmacol. Sci. 2007, 28 (8), 397− 406. (9) Cherezov, V.; Rosenbaum, D. M.; Hanson, M. A.; Rasmussen, S. G. F.; Thian, F. S.; Kobilka, T. S.; Choi, H. J.; Kuhn, P.; Weis, W. I.; Kobilka, B. K. High-resolution crystal structure of an engineered human β2-adrenergic G protein coupled receptor. Science 2007, 318 (5854), 1258−1265. (10) Katritch, V.; Cherezov, V.; Stevens, R. C. Structure-function of the G protein-coupled receptor superfamily. Annu. Rev. Pharmacol. Toxicol. 2013, 53, 531−556. (11) Rasmussen, S. G.; Choi, H.-J.; Fung, J. J.; Pardon, E.; Casarosa, P.; Chae, P. S.; DeVree, B. T.; Rosenbaum, D. M.; Thian, F. S.; Kobilka, T. S. Structure of a nanobody-stabilized active state of the β2adrenoceptor. Nature 2011, 469 (7329), 175−180. (12) Li, J.; Jonsson, A. L.; Beuming, T.; Shelley, J. C.; Voth, G. A. Ligand-Dependent Activation and Deactivation of the Human Adenosine A2A Receptor. J. Am. Chem. Soc. 2013, 135 (23), 8749− 8759. (13) Miao, Y.; Nichols, S. E.; Gasper, P. M.; Metzger, V. T.; McCammon, J. A. Activation and dynamic network of the M2 muscarinic receptor. Proc. Natl. Acad. Sci. U.S.A. 2013, 110 (27), 10982−10987. (14) Nygaard, R.; Zou, Y.; Dror, R. O.; Mildorf, T. J.; Arlow, D. H.; Manglik, A.; Pan, A. C.; Liu, C. W.; Fung, J. J.; Bokoch, M. P. The Dynamic Process of β2-Adrenergic Receptor Activation. Cell 2013, 152 (3), 532−542. (15) Rasmussen, S. G.; DeVree, B. T.; Zou, Y.; Kruse, A. C.; Chung, K. Y.; Kobilka, T. S.; Kobilka, B. K. Crystal structure of the β2 adrenergic receptor−Gs protein complex. Nature 2011, 477 (7366), 549−555. (16) Dror, R. O.; Arlow, D. H.; Maragakis, P.; Mildorf, T. J.; Pan, A. C.; Xu, H.; Borhani, D. W.; Shaw, D. E. Activation mechanism of the β2-adrenergic receptor. Proc. Natl. Acad. Sci. U.S.A. 2011, 108 (46), 18684−18689. (17) 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. Structural basis for modulation of a G-protein-coupled receptor by allosteric drugs. Nature 2013, 503 (7475), 295−299. J

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (38) Hemley, C.; McCluskey, A.; Keller, P. Corticotropin releasing hormone-a GPCR drug target. Curr. Drug Targets 2007, 8 (1), 105− 115. (39) Bortolato, A.; Doré, A. S.; Hollenstein, K.; Tehan, B. G.; Mason, J. S.; Marshall, F. H. Structure of Class B GPCRs: New Horizons for Drug Discovery. Br. J. Pharmacol. 2014, 171, 3132−3145. (40) Trzaskowski, B.; Latek, D.; Yuan, S.; Ghoshdastider, U.; Debinski, A.; Filipek, S. Action of molecular switches in GPCRstheoretical and experimental studies. Curr. Med. Chem. 2012, 19 (8), 1090−1109. (41) Hamelberg, D.; Mongan, J.; McCammon, J. A. Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules. J. Chem. Phys. 2004, 120, 11919−11929. (42) Pierce, L. C.; Salomon-Ferrer, R.; Augusto F. de Oliveira, C.; McCammon, J. A.; Walker, R. C. Routine access to millisecond time scale events with accelerated molecular dynamics. J. Chem. Theory Comput. 2012, 8 (9), 2997−3002. (43) Eswar, N.; Eramian, D.; Webb, B.; Shen, M.-Y.; Sali, A. Protein structure modeling with MODELLER. Methods Mol. Biol. 2008, 426, 145−159. (44) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79 (2), 926−935. (45) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781−1802. (46) Wang, Y.; Harrison, C. B.; Schulten, K.; McCammon, J. A. Implementation of accelerated molecular dynamics in NAMD. Comput. Sci. Discovery 2011, 4 (1), 015002. (47) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R.; Evanseck, J.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. Allatom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102 (18), 3586−3616. (48) MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved treatment of the protein backbone in empirical force fields. J. Am. Chem. Soc. 2004, 126 (3), 698−699. (49) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr; Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J. Chem. Phys. B 2010, 114 (23), 7830−7843. (50) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31 (4), 671−690. (51) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103 (19), 8577−8593. (52) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23 (3), 327−341. (53) Miao, Y.; Nichols, S. E.; McCammon, J. A. Free energy landscape of G-protein coupled receptors, explored by accelerated molecular dynamics. Phys. Chem. Chem. Phys. 2014, 16 (14), 6398− 6406. (54) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29 (7), 845−854. (55) Bhattacharya, S.; Vaidehi, N. Differences in Allosteric Communication Pipelines in the Inactive and Active States of a GPCR. Biophys. J. 2014, 107 (2), 422−434. (56) Eargle, J.; Luthey-Schulten, Z. NetworkView: 3D display and analysis of protein· RNA interaction networks. Bioinformatics 2012, 28 (22), 3000−3001.

(57) Sethi, A.; Eargle, J.; Black, A. A.; Luthey-Schulten, Z. Dynamical networks in tRNA: protein complexes. Proc. Natl. Acad. Sci. U.S.A. 2009, 106 (16), 6620−6625. (58) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graph. 1996, 14 (1), 33−38. (59) 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. Structure and dynamics of the M3 muscarinic acetylcholine receptor. Nature 2012, 482 (7386), 552−556. (60) Heller, R. S.; Kieffer, T. J.; Habener, J. F. Point mutations in the first and third intracellular loops of the glucagon-like peptide-1 receptor alter intracellular signaling. Biochem. Biophys. Res. Commun. 1996, 223 (3), 624−632. (61) Hoare, S. R.; Brown, B. T.; Santos, M. A.; Malany, S.; Betz, S. F.; Grigoriadis, D. E. Single amino acid residue determinants of nonpeptide antagonist binding to the corticotropin-releasing factor1 (CRF1) receptor. Biochem. Pharmacol. 2006, 72 (2), 244−255. (62) Wootten, D.; Simms, J.; Miller, L. J.; Christopoulos, A.; Sexton, P. M. Polar transmembrane interactions drive formation of ligandspecific and signal pathway-biased family BG protein-coupled receptor conformations. Proc. Natl. Acad. Sci. U.S.A. 2013, 110 (13), 5211− 5216. (63) Park, J. H.; Scheerer, P.; Hofmann, K. P.; Choe, H.-W.; Ernst, O. P. Crystal structure of the ligand-free G-protein-coupled receptor opsin. Nature 2008, 454 (7201), 183−187. (64) Standfuss, J.; Edwards, P. C.; D’Antona, A.; Fransen, M.; Xie, G.; Oprian, D. D.; Schertler, G. F. The structural basis of agonist-induced activation in constitutively active rhodopsin. Nature 2011, 471 (7340), 656−660. (65) Couvineau, A.; Lacapère, J.-J.; Tan, Y.-V.; Rouyer-Fessard, C.; Nicole, P.; Laburthe, M. Identification of cytoplasmic domains of HVPAC1 receptor required for activation of adenylyl cyclase crucial role of two charged amino acids strictly conserved in class II G protein-coupled receptors. J. Biol. Chem. 2003, 278 (27), 24759− 24766. (66) Lee, N. H.; Geoghagen, N.; Cheng, E.; Cline, R. T.; Fraser, C. M. Alanine scanning mutagenesis of conserved arginine/lysinearginine/lysine-XX-arginine/lysine G protein-activating motifs on m1 muscarinic acetylcholine receptors. Mol. Pharmacol. 1996, 50 (1), 140−148. (67) Punn, A.; Chen, J.; Delidaki, M.; Tang, J.; Liapakis, G.; Lehnert, H.; Levine, M. A.; Grammatopoulos, D. K. Mapping structural determinants within third intracellular loop that direct signaling specificity of type 1 corticotropin-releasing hormone receptor. J. Biol. Chem. 2012, 287 (12), 8974−8985. (68) Coin, I.; Katritch, V.; Sun, T.; Xiang, Z.; Siu, F. Y.; Beyermann, M.; Stevens, R. C.; Wang, L. Genetically Encoded Chemical Probes in Cells Reveal the Binding Path of Urocortin-I to CRF Class B GPCR. Cell 2013, 155 (6), 1258−1269. (69) Spyridaki, K.; Matsoukas, M.-T.; Cordomi, A.; Gkountelias, K.; Papadokostaki, M.; Mavromoustakos, T.; Logothetis, D. E.; Margioris, A. N.; Pardo, L.; Liapakis, G. Structural-Functional Analysis of the Third Transmembrane Domain of the CRF1 Receptor: Role in Activation and Allosteric Antagonism. J. Biol. Chem. 2014, 289 (27), 18966−18977. (70) Chugunov, A. O.; Simms, J.; Poyner, D. R.; Dehouck, Y.; Rooman, M.; Gilis, D.; Langer, I. Evidence that interaction between conserved residues in transmembrane helices 2, 3, and 7 are crucial for human VPAC1 receptor activation. Mol. Pharmacol. 2010, 78 (3), 394−401. (71) Gaudin, P.; Maoret, J.-J.; Couvineau, A.; Rouyer-Fessard, C.; Laburthe, M. Constitutive activation of the human vasoactive intestinal peptide 1 receptor, a member of the new class II family of G proteincoupled receptors. J. Biol. Chem. 1998, 273 (9), 4990−4996. (72) Schipani, E.; Langman, C.; Parfitt, A.; Jensen, G.; Kikuchi, S.; Kooh, S.; Cole, W.; Jüppner, H. Constitutively activated receptors for parathyroid hormone and parathyroid hormone−related peptide in Jansen’s metaphyseal chondrodysplasia. N. Engl. J. Med. 1996, 335 (10), 708−714. K

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (73) Conner, A. C.; Hay, D. L.; Simms, J.; Howitt, S. G.; Schindler, M.; Smith, D. M.; Wheatley, M.; Poyner, D. R. A key role for transmembrane prolines in calcitonin receptor-like receptor agonist binding and signalling: implications for family B G-protein-coupled receptors. Mol. Pharmacol. 2005, 67 (1), 20−31. (74) Durrant, J. D.; de Oliveira, C. A. F.; McCammon, J. A. POVME: an algorithm for measuring binding-pocket volumes. J. Mol. Graph. Model. 2011, 29 (5), 773−776. (75) Durrant, J. D.; Votapka, L.; Sørensen, J.; Amaro, R. E. POVME 2.0: An Enhanced Tool for Determining Pocket Shape and Volume Characteristics. J. Chem. Theory Comput. 2014, 10 (11), 5047−5056. (76) Ganguli, S. C.; Park, C.-G.; Holtmann, M. H.; Hadac, E. M.; Kenakin, T. P.; Miller, L. J. Protean effects of a natural peptide agonist of the G protein-coupled secretin receptor demonstrated by receptor mutagenesis. J. Pharmacol. Exp. Ther. 1998, 286 (2), 593−598. (77) Rovati, G. E.; Capra, V.; Neubig, R. R. The highly conserved DRY motif of class A G protein-coupled receptors: beyond the ground state. Mol. Pharmacol. 2007, 71 (4), 959−964.

L

DOI: 10.1021/jp509814n J. Phys. Chem. B XXXX, XXX, XXX−XXX