Exploring the deactivation mechanism of human β 2 adrenergic receptor by accelerated molecular dynamic simulations.

The β 2 adrenergic receptor (β2AR), one of important members of the G protein coupled receptors (GPCRs), has been suggested as an important target for cardiac and asthma drugs. Two replicas of Gaussian accelerated molecular dynamics (GaMD) simulations are performed to explore the deactivation mechanism of the active β2AR bound by three different substrates, including the agonist (P0G), antagonist (JTZ) and inverse agonist (JRZ). The simulation results indicate that the Gs protein is needed to stabilize the active state of the β2AR. Without the Gs protein, the receptor could transit from the active state toward the inactive state. During the transition process, helix TM6 moves toward TM3 and TM5 in geometric space and TM5 shrinks upwards. The intermediate state is captured during the transition process of the active β2AR toward the inactive one, moreover the changes in hydrophobic interaction networks between helixes TM3, TM5, and TM6 and the formation of a salt bridge between residues Arg3.50 and Glu6.30 drive the transition process. We expect that this finding can provide energetic basis and molecular mechanism for further understanding the function and target roles of the β2AR.

Introduction

The GPCRs represent the largest protein family correlating with signal transduction across membranes (Rasmussen et al., 2011a) and are responsible for the most of cellular responses to hormones and neurotransmitters as well as the senses of sight, olfaction and taste (Isberg et al., 2015). The β 2 adrenergic receptor (β2AR) is one of the best characterized model among the GPCRs because of its three-type identified ligands, namely agonist, antagonist and inverse agonist (Li et al., 2013). Similar to the other members of GPCRs, binding of an agonist leads to a transition of β2AR conformation toward the active state, which further induces association of β2AR with the G-proteins and various intracellular responses (Trincavelli et al., 2010). Differently, binding of an inverse agonist causes a shift of the receptor toward the inactive state, leading to unfavorable binding with the G-proteins thus blocking the signal transduction (Li et al., 2013). The presence of an antagonist blocks substrate binding site and in turn affects the activity of β2AR. The previous studies verified the active or inactive states of β2AR consist of various receptor conformations with different signaling implications (Galandrin and Bouvier 2006; Kobilka and Deupi 2007). Hence it is critical to probe dynamics behavior of conformational changes induced by ligand bindings for understanding functions of β2AR and its target roles in drug design.

Based on significance of β2AR’s dynamics information in drug design toward treatment of ischemic, coronary artery, and inflammatory diseases, many experimental and theoretical works have focused on insights into conformational changes of β2AR (Miao and McCammon 2018; Liu et al., 2019a; Ishchenko et al., 2019). The crystal structure of the β2AR-Gs complex determined by Rasmussen et al. (2011b) suggest that the largest conformational change of β2AR involves an outward motion of 14 Å at the cytoplasmic end of transmembrane segment 6, termed as TM6. Human β2ARs bound by inverse agonists and antagonist exhibit that binding pocket of β2AR can accept compounds of various chemical and pharmacological properties with only minor local conformational changes (Wacker et al., 2010). The crystal structures of the ligand-bound β1AR and β2AR from Xu et al. (2021) uncover that the catecholamine binding pockets are identical between β1AR and β2AR, but the extracellular vestibules have different shapes and electrostatic properties. Liu et al. (2020) solve the structure of the AS408-bound β2AR and their study reveals molecular mechanism of AS408 stabilizing the inactive conformation of the β2AR. The structures of the β2AR adrenergic receptor bound by an orthosteric agonist and a compound 6FA from the work of Liu et al. (2019b) unveil mechanism of the β2AR regulation by an intracellular positive allosteric modulator. These experimental works provide significant structural basis for further investigating the function of the β2AR.

