Simulations of biased agonists in the β(2) adrenergic receptor with accelerated molecular dynamics.

The biased agonism of the G protein-coupled receptors (GPCRs), where in addition to a traditional G protein-signaling pathway a GPCR promotes intracellular signals though β-arrestin, is a novel paradigm in pharmacology. Biochemical and biophysical studies have suggested that a GPCR forms a distinct ensemble of conformations signaling through the G protein and β-arrestin. Here we report on the dynamics of the β2 adrenergic receptor bound to the β-arrestin and G protein-biased agonists and the empty receptor to further characterize the receptor conformational changes caused by biased agonists. We use conventional and accelerated molecular dynamics (aMD) simulations to explore the conformational transitions of the GPCR from the active state to the inactive state. We found that aMD simulations enable monitoring of the transition within the nanosecond time scale while capturing the known microscopic characteristics of the inactive states, such as the ionic lock, the inward position of F6.44, and water clusters. Distinct conformational states are shown to be stabilized by each biased agonist. In particular, in simulations of the receptor with the β-arrestin-biased agonist N-cyclopentylbutanepherine, we observe a different pattern of motions in helix 7 when compared to simulations with the G protein-biased agonist salbutamol that involves perturbations of the network of interactions within the NPxxY motif. Understanding the network of interactions induced by biased ligands and the subsequent receptor conformational shifts will lead to development of more efficient drugs.

G protein-coupled receptors (GPCRs) are the largest group of cell surface proteins, which translate chemical information of various extracellular stimuli to a specific biological response in cells. Clinically, GPCR ligands regulate many physiological and pathological conditions and represent a remarkable class of pharmaceutical agents, with 26.8% of FDA-approved drugs targeting group A of GPCRs.[1]

Classically, the binding of an agonist to an inactive state of a GPCR causes conformational changes that lead to an active state of a GPCR, which is capable of activating a heterotrimeric G protein.[2,3] The G protein, in turn, promotes signaling by modulating the activity of effector enzymes, leading to the generation of second messengers. Termination of the signal, which is known as receptor desensitization, results in uncoupling of the receptor from the G protein through phosphorylation of the receptor C terminus and recruitment of β-arrestin to the phosphorylated C terminus of GPCRs, where β-arrestin binding stimulates receptor endocytosis.[4,5]

Recent pharmacological studies on a large number of GPCRs have shown that agonist-bound GPCRs promote signaling mechanisms in cells independent of G protein.[6−12] Thus, it has been found that β-arrestin, in addition to receptor desensitization, transmits unique signals to catalytically active proteins.[9−11] The biased receptor signaling phenomenon, also known as functional selectivity, has classified GPCRs binders into biased and unbiased agonists or antagonists. A biased ligand is capable of activating or deactivating one specific signaling pathway, whereas an unbiased ligand has relatively equal impact on multiple signaling pathways. β-Arrestin and G protein-biased signaling has been linked to side effects of several GPCR drugs.[13] For example, tachyphylaxis and drug tolerance of GPCR agonists are thought to be due to β-arrestin stimulation.[13] The design of biased agonists may therefore be a novel strategy to improve the efficacy of GPCR agonists in clinics.

Structurally, GPCRs are composed of seven transmembrane spanning helices along with connecting loops. Recently, the β2 adrenergic receptor (β2AR) has been crystallized in both inactive and active states.[14−17] The crystal structures of the active state were first produced in complex with the G protein-like antibody[16] and subsequently in complex with the Gs protein.[17] Both of these structures represent the active state of the receptor stimulating G protein-dependent signaling. Comparison of structures of the active and inactive states shows small changes within the ligand binding site and more profound changes on the intracellular side of the receptor, including the exterior movement of helix 6 and rearrangement of helices 5 and 7.[16,17]

The ability of GPCRs to initiate pharmacologically distinct signaling pathways raises the hypothesis of the existence of distinct active conformations of the receptor. Kinetics measurements of several labeled lysines and cysteines of β2AR in the presence of the agonist (THRX-144877) with high equal stimulation of β-arrestin and G protein showed higher reactivity of Lys3057.32 (the Ballesteros–Weinstein nomenclature[18] in subscription, where the most conserved residue in a given helix, X, is assigned the index X.50 and the other residues of the helix are numbered relative to the 50 position) than in the presence of the agonists with domination of G protein activation.[19] This suggests a different role of helix 7 in the active receptor state for stimulating G protein and β-arrestin.[19] Furthermore, 19F NMR studies of β2AR in complex with various ligands has indicated that ligands with β-arrestin bias affect helix 7 to a greater extent than they do helix 6.[20] Fluorescence spectroscopy studies in the arginine–vasopressin type 2 receptor have shown opposite effects in tryptophan fluorescence, with fluorophores introduced in helices 6, 7, and 8, in the presence of G-protein and β-arrestin-biased agonists.[21] These indirect experiments provide evidence that the active state recruiting β-arrestin is structurally distinct from the active state activating G-protein. However, they cannot reveal what the specific changes in the receptor domains are that are caused by a biased agonist recruiting β-arrestin. Recent work involving mutagenesis and molecular modeling has provided the first attempt in finding such changes. In particular, a sulfur–aromatic interaction between M1343.32 and Y3807.43 has been proposed to stabilize the active state of the cholecystokinin 2 receptor recruiting β-arrestin.[22]

