Investigation of allosteric coupling in human β₂-adrenergic receptor in the presence of intracellular loop 3

Protein-membrane system preparation

The system under study was adopted from Ozcan’s work [6] where the initial conformation was the x-ray crystallographic structure of human β₂AR in complex with T4 lysozyme and the partial inverse agonist carazolol (PDB id: 2RH1) at 2.4 Å resolution [2]. The anchor protein T4L was covalently attached to two ends of helices 5 and 6 in order to facilitate crystallization. After removal of T4L, the missing intracellular loop region ICL3 was added between residues 230 and 262 after being modeled as an unstructured loop of 32 residues long via homology modeling tool, MODWEB [16]. Also, the bound ligand at the ligand-binding site was removed and the apo form of the receptor was used as an initial conformation.

The receptor was then embedded in the double-layered 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC) phospholipid cell membrane generated with VMD’s Plug-in tool [17]. The receptor was positioned with an oblique angle of 8° between its main principal component along the membrane and the z-axis [18]. After solvating the protein-lipid system with VMD’s solvate module, Na⁺ and Cl^- ions were added to neutralize the total charge of the system, which is necessary for Particle-Mesh Ewald summation method used in electrostatic energy calculations [19]. The resulting periodic box dimensions were (86x86x100) in Angstrom. The total number of atoms in the system was 68.001, of which 5.055 belong to protein, 20.770 to lipids, 42.135 to water molecules, and 41 to ions.

The system’s equilibration consisted of several stages such as melting of lipid tails, minimization and equilibration with protein restrained, equilibration with protein released and lastly the production runs [6]. Nanoscale Molecular Dynamics (NAMD) v2.8 software tool was used for all our MD runs [20]. The force fields used were CHARMM27 [21, 22] for lipids, CHARMM22 [23, 24] for proteins and TIP3P model for water in the system [25]. In this work, the isothermal-isobaric (NPT) ensemble was employed using Langevin dynamics in order to keep the temperature constant with a Langevin damping coefficient (gamma) of 5/ps for all non-hydrogen atoms. The pressure was kept constant at 1 atm using a Nose-Hoover Langevin piston with 100 fs period and 50 ps damping timescale [26]. Long-range electrostatic interactions were treated by particle mesh Ewald (PME) method, with a grid point density of over 1 Å. A cutoff of 12 Å was used for van der Waals and short-range electrostatics interactions with a switching function. Time step was set to 2 fs by using SHAKE algorithm for bonds involving hydrogens [27] and the data was recorded at every 200 ps. The number of time steps between each full electrostatics evaluation was set to 2. Short-range non-bonded interactions were calculated at every time step.

Independent MD runs with and without restraints

The restraint distances that we selected for restraint MD runs are provided in Table 2. Seven of these distances are between Asp113 on TM3 and residues on TM5 (Ser203, Ser204, Ser207), TM6 (Phe289, Asn293) and TM7 (Asn312). One other is between residues on TM6 and TM7. Mainly, the size of the ligand-binding site is determined by the position of TM5, TM6 and TM7 with respect to the more stationary TM3 as depicted in Fig. 1. The distance ranges previously observed for the active state of the receptor in several different experimental studies [28–33] are also listed in the second column.

Table 2 also lists the same distances in both active and inactive x-ray structures (PDB ids: 3SN6 and 2RH1). For the inactive state, all distances are higher than the experimentally observed ranges, as expected. Surprisingly, for the active state, almost all distance values are slightly out of range. Yet, most of them are smaller than those in the inactive state, especially those between three serine residues on TM5 and Asp113 on TM3. Therefore, we will focus on the first three distances (rows in Table 2), which clearly distinguish the active states from the inactive ones.

In the first restrained MD run (rstr1), Ser203-Asp113 distance was set to 16 Å for the first 300 ns and then increased to 17 Å for the remaining 200 ns. Other distances were set to values closer to those observed in the inactive crystal structure. Here, the value of 17 Å was explicitly selected as it was recorded in the last frame of the original MD run (MD1μs), when ICL3 closed upon G-protein binding site.