Apart from experimental studies, theoretical computations are also applied to probe the function and target roles of the β2AR. Supervised machine learning and dynamic network analysis from Chen et al. (2021c) show that binding of nanobody produces an allosteric effect on ligand-specific active states and triggers tighter and stronger local communication networks between the Nb80 and the ligand-associated sites. Microsecond-timescale molecular dynamics (MD) simulations on β2ARs in multiple wild-type and mutant states indicate that conformations of inactive β2AR reach an equilibrium between the lock formed state and the lock broken one, whether or not the cocrystallized ligand is present (Dror et al., 2009). Recently, experimental and theoretical insights into conformational changes of β2AR induced by ligand bindings are ongoing (Li et al., 2013; Staus et al., 2016; Liu et al., 2019a; Stanek et al., 2019; Wang and Miao 2019). Despite great success obtained in the previous works, dynamic information relating with the state changes between the active and inactive β2AR is still insufficient. Therefore, it is highly essential to further probe conformational changes of the β2AR for understanding the activity regulation and function of the β2AR.

Conventional molecular dynamics (cMD) simulations (Xue et al., 2018a; Xue et al., 2018b; Sun et al., 2021a) and predictions of binding free energies (Sun et al., 2014; Sun et al., 2021b) can provide useful dynamic information and energetic basis for elucidating functions of targets. Recently, a more efficient sampling technology, Gaussian accelerated molecular dynamics (GaMD) simulation (Miao et al., 2015; Wang et al., 2021), is proposed by Miao and McCammon (2016) to improve conformational sampling of targets, moreover this technology has obtained great successes in exploration of ligand-induced conformational alterations of targets (Wang and Miao 2020; Chen et al., 2021a; Chen et al. 2021b; Chen et al. 2022). More importantly, Draper-Joyce et al. (2021) use GaMD and cMD simulations to probe positive allosteric mechanisms of adenosine A1 receptor-mediated analgesia and their results verify that GaMD simulations better capture conformation transitions between the active and inactive states of adenosine A1 receptor than cMD simulations. In this work, three different ligands, namely agonist P0G (Rasmussen et al., 2011a), antagonist JTZ (Wacker et al., 2010) and (Wacker et al., 2010) are selected to probe deactivation mechanism of the β2AR and decipher free energy profiles that affect conformational transition of the β2AR. The structures of P0G, JTZ and JRZ are displayed in supporting information (Supplementary Figure S1). To obtain rational conformational samplings, Two replicas of GaMD simulations are performed on the P0G-bound active β2AR with the Gs protein (Rasmussen et al., 2011b), P0G-, JTZ- and JRZ-bound active β2AR without the Gs protein. Free energy landscapes (FELs) are built by using reaction coordinates from the GaMD trajectories to reveal energetic basis and clarify dynamics information underlying the deactivation process of the β2AR. This work is expected to provide molecular mechanism for deeply understanding the function of the β2AR.

Methods

System preparation

The initial coordinates of the P0G-bound active β2AR with the Gs protein used for GaMD simulations are taken from protein data bank (PDB) and its entry is 3SN6 (Rasmussen et al., 2011a). The active β2AR removing the Gs protein from 3SN6 is superimposed with the crystal structure 3P0G (Rasmussen et al., 2011b), 3NYA (Wacker et al., 2010) and 3NY8 (Wacker et al., 2010) to ,respectively generate the P0G-, JTZ-, and JRZ-bound active β2AR without the Gs protein. Among four complexes, the active β2AR and the Gs protein together with ligands P0G, JRZ and JTZ are retained at the starting model. All missing hydrogen atoms are bonded to their corresponding heavy atoms by using the Leap module in Amber 20 (Salomon-Ferrer et al., 2013a). The protonation states of residues in β2AR are checked by using the program PROPKA (Li et al., 2005; Bas et al., 2008) and the rational protonation states are assign to the β2AR. The active β2AR is inserted into a palmitoyl-oleoylphosphatidyl-choline (POPC) bilayer with all overlapping lipid molecules removed. The above mentioned complexes are solved at the water box consisting of TIP3P water model (Jorgensen et al., 1983). The amber force field ff19SB (Tian et al., 2020) and amber lipid force field LIPID 14 (Dickson et al., 2014) are used to respectively produce force filed parameters of the Gs protein and β2AR as well as lipid POPC. The atomic partial charges in general AMBER force field (GAFF) and the second generation of GAFF (GAFF2) can generate reliable force field parameters of small molecules and it is used to accurately predict the affinity and binding selectivity of ligands to receptors (Wang et al., 2004; He et al., 2020). Thus, GAFF2 is adopted to yield the force field parameters of ligands P0G, JTZ, and JRZ. The method of the Austin Model 1 with bond charge correction (AM1-BCC) (Jakalian et al., 2000; Jakalian et al., 2002) method is applied to produce the atomic charges of P0G, JTZ, and JRZ through the Antechamber tool in Amber (Wang et al., 2006). The system charges are then neutralized at salt environment of 0.15 M NaCl. All of the aforementioned operations are finished by using the web-sever CHARMM-GUI (Wu et al., 2014; Lee et al., 2019).