In this study we probe and compare receptor conformational changes caused by biased agonists binding using molecular dynamics (MD) simulations on the example of the available crystal structure of the β2AR active state. Exploration of the β2AR conformational ensemble has been initiated with various MD simulation protocols, involving a combination of coarse-grained and all-atom simulations, adiabatic biased and metadynamics simulations, and microsecond conventional MD simulations.[23−25] Simulations from the inactive-to-active state of the β2AR, using the available crystal structures, in the presence of agonists, partial agonists, and neutral antagonists have demonstrated the different abilities of the ligands to stabilize the active state activating G-protein.[23,24] The active-to-inactive state microsecond simulations have provided details of the β2AR activation mechanism stimulating the G protein.[25] Hereby, we extend the simulations of β2AR in the empty form and bound to biased agonists by employing intensive conventional and accelerated MD simulations (cMD and aMD, respectively) to initiate understanding of functional selectivity at the molecular level.

Salbutamol (Sal) was chosen as an agonist biased to the G-protein and N-cyclopentylbutanepherine (CPB) as an agonist with the domination of β-arrestin bias.[26] Both ligands are partial agonists compared to isoproterenol.[26] Interestingly, the β2AR agonists used in the treatment of asthma cause tachyphylaxis, which is linked to β-arrestin recruitment.[13]

The aMD simulations, in which an external potential is applied to dihedral angles of a protein to overcome energetic barriers, has been recently successfully applied to study large conformational changes of soluble proteins in a shorter time.[27−30] Unlike many other computational methods for enhanced sampling, aMD increases sampling without the use of end states or a predefined “reaction coordinates” and, therefore, allows for greater freedom in sampling diverse reaction pathways.[31,32] In our study we applied aMD simulations to the membrane protein for the first time to compare the dynamics of the β2AR bound to biased agonists and identified conformational changes that lead to activation of distinct signaling pathways in GPCRs.

Application of aMD has allowed the observation of global and local changes during the active-to-inactive state transition of the empty receptor and bound to Sal or CPB within the nanosecond time scale. In addition, aMD has identified distinct conformational ensembles stabilized by the two agonists and highlighted specific changes caused by CPB. The generated receptor conformational ensembles can be helpful in the development of structure-based drug design methodologies to design novel efficient and selective drugs.

Methods

System Preparation

The crystal structure of the active state of the human β2AR (3P0G)[15] was used for molecular docking and simulations. The G protein-like antibody, T4-lysozyme, in the crystal structure was removed, and the intracellular loops were joined together using the Prime 2.2 program (Schrödinger, LLC, New York, NY, USA; Prime., 2008). The receptor structure and ligands (Sal and CPB) for molecular docking and simulations were prepared using the Schrodinger suite with the OPLS2005 force field.[33,34] The agonists were docked to the receptor with the Glide 5.6 program (Schrödinger, LLC; Glide 5.6, 2009), using the extra-precision scoring setting. The docking pose of Sal was selected after comparison with the available crystal structure of the Sal−β1AR complex (PDB code 2Y04). This pose was also the top scored pose in the docking studies. The docking pose of CPB with top scoring was selected for the simulations. In fact, in all poses CPB had contacts with two Ser, providing similar orientations of the catecholamine group and none of CPB docking orientation allows a hydrogen bond with N293. Note, we observed the hydrogen bond interactions with N293 during the MD simulations (Supporting Information, Figure 3S). System Builder of Maestro GUI 9.0 (Schrodinger, LLC; 2010) was used to embed the agonist–receptor complexes in a water–lipid box with dimension 83 × 90 × 106 Å and neutralized with sodium and chloride ions. The resulting systems had 154 lipid molecules, 52 chloride ions, and 40 sodium ions. All atom 1-palmitoyl-2-oleoylphosphatidylcholine bilayer (POPC) and SPC water were used for the lipid and water models, respectively. The average size of simulations systems was 70 000 atoms.

Conventional and Accelerated Molecular Dynamics

All MD simulations were performed using Desmond 2.2[35] with the OPLS2005 force field[33,34] and full particle mesh Ewald electrostatics.[36] Operational parameters have included a 2 fs time step and a 9 Å cutoff for the truncation of nonbonded interactions. The prepared complexes in the hydrated phospholipid bilayer were first minimized with the fixed receptor–ligand complex, applying the conjugate gradient algorithm up to a convergence threshold of 0.5 kcal/mol/Å. In the next step, the systems with the fixed receptor–ligand complexes were heated to 300 K in the NVT ensemble, where each 10 K increase was during 12 ps. Lipids and water were equilibrated over 8 ns at 300 K, followed by 2 ns of equilibration with the restrained protein backbone and 2 ns of equilibration with the fully released system. The equilibrated biomolecular systems were subjected to constant temperature (300 K) and pressure (1 atm) production molecular dynamics with cMD and aMD protocols. Three simulation runs for cMD and aMD protocols were performed, and each run had a random initial velocity assignment (Table 1). The SHAKE algorithm was used to constrain all covalent bonds involving hydrogen atoms.[37]