In the second restrained MD run (rstr2), the same three critical distances between serine residues and Asp113 were restrained to 8 Å, 10 Å and 8 Å, which all fall within the experimental range of the active state. Similarly, the other distances were set closer to those in the crystal structures. The third run (rstr3) almost share the same set of restraints as in rstr1 and rstr2, except that it starts with a different initial conformation of the receptor (See Table 1). In the fourth run (rstr4), the same initial conformation was used as in rstr3, with an additional restraint distance set between two backbone alpha-carbon atoms in Ser207 and Asp113. Finally, in the last run (rstr5), the final frame of the original MD run (MD1μs) was used as a starting conformation and the ligand-binding site was severely constricted.

Analysis of two continuation runs indicates stability of packed ICL3

Two 500 ns long MD runs were performed as an extension of the original MD1μs simulation (Ozcan et al., [6]), as previously described in Methods section. Initially in MD1μs run, ICL3 was in an extended conformation, and highly mobile as illustrated in the upper portion of Fig. 2c (red ribbon). At around 600 ns of MD1μs, it started to pack under β₂AR and kept its stationary state until the end of the simulation. The aim here was to investigate how long this stationary, and relatively restricted state would carry on. Both extended simulations, the so-called MD1μs_ctd1 and MD1μs_ctd2, selected the final snapshot of the original run as their initial conformation. MD1μs_ctd1 started with the same velocities as in MD1μs’s final state, whereas MD1μs_ctd2 was carried out with a different velocity distribution, in order to enhance the sampling.

During the first continuation run, ICL3 stayed in its close form with only minor fluctuations (~2 Å) as shown in the RMSD profiles illustrated in Fig. 2a. All the RMSD calculations were carried out after the alignment of each MD snapshot in the trajectory to the initial snapshot based on the transmembrane region, as it is the least mobile part of the receptor. The RMSD profiles labeled as All, Core, and Tmemb represent the RMSD of the whole receptor, the core which is the receptor without ICL3 and the transmembrane region which consists of seven alpha helices located inside the membrane, respectively (See Additional file 1: Figure S1). This complete blockage of the G protein’s binding site suggests an inhibition of the receptor’s activity.

In the second continuation run using a different velocity distribution, a temporary increase up to 5 Å in the RMSD value was observed halfway through the trajectory, which led to a slight opening of ICL3 as reflected by the white colored snapshot in Fig. 2b. However, this opening was only temporary and lasted for about 20 ns. This was followed by a sharp decrease in RMSD to 4 Å caused by the closure of ICL3 to a position slightly different than the initial state. ICL3 stayed there for the rest of the simulation.

In both continuation runs, the closed state of ICL3 was preserved, representing an extreme inactive state, where the G-protein binding site was completely blocked. Experimental studies revealed both active and inactive states of the receptor, but none of these structures incorporated the ICL3 region. Here, it is the presence of ICL3 that caused the receptor to adopt such a novel inactive state, which was found to be noticeably stable.

The fluctuation of each residue averaged over the whole trajectory was determined for each simulation and illustrated in Fig. 2d. As expected, in both continuation runs, the mobility of ICL3 stayed at a much lower level than in the original MD simulation. Moreover, a slight decrease of mobility was observed in every part of the receptor, especially on two important loop regions, ICL2 (intracellular) and ECL2 (extracellular), in conjunction with the decrease in ICL3’s mobility.

The stability of ICL3 was further investigated by a detailed analysis of the hydrogen bond network in the loop conformation. Figure 3 illustrates the profile of the residues involved in hydrogen bonding along the trajectory, which was focused on ICL3 and its neighborhood region. By the time ICL3 closure is completed at around 600–700 ns, a total of 8 stable hydrogen bonds has been observed between ICL3 and the rest of the receptor (core region), which was maintained throughout the simulation. It is noteworthy that this stable network of hydrogen bonds was located mostly at the two junctions of ICL3-TM5 and ICL3-TM6. In the first continuation run, nearly all 8 hydrogen bonds were preserved, whereas in the second continuation run, half of them was lost when a slight opening was observed, but still, an alternative close state of ICL3 was adopted later towards the end of the simulation with most of the hydrogen bonds recovered.

In our simulation studies, the closure of ICL3 was strongly coupled with the lower part of TM6, which exhibited an inward movement of 7.5 Å, in the opposite direction of the outward movement of 14 Å observed during activation (experimentally measured at the Cα carbon of Glu 268 [1]). The RMSD profiles of the intracellular part of TM6 with respect to the active state (PDB id:3SN6) illustrated at the top section of Additional file 2: Figure S2 show that as ICL3 started to change its conformation to a closed state, TM6’s intracellular part shifted to the opposite direction of activation and stayed there for both continuation runs. On the other hand, the intracellular part of TM5 attached to ICL3 at the other end, seemed unaffected by these conformational variations. As illustrated at the lower section of Additional file 1: Figure S1, TM5 stabilized at around 2 Å during the original MD as well as both continuation runs.