GaMD simulations

To remove high-energy contacts between atoms of simulated systems, each system is optimized using the steepest descent minimization of 50,000 steps and the conjugate gradient one of another 50,000 steps. Subsequently, the systems endure a 2-ns soft heating process from 0 to 310 K by restraining heavy atoms of the ligand-bound active β2AR with or without the Gs proteins with 1 kcal/(mol Å2) harmonic constant in a constant number, volume and temperature (NVT) ensemble. Then, four systems are further equilibrated for 2 ns in a constant number, pressure and temperature (NPT) ensemble at 1 atm and 310 K by using the same restraints as in the previous NVT simulation. After that, a 3-ns cMD simulation is performed to examine potential energy statistics, involving the maximum, minimum, average, and standard deviation of four systems. Then, a 30-ns GaMD equilibration with the boost potential is run on the P0G-, JTZ-, and JRZ-bound active β2AR without the Gs proteins, while a 40-ns GaMD equilibration is done on that with the Gs proteins. Finally, two separate 1-μs GaMD simulations are conducted with randomized initial atomic velocities to relax each system. During all current cMD and GaMD simulations, chemical bonds between hydrogen atoms and heavy ones are restrained with the SHAKE algorithm (Ryckaert et al., 1977). The temperatures of the systems are regulated with the Langevin dynamics with a collision frequency of 2.0 ps−1 (Izaguirre et al., 2001). An appropriate cutoff value of 12 Å is adopted to perform calculations of electrostatic interactions with the particle mesh Ewald method (Essmann et al., 1995) and estimations of van der Waals interactions. To be convenient for the post processing analysis, GaMD trajectories of two separate replicas are combined into a single joined trajectory (SJT). The PyReweighting toolkit (Miao et al., 2014) is employed to reweight the data stemming from the CPPTRAJ analysis on the SJT (Roe and Cheatham 2013) and recover the original free energy of four simulated systems. The details of GaMD simulations and principal component analysis (PCA) have been clarified in our previous works (Chen et al., 2021d). All simulations through this current study are run by mean of the program pmemd.cuda inlayed in Amber 20 (Salomon-Ferrer et al., 2013b).

Results

Structural stability of the P0G-bound β2AR with the Gs protein

To understand structural flexibility of the β2AR in the binding environment, root-mean-square fluctuations (RMSFs) of the Cα atoms from four simulated systems are computed by averaging on two replicas (Figure 1A). Compared to the P0G-bound β2AR with the Gs protein, structural flexibility of all domains in the P0G-, JRZ-, and JTZ-bound β2AR without the Gs protein are strengthened, verifying that the Gs protein is needed to stabilize the active state of the β2AR. Except for the loops linking helixes, structural flexibility of helixes TM1, TM2, TM5, and TM6 from the P0G-, JRZ-, and JTZ-bound β2AR without the Gs protein are enhanced relative to that with the Gs protein. Meanwhile, structural flexibility of TM2 and TM4 from the P0G- and JTZ-bound β2AR without the Gs protein is also increased compared to that with the Gs protein. To access the effect of the Gs protein on global flexibility of the β2AR, molecular surface areas (MSAs) of the β2AR are calculated by using the atomic coordinates (Supplementary Figure S2). It is observed that the MSAs of the P0G-, JRZ-, and JTZ-bound β2AR without the Gs protein are increased by ∼702, 1,032, and 684 Å2 relative to the P0G-bound β2AR with the Gs protein, respectively, suggesting that binding of the Gs protein decreases global flexibility of the β2AR.