Table 1

Summary of Conventional MD (cMD) and Accelerated MD (aMD) Simulations

biosystem	length, ns	no. of simulations	av Cα rmsd, Å	av ligand rmsd, Å	av helical content, %	area per lipid headgroup, Å2	lipid bilayer thickness, Å	transition time,a ns
cMD
empty β2AR	∼91	3	2.9 ± 0.1		86 ± 2	51.2 ± 11	45.7 ± 3	41 ± 7
β2AR- CPB	∼220	3	3 ± 0.2	0.6 ± 0.4	89 ± 2	50.2 ± 9	43 ± 3
β2AR- Sal	∼300	3	2.7 ± 0.3	0.7 ± 0.3	82 ± 4	55.4 ± 11	43.4 ± 2
aMD
empty β2AR	∼122	3	3.2 ± 0.2		85 ± 4	60.8 ± 9	43 ± 6	17 ± 4
β2AR- CPB	∼213	3	2.7 ± 0.2	0.5 ± 0.2	87 ± 2	56.5 ± 7	44.4 ± 8	156 ± 7
β2AR- Sal	∼370	3	2.9 ± 0.3	0.7 ± 0.2	86 ± 5	57.4 ± 13	42 ± 6	274 ± 8

The transition time indicates the transition from the active-like conformation to the inactive-like conformation.

The aMD protocol[31,32] with the alteration of the potential energy surface of the biomolecular system was applied to enhance conformational sampling of the systems and, therefore, to enlarge the accessible time scale of cMD. Acceleration comes from a “boost potential”, ΔV(r), which is added to the original dihedral potential, V(r), that increases the energy to V*(r) within basins, using the equationsandwhere E is the threshold dihedral energy specified by the user, which controls the level of the potential surface affected by bias, and α is the acceleration parameter establishing the shape of the modified potential.

Thus, a transition from one state to another occurs with an accelerated rate as a result of propagation of a trajectory on this modified energy surface. The boost dihedral potential was applied on the receptor alone as well as the receptor and lipids together. The average dihedral energy of the system during cMD simulations was used as a reference for computing E and α applying the following equations: E = Vav + Vav × c and α = E/5, where Vav is the average dihedral energy of the system in the initial 5 ns of the cMD and c is constant (0.2, 0.3, 0.4, and 0.5). Several acceleration E and α values were probed (Supporting Information, Table 1S), and E0.3 = Vav + Vav × 0.3 and α = E0.3/5 were chosen for protein and lipids to conduct multiple runs with random velocities for production simulations in all biosystems.

Principal Component Analysis (PCA)

The calculation of the covariance matrix and its diagonalization to get the associated principal components describing the direction and amplitude of receptor motions was performed with Gromacs 4.5.5[38] employing the g_covar module. The principal components (PC) obtained from mapping the receptor Cα atoms of all the simulated trajectories were used to project the trajectories of a particular receptor system. This has allowed comparison of simulation data within a common subspace. To compare cMD and aMD simulations, equal length of the simulations was used for the PCA. The length of the simulated trajectories for the PCA analysis was defined by the active-to-inactive transition time. The sampling frequency of these trajectories was 30 ps. The principal components, PC1 and PC2, provide a distinct separation of the receptor active and inactive states. The two first principal components, representing the largest contribution to the overall fluctuation of the receptor (60–70%), were used to project the simulated ensemble of receptor conformations from all simulated trajectories to visualize the receptor essential motions. The projection of the simulated structures was conducted on individual trajectories using the same Cα atoms. The g_anaeig module of Gromacs 4.5.5 was used for the projection. To avoid general translation and rotation of the receptor, the trajectories were aligned on the basis of the Cα atoms. The available crystal structures of the β2AR with PDB codes 2RH1, 3D4S, 3NY8, 3NY9, 3NYA, 3PDS, 3P0G, and 3SN6 were used to monitor whether the simulated conformational space samples the projection of available experimental structures. The porcupine plot was built on the basis of the extreme conformations derived from the first principal component using the g_anaeig module. The porcupine plot was constructed using the most stable ensemble of the receptor conformations (the third ensemble of conformations was used for the receptor complexes) obtained from the sampling probabilities studies. The porcupine plot was plotted using the VMD script.

Sampling Probabilities

We used the g_sham module of Gromacs 4.5.5 to compute probability sampling (υ) landscape in the PC1–PC2 conformational space from the simulated trajectories identified by g_anaeig. The probability sampling landscape was calculated using the equationwhere kb is the Boltzmann constant, T is 300 K, P(xi) is an estimate of the probability density function obtained from a histogram of the MD data (eigenvalues and eigenvectors, in our case), and Pmax(x) is the probability of the most probable state.

Other Analyses