Two of the key residues at the binding site are Asp113 on TM3 and Ser207 on TM5, which are known to interact both with agonists and antagonists, via hydrogen bonds or close contacts. They are situated at the two distant corners of the binding site and when the ligand is favorably bound, Ser207 is near the ligand’s aromatic moiety, while Asp113 usually makes multiple hydrogen bonds with the ligand’s polar end group (See Fig. 1b). Therefore, the distance between these two residues directly controls the binding capability of the ligand. Experimental measurements already determined an approximate distance range of [8 Å -10 Å] between the two side chain atoms, O^γ of Ser207 and C^γ of Asp113, when the receptor was found in its active state [31, 32]. As the receptor passes from an active state to an inactive one, the same distance also increases and stabilizes roughly at around 11 Å -12 Å. Thus, as an indicator of activation/inactivation, the same distance was monitored for both continuation runs.

In our original MD run, a close correspondence between this value and the conformational state of the lower part of TM6 was established as shown in Fig. 4; as ICL3 exhibited its major conformational shift from an open to a closely packed state, the lower part of TM6 shifted towards the core of the receptor (See Fig. 4d) and at the same time, the Ser207-O^γ and Asp113-C^γ distance started to increase up to 17 Å - 18 Å, which is majorly caused by the outward shift of TM5 (See Fig. 4c). In both continuation runs, the same distance fluctuated within a range of 13 Å - 18 Å, which is still above 11 Å -12 Å of the crystal structure of the inactive state [1, 2] (See Fig. 4a and b). Moreover, no significant conformational change in the lower part of TM6 was observed, which is mainly caused by the stationary ICL3.

Rapid closure of ICL3 is observed as restraints expand the ligand-binding site

The goal here was to investigate the allosteric coupling between the intra- and extracellular parts of the receptor, by applying some distance restraints to several key residues at the extracellular ligand-binding site region (See Fig. 1). As listed in Table 1 (See Methods section), there exist seven distances between side chain atoms that were experimentally observed within a certain range when the receptor adopted an active state [28–33]. In our first constrained simulation (rstr1), one of those distances which exists between S203O^γ and D113C^γ, was restrained to 16 Å for 300 ns and later increased to 17 Å for another 200 ns, while the remaining six were restrained to those observed in the inactive crystal structure (PDB id: 2RH1) for the whole 500 ns long simulation (See Table 2).

The high value of 17 Å was especially selected for S203O^γ-D113C^γ distance in order to reveal the same allosteric response of the intracellular part of TM6 and ICL3 observed previously in the original MD run. As expected, a close correspondence was observed between the extracellular and intracellular parts of the receptor, as ICL3 packed towards the core of the receptor by the end of 200 ns, which is about 400 ns earlier than in the original MD1μs.

The closure of ICL3 was monitored through the x and y coordinates of its center of mass, as illustrated with colored points corresponding to different time regimes in Fig. 5a. The interacting alpha helical part of G protein was shown as a straight line connecting all its x-y coordinates, simply to give an idea about its position with respect to ICL3. Also, in Fig. 5b and c, a total of 20 snapshots gradually changing color from red to white and finally to blue well demonstrate the closure of ICL3 towards the core of the receptor during simulation in different angles.

ICL3 preserves its open conformation as restraints narrow the ligand-binding site

A second restrained MD run (rstr2) was performed with the same initial frame as used in the first run. This time, the ligand-binding site region was narrowed down via bond restraints between three serines (Ser203-O^γ, Ser204-O^γ, Ser207-O^γ) and Asp113-C^γ to 8 Å, 10 Å and 8 Å, respectively. The simulation was performed for a total of 500 ns. ICL3 preserved its initially open conformation throughout the simulation, in agreement with the allosteric coupling behavior between the intra- and extracellular parts. Similar to the first restrained run, the position of ICL3’s center of mass was monitored and all 20 snapshots were illustrated from the side and the intracellular views as in Fig. 6.

ICL3 closure necessitates the outward tilt of TM5