FIGURE 1

Structural fluctuation and free energy information: (A) RMSFs of the Cα atoms from β2AR in four simulated systems, (B) RMSDs of non-hydrogen atoms from P0G, JRZ, and JTZ in four simulated systems, (C) free energy landscape of the P0G-bound β2AR with the Gs protein constructed using the distances between the Cα atoms of residues Arg3.50 and Glu6.30 and the RMSD of P0G, and (D) the evolution of the distance between the Cα atoms of residues Arg3.50 and Glu6.30 as the simulated time. PMF is scaled in kcal/mol.

To unveil structural fluctuations of ligands P0G, JRZ, and JTZ in binding pocket of the β2AR, root-mean-square deviations (RMSDs) of non-hydrogen atoms from these three ligands are calculated and their frequency distribution are displayed at Figure 1B. Similar fluctuation tendency is observed in the RMSD plot of P0G in the β2AR with/without the Gs protein, which includes two peak values although larger RMSD values are identified in the simulations of system without the Gs protein. Meanwhile, the RMSDs of JRZ and JTZ in the β2AR without the Gs protein are larger than that in the P0G-bound β2AR with the Gs protein. The aforementioned results indicate that the structural fluctuations of P0G, JRZ, and JTZ in the binding pocket are strengthened without binding of the Gs protein, suggesting there are allosteric effects between the agonist and the G protein binding pockets. These results are consistent with the result that the G protein binding can increases the binding affinity of agonist (Rasmussen et al., 2011a).

To reveal the effect of an agonist binding on free energy profile of the active β2AR with the Gs protein, the distance of the Cα atom in Arg3.50 away from that in Glu6.30 and the RMSD of non-hydrogen atoms from P0G are used as reaction coordinates to build the FEL. The distance between the Cα atoms of Arg3.50 and Glu6.30 can reflect conformational transition between the active state and inactive one of the β2AR, while the RMSD of P0G can exhibit its structural fluctuation in the β2AR. The constructed FEL is displayed in Figure 1C. For the P0G-bound active β2AR with the Gs protein, GaMD simulations only capture a low energetic state (Figure 1C) with a RMSD of 0.92 Å relative to the crystal structure 3SN6, agreeing with the X-ray structure (Supplementary Figure S3). The evolution of the distance between the Cα atoms of Arg3.50 and Glu6.30 as the simulation time is calculated (Figure 1D) and the information shows that the β2AR always keeps the active state due to binding of the β2AR, further verifying that the Gs protein is needed for stabilizing the active state of the β2AR.

Free energy landscapes of P0G, JRZ, and JTZ-bound active β2AR without the Gs protein

To investigate the effect of binding of P0G, JRZ, and JTZ on the activity of the β2AR, the same reaction coordinates as that used in the P0G-bound β2AR with the Gs protein are utilized to construct the FELs of the P0G-, JRZ-, and JTZ-bound active β2AR without the Gs protein (Figure 2). The time courses of the distance between the Cα atoms of Arg3.50 and Glu6.30 and structural superimposition of the β2AR from different energetic states are also shown in Figure 2.

FIGURE 2

Free energy landscapes and conformational changes: (A,D,G) corresponding to free energy landscapes of the P0G-, JRZ-, and JTZ-bound β2AR without the Gs protein constructed using the distances between the Cα atoms of residues Arg3.50 and Glu6.30 and the RMSDs of ligands, respectively, (B,E,H) representing structural superimposition of the P0G-, JRZ, and JTZ-bound active, intermediate and inactive β2AR without binding of the Gs protein, individually, and (C,F,I) indicating the evolution of the distance between the Cα atoms of residues Arg3.50 and Glu6.30 as the simulated time in the P0G-, JRZ-, and JTZ-bound β2AR without the Gs protein, separately. PMF is scaled in kcal/mol.