Most of the trajectory analysis was performed with VMD 1.9.1.[39] The distance and angle cutoffs of 3 Å and 35°, respectively, were used for hydrogen bonds. The helix bending was calculated and plotted using the Bendix application[40] in VMD 1.9.1. Images were captured with VMD 1.9.1, Maestro GUI 9.0, and PYMOL 0.99 (DeLano WL, 2002).

Results

Conventional and Accelerated Molecular Dynamics Simulations

Agonists stabilize the active state of GPCRs. The activated form of GPCRs is then phosphorylated by kinases at the intracellular side of receptors. In turn, the active phosphorylated GPCRs are capable of binding to β-arrestin,[41] which sterically blocks interactions with the G protein and thus suppresses classical G- protein-mediated signaling. In parallel, the receptor in the complex with β-arrestin is active to interact with various signaling proteins, leading to G protein-independent signaling.[11,42,43] It is likely that the G-protein active state of GPCRs transforms to the β-arrestin active state for initiating of G protein-independent signaling rather than the inactive state of the receptor converting directly to the β-arrestin active state.[22] We, therefore, study the dynamics of the available active state of β2AR in the empty form and in the presence of the biased agonists, hoping to capture characteristics of the β-arrestin active state. To initiate our study, we have docked Sal and CPB to the crystal structure of the active state that was available at the time (PDB 3P0G, see details under Methods) and used the obtained complexes for the MD simulations in the water–lipid bilayer. Figure 1 shows the binding mode of Sal and CPB in the β2AR from the docking studies. The nanobody, T4-lysozyme, in the crystal structure was removed, and the intracellular loops were joined together.

Figure 1

Binding mode of salbutamol (Sal) and N-cyclopentylbutanepherine (CPB) in the β2 adrenergic receptor (β2AR).

Conformational dynamics of the β2AR active state in the empty form and bound to the agonists were explored using intensive cMD and aMD simulations, which are summarized in Table 1. aMD simulations with an external potential applied to the dihedral angles of the receptor and lipid atoms (see Methods) were used to enhance conformational changes.[31,32] Wang et al.[44] have shown that aMD simulations reproduce the structural and dynamic properties of lipids and accelerate the equilibration of various membrane bilayers, indicating the applicability of the external potential to the dihedral angles of lipids, in addition to a protein, in the simulations of membrane proteins. To identify a suitable threshold dihedral energy (E) and acceleration parameter (α) for the protein and lipids in the aMD simulations, several E and α values were tested (Supporting Information, Table 1S) and those that reproduced the known conformational changes (see below and the next section), while retaining the receptor fold, were chosen (see Methods).

Biomolecular systems were stable in the cMD and aMD simulations with the chosen E and α as shown by the average root-mean-square deviation of the Cα atoms of the protein and heavy atoms of ligands, as well as the average helical content of the helices during the simulations (Table 1). Notably, the chosen boost potential (see Methods and Table 1S in the Supporting Information) in the aMD did not cause any loss of the secondary structure, suggesting that the conformational changes in the receptor structure discussed below are the result of intrinsic receptor mobility. The area per lipid headgroup and lipid bilayer thickness were stable in the cMD and aMD simulations but were underestimated relative to the experimental values for POPC, that is, 62.7 ± 1.3 Å2 and 39.8 ± 0.8 Å.[45] This could be due to the measurements in an inhomogeneous system containing protein.

To compare the conformational sampling of cMD and aMD simulations, we conducted a PCA[46] by projecting the simulated trajectories of each receptor system onto the plane defined by the two lowest PCs obtained from mapping all of the simulated trajectories of β2AR (See Methods). The available crystal structures were also projected onto the simulated space to compare the performance of cMD and aMD methods in sampling of crystal structure projections. The motions along these two lowest PCs separate the receptor crystal structures of the active state from the inactive state as shown in Figure 2. In particular, the major change between the inactive and active states, which is the outward movement of helix 6, is captured by both PCs. The PC plots indicate that all aMD simulations sample the projections of the crystal structures, unlike cMD simulations. Thus, the active state of the receptor in the empty form and bound with the agonists reaches the projections of the crystal structures of the inactive state along aMD simulations trajectories. In particular, the transition of the active to inactive conformation requires 17 ± 4 ns in the receptor empty form, whereas the transition occurs in a longer time in the complexed receptor, that is, 156 ± 7 and 274 ± 8 ns bound to CPB and Sal, respectively. The absence of a signaling protein at the intracellular side of the receptor that stabilizes the conformational changes during activation led to a spontaneous transition of the agonist–receptor complexes from the active to inactive-like state.[16,17] Nevertheless, we hope to capture conformational changes in the receptor bound to CPB that might lead to the active state signaling through β-arrestin. Interestingly, the receptor bound to Sal sampled the displacement defined by the crystal structure of the active state bound to the G protein in cMD and aMD simulations.

Figure 2

Projection of the receptor conformational space simulated by cMD and aMD onto the two lowest principal components produced from the analysis of all simulated trajectories. The projection of the crystal structures is in black. The two black dots with positive values along PC1 correspond to the projections of the active states (the β2AR crystal structures with the G protein-like antibody and the Gs protein).