One important finding about the closure of ICL3 in the first restrained run and also the original run was the simultaneous outward tilt of TM5 towards the lipid bilayer, which is crucial in initiating the conformational changes along TM5 and TM6 and consequently on ICL3 (See Fig. 4c). In both runs, the distance restraints applied to residues on TM3 and TM5 shifted TM5 but not the more stationary TM3. Consequently, this desired outward tilt in the extracellular part of TM5 was followed by the inward tilt of the intracellular part of TM5 and also of TM6, which induced the expected ICL3 closure (See Fig. 4d).

The necessity of TM5’s outward tilt was demonstrated in another 500 ns long restrained run (rstr3 in Table 2) that used similar restraints as in the first restrained run, but an alternative initial conformation, in which ICL3 was in an extended form but slightly packed and oriented towards the core of the receptor (See Additional file 3: Figure S3). The applied restraints simply necessitated an expanded ligand-binding site, which was expected to induce the closure of ICL3. However, no closure was observed in ICL3, which covered a wide range of alternative states nearby G-protein binding site and towards the end of 500 ns, ended up close to its initial position (See Additional file 4: Figure S4). When the conformational change in TM5 was observed, it was clear that as a result of the distance restraint, the outward tilt in TM5 was not notable since both TM3 and TM5 moved apart at the extracellular side (See Additional file 5: Figure S5). Furthermore, no major conformational change in the intracellular part of TM5 was observed. Consequently, ICL3’s motion stayed random between open and close states, and no closure was observed. This result shows that the inward tilt of TM6 at the intracellular side was not enough to induce the closure of ICL3, which necessitates the inward tilt in both TM5 and TM6.

Our next attempt in rstr4 was to impose an additional bond restraint that will bring out the desired outward tilt in TM5, which was not obvious in our previous run (See Table 2). Since the backbone atoms’ fluctuations are usually minor compared to those of side chain atoms, the bond restraint of 17 Å was imposed between two backbone atoms, C^α atom of Ser207 and C^α atom Asp113. This time, the new additional restraint was expected to cause the important outward tilt in the extracellular part of TM5. Indeed, both the expected ICL3 closure and the desired outward tilt in TM5 were observed. In addition, ICL3 closure was accomplished under 100 ns, which was two times faster than the first restrained run (See Fig. 7). Another difference was the final position of ICL3, which was shifted about 5 Å in the x-axis with respect to the previously observed positions and located towards the middle of the G-protein binding cavity. In order to further investigate the stability of ICL3 in this alternative closed state, another 500 ns long MD run (MD500ns) was performed with all restraints removed (run #8 in Table). ICL3 preserved its closed state as illustrated with the center of mass profile in Additional file 6: Figure S6.

Packed ICL3 could not be opened by constricting the ligand-binding site

The final restrained run (rstr5) was set up to observe the allosteric effect caused by narrowing the ligand-binding site region. The final snapshot of the original MD run (MD1μs) was taken as the initial conformation. Here, the ICL3 was fully packed, blocking the G-protein binding site. The ligand-binding site was severely constricted with distance values of 8 Å between almost all pairs of atoms (See Table 2), yet any attempt failed to free the ICL3, which only covered a very confined space during 500 ns long MD run (See Fig. 8). This last result simply point to an important aspect of the receptor’s dynamics. It is rather easy to induce the packing of a loose ICL3 by expanding the extracellular binding site region. Yet, it is almost impossible to unpack an already packed or a half packed ICL3 by simply narrowing the extracellular binding site region. Clearly, the energetic barrier to unpack the ICL3 and consequently open the G-protein binding site is too high to be overcome by a few restraints applied at a far region of the receptor. This energetic barrier is most likely due the existence of several hydrogen bonds that exist between ICL3 and the adjacent ends of TM5 and TM6. Thus, the outward tilt of TM6 including the ICL3, which exposes the G-protein binding site, needs to be induced by some exterior forces acting directly on that specific region only.

One activation mechanism proposed by Dror et al. [34] also supports this finding. They have shown that in its basal form, the receptor’s intracellular part of TM6 fluctuated between open or half open (intermediate) states, and adopted a fully open active state only when a G protein was bound from the intracellular region and pushed the binding site to an open form. If an agonist was bound at the extracellular binding site, then this active state was stabilized. On the other hand, when G protein was released from this agonist-bound state, it was observed that the receptor’s intracellular part of TM6 quickly returned to its inactive state obstructing the G-protein binding site. This finding indicated that the active state cannot be induced by some agonists alone and can only be reached by some exterior forces.