Without the Gs protein, four low energetic states are identified in the P0G-bound receptor simulations, including an active state, an intermediate state and two inactive states (Figure 2A), suggesting that GaMD simulations have successfully captured the deactivation process of the active β2AR. The structural superimposition of the active, intermediate and inactive states reveals the transition pathway of the active β2AR toward the inactive one (Figure 2B; Supplementary Figure S4). Two interesting phenomena are observed: 1) the bottom of helix TM6 goes away from the position of the active β2AR and is close to the helixes TM3 and TM5, 2) the bottom of helix TM5 shrink upwards, especially for the inactive sate (Figure 2B; Supplementary Figure S4). The time course of the distance between the Cα atoms of Arg3.50 and Glu6.30 in two separate replicas is provided in Figure 2C, indicating that GaMD simulations capture the transition of the active β2AR into the inactive state. The time evolution of the distance of the nitrogen atom (N) in Arg3.50 away from the oxygen atom (O) in Glu6.30 in two separate replicas is depicted in Supplementary Figure S5 and the average value of this distance in two simulations Rep1 and Rep2 are ∼5.1 and 10.8 Å, respectively, indicating that a strong salt bridge interaction is formed between residues Arg3.50 and Glu6.50 during the simulation Rep1. The above presented conformational changes of helixes TM5 and TM6 are also observed at the work of Rasmussen et al. (2011b).

In the JRZ-bound active β2AR without the Gs protein, two low energetic states, namely an intermediate state and an inactive one, appear at the FEL (Figure 2D). To clarify the transition process, the intermediate and inactive structures detected by GaMD simulations are superimposed with the initial optimized structure (initial) of the JRZ-bound active β2AR without the Gs protein (Figure 2E; Supplementary Figure S6). The superimposed results show similar phenomena to the P0G-bound active β2AR without the Gs protein: 1) the bottom of helix TM6 leaves the position of the active β2AR and is close to two helixes TM3 and TM5, 2) the bottom of helix TM5 shrinks upwards, which verifies that the transition process of the active β2AR toward the inactive state exists in our current GaMD simulations. To further confirm this transition process, the equilibrated phases in two separate replicas are also used to compute the evolution of the distance between the Cα atoms of Arg3.50 and Glu6.30 as the simulation time and the results are depicted in Figure 2F and Supplementary Figures S7A,B. It is found that the JRZ-bound active β2AR without the Gs protein quickly realizes the transition of the active β2AR toward the inactive state before the first 20 ns of the equilibration in two separate simulations (Supplementary Figure S7B). Meanwhile, time evolution of the distance of the nitrogen atom (N) in Arg3.50 away from the oxygen atom (O) in Glu6.30 in two separate replicas are also estimated (Supplementary Figure S7C) and the average values of this distance in two separate simulations Rep1 and Rep2 are ∼5.1 and 12.0 Å, respectively, indicating that a stable salt bridge interaction is formed through the simulation Rep1.

In the JTZ-bound active β2AR without the Gs protein, GaMD simulations detect four different energetic states, including an active state, an intermediate state and two inactive states, which captures a transition process of the active β2AR toward the inactive state (Figure 2G). The alignment of the β2AR located at the active, intermediate and inactive states display a transition pathway of the active β2AR toward the inactive one (Figure 2H; Supplementary Figure S8). Meanwhile, this structural alignment also reveals two interesting results: 1) the bottom of helix TM6 leaves the active position of the β2AR and moves toward helixes TM3 and TM5, 2) the bottom of TM5 generates an upward shrinkage. The function of the distance between the Cα atoms of Arg3.50 and Glu6.30 as the simulation time is computed with the SJT (Figure 2I) and the information suggests that both the simulations Rep1 and Rep2 capture the transformation of the active β2AR into the inactive state. The time evolution of the distance between the nitrogen atom (N) in Arg3.50 and the oxygen atom (O) in Glu6.30 is exhibited at Supplementary Figure S9. The average values of the Arg3.50-Glu6.30 NO distance in two separate simulations Rep1 and Rep2 are 5.1 and 5.0 Å, respectively, indicating that a stable salt bridge interaction between Arg3.50 and Glu6.30 is produced during the simulations Rep1 and Rep2. The study of Dror et al. (2009) also captures the transition process of the active β2AR toward the inactive state and the formation of the salt bridge between Arg3.50 and Glu6.30, which supports our current results.