In our cMD simulations, the active-to-inactive state transition was monitored only in the empty receptor, where the transition time was 41 ns. Dror et al.[25] have shown that the active state of the β2AR bound to the full agonist requires a substantially longer time (>6 μs) to move to the inactive state in the cMD simulations.

Overall, PCA shows enhanced sampling in aMD simulations compared to cMD simulations and relevance of aMD for investigating conformational changes in the empty receptor and bound to different agonists. We thus use aMD trajectories for further analysis of the receptor conformational changes.

aMD Simulations Capture the Inactive-like State Characteristics

To further characterize the aMD simulated conformational ensembles of β2AR in the empty form and bound to biased agonists in a detail, we focus on the appearance or disappearance of interhelical interactions and water clusters in the simulation trajectories, which are distinctive for the crystal structures of the receptor in the inactive and active states. Furthermore, the study of these receptor local changes is a way to validate the performance of aMD, whether aMD simulations reproduce specific contacts, in addition to helix movements captured by the PCA.

Y2195.58 forms a hydrogen bond with Y3267.53 of the highly conserved NPxxY motif in the crystal structures of the active state, holding helix 6 in the outward position.[16,17] In the representative trajectory of our aMD simulations, this hydrogen bond was broken at 20, 32, and 97 ns in the empty receptor and bound to CPB and Sal, respectively, followed by the movement of Y3267.53 to a position similar to one in the inactive state position (Figure 3A). The breakage is accompanied by the movement of helix 6 into the helical bundle. A similar pattern of the hydrogen bond breakage was observed in other simulated replicates.

Figure 3

Evolution of the structural characteristics of the β2AR inactive state: (A) breakage of the Y2195.58–Y3267.53 hydrogen bond monitored in the form of the distance between oxygen atoms of the hydroxyl group and the number of hydrogen bonds; (B) inward movement of the F2826.44 side chain monitored as a distance between the center of mass of the aromatic ring (carbon atoms) and the Cα of I1213.40; (C) formation of the salt bridge between R1313.50 and E2686.30 (known as an ionic lock) monitored as the distance between the CZ atom of R1313.50 and CD atom of E2686.30 and the appearance of the salt bridge. The evolution of distances and hydrogen bonds is shown with solid and dashed lines, respectively. MD analyses for Sal, CPB bound, and receptor in the empty form are shown in black, red, and green, respectively.

The movement of helix 6 in simulations is accompanied by the movement of F2826.44 into the region that this residue occupies in the inactive state. The movement is shown by measuring the distance between the center of mass of the aromatic atoms of F2826.44 and the Cα atoms of I1213.40 over the simulations (Figure 3B). Thus, F2826.44 faces the membrane in the active state and is buried back into the helical bundle by passing I1213.40 to reach the inactive state in our simulations. The movement of F2826.44 to the inactive state requires 44, 112, and 216 ns in the empty receptor and bound to CPB and Sal, respectively.

The key interaction that characterizes the inactive state of the receptor is the salt bridge (known as an “ionic lock”) between R1313.50 of the E(D)RY motif and E2686.30. During the activation this interaction breaks, leading to the outward movement of helix 6.[47] We have observed in our simulations the formation of this interaction within 32, 54, and 290 ns in the empty and CPB- and Sal-bound receptors, respectively. Formation of these interactions locks helices 6 and 3 back to the inactive state (Figure 3C). A movie of the transition from the active to inactive states for the empty receptor is available in the Supporting Information. In cMD simulations, the formation of this ionic lock during the transition from the active to inactive state required ∼6 μs.[25]

We further examined the distribution of water molecules in the helical bundle during the simulations. Because the crystal structures of the active state contain no water molecules in the interhelical cavity, the simulations were initiated with the empty cavity. Notably, we found that few water molecules entered the interhelical cavity from the intracellular side during the equilibration period of the simulations, but none of them formed stable hydrogen bond contacts. However, in the long production period of the simulations we identified three stable water clusters with a significant residence time, which are summarized in Table 2 for the representative trajectories and shown in Figure 4. These clusters resemble the water clusters found in the crystal structure of the inactive state of β2AR and other GPCRs[48] (Supporting Information, Figure 1S).

Table 2

Water Binding Sites during the Active-to-Inactive State Transition of the β2AR Empty Form and Bound to CPB or Sal in aMD Simulationsa

water cluster	water binding site residues	Nwat (X-ray)	Nwat/τo /τr (empty), ns	Nwat/τo/τr (CPB), ns	Nwat/τo/τr (Sal), ns
1	N511.50 D792.50 I2776.39	6	5/51/161	4/101/111	5/28/272
	T2816.43 W2866.48 G3157.42
	N3187.45 S3197.46 N3227.49
2	C2856.47 F2896.51 L3117.38	1	1/3/119	1/43/169	1/130/225
3	A461.45 G501.49 G3267.47	1	1/1/17		1/20/50
					1/301/25