As shown in Figures 2B,E,H, bindings of P0G, JRZ, and JTZ all induce the transition of the active β2AR without the Gs protein toward the inactive state. However, the transition process endures different time in the P0G-, JRZ-, and JTZ-bound β2AR without the Gs protein. For the P0G-bound β2AR without the Gs protein, the Rep1 and Rep2 simulations, respectively spend 152 and 339 ns in realizing the transition process and their averaged time is 245.5 ns. As for the JRZ-bound β2AR, the Rep1 and Rep2 simulations individually take 12 and 9 ns to capture the transition process and their average time is 10.5 ns. In the case of the JTZ-bound β2AR without the Gs protein, the Rep1 and Rep2 simulations separately spend 225 and 48 ns in achieving the transition process and their averaged time is 136.5 ns. Thus, the agonist P0G and the inverse agonist JRZ, respectively take the longest time and the shortest time to realize the transformation of the active β2AR into the inactive one, which agrees with the biology process of three ligands (Michel et al., 2020).

To clarify main forces driving the deactivation process of the active β2AR, the interactions of key residues from helixes TM3, TM5, and TM6 in the active, intermediate and inactive states of the β2AR are analyzed (Supplementary Figure S10). In the active β2AR, a cation-π interaction appears between residues Arg3.50 and Tyr7.53, meanwhile residues Ile6.78 and Leu6.37, respectively produce the CH-π interactions with residues Tyr5.58 and Phe5.62 (Supplementary Figure S10A). In the P0G-bound intermediate state of the β2AR, the CH-π interaction networks of Met5.54 with Phe6.44, Met6.41 and Leu6.37 with Tyr5.58 and Leu6.34 with Phe5.62 are detected, separately (Supplementary Figure S10B), while in the P0G-bound inactive state of the β2AR, the CH-π interaction of Leu3.43 with Phe6.44 and Met6.41 with Tyr5.58 together the CH-CH interactions between Val5.61 and Leu6.34 are identified (Supplementary Figure S10C). The CH-π interactions of Leu5.51 and Met5.54 with Phe6.44, Leu6.37 with Tyr5.58 and Phe5.62 and Leu6.34 with Phe5.62 appear at the JRZ-bound intermediate state of the β2AR (Supplementary Figure S10D), but that of Leu3.43 with Phe6.44, Leu6.49 with Phe5.47, Met6.41 with Tyr5.58 and Leu6.34 with Phe5.62 are recognized at the JRZ-bound inactive state of the β2AR (Supplementary Figure S10E). Except for the π-π interaction between His6.31 and Phe5.62 from the JTZ-bound inactive β2AR, the CH-π interaction networks in the JTZ-bound intermediate and inactive states of the β2AR are similar to that in the JRZ-bound intermediate and inactive states of the β2AR (Supplementary Figures S10F,G). By comparison with the active β2AR, the hydrophobic interaction networks in the substrate-bound intermediate and inactive states of the β2AR change, and finally two common interactions between TM3 and TM6, including the CH-π interaction of Leu3.34 with Phe6.44 and a salt bridge between residues Arg3.50 and Glu6.30, are formed in the inactive state of the β2AR, which stabilizes the inactive state of the β2AR, thus the changes in hydrophobic interaction networks among TM3, TM5 and TM6 drive the deactivation process of the active β2AR.

Concerted motions of the β2AR revealed by principal component analysis