Nwat is the number of water molecules in the X-ray structure of the β2AR inactive state, in the aMD simulations of the β2AR empty form and bound to CPB or Sal. τo is the occupancy time that the water molecule needs to occupy the water binding site. τr is the residence time that the water molecule spends within the water binding site.

Figure 4

Three stable water clusters identified in the aMD simulations. The snapshot was made from the simulations of the empty form. Cluster 1 involves up to five water molecules, which are localized within helices 1, 2, and 7, shown in yellow. Cluster 2 includes one water molecule in the water binding site formed by the backbone of C285, F289, and L311, shown in purple. Cluster 3 contains one water molecule forming the hydrogen bonds with the backbone of A46, G50, and G326 shown in pink. The identified clusters correspond to the water clusters found in the crystal structure of GPCRs in the inactive form.

Several water molecules, mainly from the intracellular side, entered the space between helices 1, 2, and 7 in the agonist-bound receptor and formed hydrogen bond interactions with W286, D79, and N51 of the NPxxY motif (Figure 4) with residence times of >100 ns (cluster 1, in yellow, Table 2). The water molecules filled the interhelical cavity in the empty receptor from both extracellular and intracellular sides. This network of interactions and its rearrangement involving the NPxxY motif are known to be important in the GPCR activation and desensitization processes.[49,50]

Furthermore, we found a stable water molecule in the site formed by the backbone of C2856.47, F2896.51, and L3117.38 (cluster 2, in purple) as shown in Figure 4. The water molecule at this position is found in the crystal structure of the β2AR inactive state, where this water stabilizes the proline kink of the WxPF/Y motif, preventing the outward movement of W2866.48.[48] The G protein active state of the receptor likely does not have this water binding site[16,17] as the helical turn after Pro has a little unwinding such that the distance between the residues in the inactive sate increases, making it difficult to form several hydrogen bonds with a water molecule at the same time to keep the water molecule around these residues.

We observed a stable water molecule in the site formed by the carbonyl group of the backbone of A461.45, G501.49, and G3267.47 (cluster 3, in pink) in the empty form and Sal-bound receptor (Figure 4). The CPB-bound receptor does not have this water site in the simulated replicates. This is because of an increase of distance between helices 1 and 7 (see the last section) that brings A461.45 and G501.49 away from G3267.47; thus, the water molecule cannot be held by these residues anymore.

In addition, in our simulations we have monitored the water molecules in the positions of other structural water clusters (Supporting Information, Figure 1S) found in the crystal structure of the β2AR inactive state, but these water molecules had a short residence time (<1 ns). Interestingly, in our 50 ns simulations of the empty or ligand-bound β2AR inactive state with crystal water molecules, the water molecules in these positions were not stable either.[51]

In summary, our aMD simulations show that the receptor empty form and bound to two agonists converts to the inactive-like state with recovery of known microscopic interactions of the inactive state in the nanosecond time scale within three conducted replicates of simulations, indicating the relevance of aMD simulations in the active-to-inactive state transition.

aMD Simulations Show Stabilization of Distinct Ensembles of Transient Conformations

To examine the receptor conformational space stabilized by the biased agonists and to characterize stable transient conformations in the simulations, we further use PC1 and PC2 to calculate the probability sampling landscape from the probability of occurrence of a transient conformation within the PC1–PC2 space, as shown Figure 5. The plot also shows the projection of the available crystal structures (black dots).

Figure 5

Sampling probabilities from simulations indicate differences in occupancies of three major stable conformational ensembles stabilized by agonists. The sampling probability landscapes are calculated from the probability of occurrence of a transient conformation within the conformational space simulated by aMD. The conformational sampling is shown by the two lowest principal components calculated from the principal component analysis of all the simulated trajectories. The projection of the crystal structures is in black dots.

From the constructed sampling probability plot, we identified one and the most populated transient ensemble of conformations in the simulations of the receptor empty form. This ensemble is in the center of the plot with the conformations being close to the inactive state (the helix Cα rmsd is 1.9 Å) (Table 3). It is characterized by the breakage of the Y2195.58–Y3267.53 hydrogen bond and the intermediate position of F2826.44 (discussed above), which brings helix 6 to the semi-inward state. Further formation of the R1313.50–E2686.30 salt bridge in this ensemble leads to pulling helix 6 to the fully inward state, forming the inactive state.

Table 3

Root-Mean-Square Deviation of the β2AR Helix Cα Atoms of the Stable Ensembles Identified in the aMD Simulations of the Receptor Empty Form and Bound to Sal and CPB

stable ensemble	Cα rmsd from the active state X-ray, Å	Cα rmsd from the inactive state X-ray, Å	ensemble 1 (Sal), Å	ensemble 2 (Sal), Å	ensemble 3 (Sal), Å
1, empty	2.7	1.9	2.2	2.1	2.1
1, CPB	1.6	2.9	2.0	2.5	3.4
2, CPB	2.8	1.5	2.8	1.8	2.3
3, CPB	3.1	1.5	3.2	2.8	2.5
1, Sal	1.4	2.8
2, Sal	1.7	2.5
3, Sal	3.0	1.5