The previous FEL analyses reveal that GaMD simulations capture the transition of the active β2AR toward the inactive one, which certainly generates vital effect on dynamics behavior of the β2AR. To check this issue, PCA is performed on the four simulated systems. The first eigenvector arising from PCA is visualized and the results are plotted at Figure 3. It is observed that structure domains in the β2AR exhibit highly concerted motions and binding of the Gs protein affects dynamics behavior of the β2AR. The helixes TM1-TM4 and the loops at the top of the β2AR have an upper right motion and they tend to leave the helixes TM5 and TM6 in the P0G-bound active β2AR with the Gs protein (Figure 3A). However, by comparison with the P0G-bound β2AR with the Gs protein, the helixes TM1-TM4 and the loops at the top of the β2AR have a high concerted motion toward the helixes TM5 and TM6 in the P0G-, JTZ-, and JRZ-bound β2AR without the Gs protein, moreover the helixes TM5 and TM6 have a tendency being close to each other (Figures 3B–D). These different dynamics behavior of structural domains from the β2AR in the states with or without the Gs protein may form a main force to drive the transition of the active β2AR bound by P0G, JRZ, and JTZ toward the inactive state.

FIGURE 3

Concerted motions of structural domains from the β2AR: (A) the P0G-bound active β2AR with the Gs protein, and (B–D) corresponding to the P0G-, JTZ-, and JRZ-bound active β2AR without the Gs protein, respectively.

Discussion

Agonists, antagonists and inverse agonists play significant roles in regulation on the activity of the β2AR. Insights into deactivation mechanism of the active β2AR are of high significance for understanding the function and target roles of the β2AR. To capture the transition process from the active β2AR to the inactive state, we construct four simulated systems by using the active β2AR, including the P0G-bound active β2AR with the Gs protein and the P0G-, JTZ-, and JRZ-bound active β2AR without the Gs protein. 2-μs GaMD simulations are performed to sample conformational space of the β2AR in different binding environment, in which the transition process from the active β2AR to the inactive state is successfully captured.

FELs are constructed by using the distances of the Cα atom in Arg3.50 away from that in Glu6.30 and the RMSDs of three ligands as reaction coordinates to reveal the changes in different energetic states of the β2AR. Our current GaMD simulations identify more low energetic states in the P0G-, JTZ-, and JRZ-bound active β2AR without the Gs protein than in the P0G-bound active β2AR with the Gs protein, indicating that binding of the Gs protein can stabilize the active state of the β2AR. Similar results are found in the RMSF analyses, which suggested that much lower flexibility is identified in the P0G-bound β2AR with Gs protein compared to the P0G-, JRZ-, and JTZ-bound systems without the Gs protein.

Both FELs and the time evolution of the distances between the Cα atoms of Arg3.50 and Glu6.30 suggest that GaMD simulations capture the transition process of the P0G-, JTZ, and JRZ-bound active β2AR toward the inactive state and identify the transition pathway. During the transition process, the helix TM6 leaves the position of the active β2AR and is close to the helixes TM3 and TM6, which is in consistence with the changes in concerted motions of the β2AR uncovered by PCA. The structures of the P0G-bound active, JRZ-bound inactive and JTZ-bound inactive β2AR without the Gs protein captured by GaMD simulations are respectively superimposed with their corresponding crystal structures 3P0G, 3NY8, and 3NYA and their RMSDs are 0.95, 1.34, and 0.94 Å, separately, indicating that our captured energetic structures are in good agreement with the experimental structures. More importantly, our GaMD simulations of the P0G-, JTZ-, and JRZ-bound active β2AR without the Gs protein capture the formation of a salt bridge interaction between residues Arg3.50 and Glu6.30 and the changes in hydrophobic interaction networks among TM3, TM5 and TM6. Thus long-range electrostatic interaction stemming from the salt bridge and the changes of hydrophobic interaction networks play an important role in the deactivation process. This study is expected to provide dynamics information for deeply understanding the function and target roles of the β2AR.