In contrast, the receptor in the presence of the biased agonists has three major ensembles of the stable transient conformations. However, localization and population of these ensembles are different, suggesting that the ligands stabilize distinct ensembles of conformations. Thus, in the receptor conformational space stabilized by CPB, the first ensemble, which is located closer to the active state, is highly populated with a helix Cα rmsd of 1.6 Å (Table 3) from the initial state and is characterized by the broken Y2195.58–Y3267.53 hydrogen bond. The second ensemble, where F2826.44 occupies the inactive position and helix 6 has the inward position, is close to the inactive state (1.5 Å). The third ensemble is characterized by the formation of the ionic lock with an rmsd of 1.5 Å from the inactive structure.

The receptor conformational space stabilized by Sal involves two ensembles with the Cα rmsd close to the active state (1.4 and 1.7 Å, respectively), whereas the last and more populated ensemble is closer to the inactive state (1.5 Å) (Table 3). The first ensemble of the conformations is characterized by the broken Y2195.58–Y3267.53 hydrogen bond, the second ensemble has, in addition, F2826.44 at the position of the inactive state, and the third ensemble has the fully inward position of helix 6.

We have also calculated the Cα rmsd between the stable ensembles of the receptor conformations bound to Sal and CPB (Table 3) and found a notable deviation between the structures, supporting the transient stable ensembles stabilized by the agonists being structurally different.

The porcupine plots calculated for the last stable ensembles of the biosystems using the PC1 show notable variation in directions and amplitudes of harmonic motions of atoms in the region of helices 2, 7, and 8 (Supporting Information, Figure 2S), indicating differences in the dynamic properties of complexes stabilized by the agonists that might reflect different pharmacological properties of CPB and Sal. We therefore further explore the structural changes in the β2AR caused by the agonists.

Effect of the Cyclopentyl Group on the Interhelical Network of Interactions

Whereas the catecholamine groups of CPB and Sal keep similar contacts with the Ser2035.42, S2075.46, and N2936.55 along the simulations (Supporting Information, Figure 3S), the cyclopentyl ring of CPB provides a major difference in binding by forming additional contacts with helices 2, 3, and 7. The involvement of these helices in binding of the biased agonists has been demonstrated in the recent crystal structures of β1AR in a complex with bucindolol or carvedilol.[53]

The immediate effect of the cyclopentyl group within the binding site is measured by the time series of the radius of gyration of the cyclopentyl group binding pocket composed by W1093.28, G902.61, and I3097.36 in Figure 6. The radius of gyration of the pocket is smaller than in the Sal complex or the empty form of the receptor, suggesting that the pocket shrinks as a result of hydrophobic interactions of the cyclopentyl ring with these residues. This pattern in the simulated system has been observed in the other two simulation replicates.

Figure 6

Dynamics of the cyclopentyl group binding pocket in the empty form of β2AR and bound to CPB and Sal. The fluctuation of the size of the pocket is calculated via the radius of gyration created by G902.61, W1093.28, and I3097.36.

To further observe the impact of the cyclopentyl group on the receptor dynamics, we calculated the distortion of helices 2 and 7 in the most populated transient ensembles of the receptor conformations identified in the PC1–PC2 conformational space, in particular, ensemble 1 of the empty state and ensemble 3 of the Sal- and CPB-bound receptor. To perform this task, we computed the angle variation along the helix length as shown in Figure 7. Thus, the distinct angle changes in the helix residues indicate different dynamics of helices in three systems. Helices 2 and 7, colored on the basis of the angle variation, are shown in Figure 7C. Among the ligand-bound receptor states, the notable angle variation is observed at the region of residues 314–321 in helix 7 involving the functional NPxxY motif and the region of residues 74–84 in helix 2. Both regions face each other and are likely involved in ligand-specific side-chain reorganization.

Figure 7

Helix angle profile for helices 2 and 7. The angle is calculated along the length of the helix using the Bendix program. The representative conformations of the last most populated ensembles are shown with enlarged helices 2 and 7 colored on the basis of the angle value.

Table 4 summarizes hydrogen bond occupancy within the network of the NPxxY motif and helix 2 in the simulations, and Figure 8 shows this network in the stable final ensembles of the receptor conformations. In general, we have observed different direct and water-mediated hydrogen bonding involving N511.50, D792.50, D1133.32, S1203.39, Y3167.43, and N3187.45 in three receptor systems. Thus, we noted tightening of interactions between helices 1 and 2 as well as between helices 3 and 7 in the CPB-bound receptor, between helices 6 and 7 in the Sal-bound receptor, and between helices 2 and 7 in the empty form of the receptor (highlighted in bold in Table 4). Although the specific impact of these residues in the stabilization of the β-arrestin conformation of β2AR must be proved experimentally, the importance of the network involving the NPxxY motif has been shown for receptor activation and internalization.[54,55] The importance of residues at positions 3.32 and 7.43 in the stabilization of the CCK2R conformation recruiting β-arrestin has been recently demonstrated in mutagenesis studies.[22] Notably, rearrangement of the NPxxY motif network of interactions changes the position of helices, which can explain the loss of the water binding site between helices 7 and 1 (cluster 3 in Table 2).

Table 4

Hydrogen Bond Occupancy in the Distinct Network of Interactions Involving Helices 1–3, 6, and 7

	occupancy, %
hydrogen bonds	empty	CPB-bound	Sal-bound
N511.50side chain–D792.50side chain	0	81	6
N3227.49side chain–D792.50side chain	82	47	38
Y3267.53side chain–D792.50side chain	63	0	0
N3227.49side chain–S1203.39side chain	0	0	23
W2866.48side chain–N3187.45side chain	15	4	1
Y3087.33side chain–N2936.55side chain	13	51	84
Y3167.43side chain–D1133.32side chain	36	64	42
N3187.45side chain–S1203.39side chain	3	39	0
T2816.43side chain–N3187.45side chain	31	5	5

Figure 8

NPxxY network of interactions in the characteristic receptor conformations of the last stable ensemble of conformations. In addition, we visualized C3277.54 and the salt bridge between K305 and D195, which are known to provide different dynamics in the presence of biased agonists from the 19F NMR and mass spectrometry studies, respectively.

In addition, we noted different dynamics of the extracellular salt bridge involving D192el2 and K3057.32 (Figure 8) in the presence of two distinct ligands; in particular, we found that residues have higher fluctuation in the presence of CPB than in the presence of Sal as the rmsf values of the side-chain heavy atoms are 2.9 ± 0.6 and 1.4 ± 0.3 Å, respectively. As expected, the highest fluctuation of the residues is in the empty form (3.4 ± 0.4 Å). The observed ligand-dependent fluctuation of the residues is consistent with the interpretation of the spectroscopic data, which suggests distinct conformational coupling of the salt bridge in the presence of various ligands.[19,56]

We have also examined the dynamics of C3277.54 (Figure 8) as 19F NMR studies have suggested a different mobility of this residue in the presence of biased and unbiased ligands.[20] We found that C3277.54 has a more solvent-exposed conformation in the receptor bound to Sal compared to the receptor bound to CPB, which correlates with the experimental data.

Our aMD simulations further suggest the different perturbations of helices 2 and 7 in the presence of the biased agonist recruiting β-arrestin compared to the biased agonist activating the G protein.

Discussion

Functional selectivity or ligand-biased signaling in which certain ligands can differentially activate signaling pathways is a new concept in molecular pharmacology developed in the past decade. Although functionally selective ligands and mutations have been described for many receptors, the molecular basis of how changes of functional groups in a ligand affect functional properties of the ligand remains to be elucidated. The recent data from mass and NMR spectroscopies, as well as mutagenesis studies, has started to provide first insights into differences in helix perturbations caused by G protein- and β-arrestin-biased ligands.[19,20,22]

In this work, we have applied the enhanced sampling method of aMD in atomic scale simulations to compare the dynamics of the biased agonists, Sal and CPB, and link them with the tendencies of these ligands to stabilize distinct functional states of the receptor. aMD simulations have been shown to reproduce many structural and dynamics properties of soluble proteins.[27−30] Here, we used aMD simulations in the study of the membrane receptor for the first time and monitored the receptor conformational changes from the active to inactive state on the nanosecond time scale. Because agonists stabilize the active state of the receptor and the phosphorylated active state is recognized by β-arrestin,[41] we have started our simulations with the agonist-bound active state of the receptor using the X-ray structure of the receptor in the complex with the G protein-like antibody,[16] which was available at the initiation of this study. In the absence of the antibody, all complexes were shifted to the inactive-like conformation of the receptor, restoring the known contacts of the inactive state, that is, the inward position of F2826.44, the ionic interaction between R1313.50 and E2686.30, and the water clusters. In addition, we have observed distinct dynamics in the receptor when bound to two different partial agonists, Sal and CPB, notably in the network of interactions involving helix 7. The distinct perturbation of helix 7 in the presence of the β-arrestin-biased ligand compared to the G protein-biased ligand is in good agreement with recent experiential studies.[19,20,22]

A theoretical background of receptor structural plasticity in the form of the energy landscape has been recently illustrated by Deupi and Kobilka.[57] The energy landscape was used in interpretation of the ligand efficacy and the ability of the ligand to trigger a particular signalling mechanism. In addition, the energy landscape was used in understanding the single molecule force spectroscopy experiments, where the effect of ligands with different efficacy on dynamics of β2AR has been investigated.[58] Here, we have explored the sampling probability landscape resembling the energy landscape in the presence of two biased agonists and shown distinct stable ensembles of receptor conformations stabilized by the agonists, which feature different behaviors of helix 7.

Although the specific network of interactions that is responsible for the stabilization of the receptor conformation recruiting the β-arrestin state is yet to be confirmed experimentally, we have attempted to further characterize the different dynamics of two biased ligands that lead to specific receptor conformations triggering distinct signaling pathways. Understanding the network of interactions induced by the biased ligand and the receptor conformational shifts under the biased ligand binding might stimulate the development of more efficient drugs.