Effect of intracellular loop 3 on intrinsic dynamics of human β2-adrenergic receptor
© Ozcan et al.; licensee BioMed Central Ltd. 2013
Received: 3 June 2013
Accepted: 6 November 2013
Published: 9 November 2013
To understand the effect of the long intracellular loop 3 (ICL3) on the intrinsic dynamics of human β2-adrenergic receptor, molecular dynamics (MD) simulations were performed on two different models, both of which were based on the inactive crystal structure in complex with carazolol (after removal of carazolol and T4-lysozyme). In the so-called loop model, the ICL3 region that is missing in available crystal structures was modeled as an unstructured loop of 32-residues length, whereas in the clipped model, the two open ends were covalently bonded to each other. The latter model without ICL3 was taken as a reference, which has also been commonly used in recent computational studies. Each model was embedded into POPC bilayer membrane with explicit water and subjected to a 1 μs molecular dynamics (MD) simulation at 310 K.
After around 600 ns, the loop model started a transition to a “very inactive” conformation, which is characterized by a further movement of the intracellular half of transmembrane helix 6 (TM6) towards the receptor core, and a close packing of ICL3 underneath the membrane completely blocking the G-protein’s binding site. Concurrently, the binding site at the extracellular part of the receptor expanded slightly with the Ser207-Asp113 distance increasing to 18 Å from 11 Å, which was further elaborated by docking studies.
The essential dynamics analysis indicated a strong coupling between the extracellular and intracellular parts of the intact receptor, implicating a functional relevance for allosteric regulation. In contrast, no such transition to the “very inactive” state, nor any structural correlation, was observed in the clipped model without ICL3. Furthermore, elastic network analysis using different conformers for the loop model indicated a consistent picture on the specific ICL3 conformational change being driven by global modes.
KeywordsICL3 Molecular dynamics simulation Transmembrane helix 6 G-protein binding site Ligand docking Essential dynamics
As the largest family of membrane proteins in the human genome, the G protein coupled receptors (GPCRs) are structurally characterized by the presence of seven membrane-spanning α-helical segments with an extracellular N terminus and an intracellular C terminus. Upon binding to agonists, a series of conformational changes propagate along transmembrane helices and reach the intracellular part of the receptor, which directly interacts with the hetero-trimeric G-protein. Consequently, G protein’s activation triggers different cascades of events depending on the type of agonists bound to the receptor. Therefore, as the initiation point to the flow of signals into cells, GPCRs are associated with a plenty of diseases that make members of this family significant pharmacological targets.
The first solved X-ray crystal structure of GPCR belongs to bovine rhodopsin [1, 2], which is followed by the crystal structure of human β2-adrenergic receptor (β2AR) in the inactive state [3, 4]. Since 2007, the cholesterol bound form of β2AR (PDB:3D4S) , the structure of turkey β1-adrenergic receptor (PDB:2VT4) , the structure of a methylated β2AR (PDB:3KJ6)  and various forms of inactive states of β2AR bound to antagonists such as ICI 118,551 and alprenolol (PDB:3NY8,3NY9,3NYA,3PDS) [8, 9] have been reported. Finally, the nanobody-stabilized active state of β2AR in complex with G-protein, has been solved by Rasmussen and his coworkers (PDB:3SN6) [10, 11]. Still, these static pictures of the receptor remain insufficient to describe the dynamic character of the receptor, which governs the function. It is a well-established concept that proteins have an intrinsic ability to sample an ensemble of distinct conformations in order to perform certain functions . The ligand simply selects the optimal receptor conformation for binding followed by an induced fit to stabilize the final conformation. Many questions remain on these multiple, ligand-specific conformational states of β2AR with different levels of activity from fully active to fully inactive, which induce distinct signaling pathways.
The ternary complex model proposed in 1980 by Lefkowitz and his coworkers  describes an allosteric mechanism for receptor activation. The agonist molecule, when bound to the extracellular part, simply promotes and stabilizes the high affinity β2AR-G protein complex. Following the laws of thermodynamics, binding of G-protein increases the receptor’s affinity for agonist binding to the same extent. Fluorescence spectroscopic studies of β2AR by Ghanouni et al.  presented a model with multiple, agonist-specific receptor states, in which the activation occurs through a sequence of conformational changes. They also suggested that the activation barrier for transition from intermediate to active state is high, and that in vivo the barrier is more likely reduced by G protein binding. The presence of an intermediate state is further supported by the fluorescence spectroscopy studies of Swaminath et al. [15, 16], suggesting a mechanistic model for GPCR activation, where agonist binding stabilizes a series of conformational states with distinct cellular functions.
In addition to experiments, several MD simulation studies have been conducted after the inactive and active states of the receptor have been solved by X-ray crystallography. One simulation study by Dror et al.  reveals that the receptor exists between two distinct inactive conformations of the receptor, one with the ionic lock intact and one with the lock broken. In 2011, Dror and his coworkers proposed a completely different activation mechanism in which the structural changes start at the G protein binding site propagating upwards as opposed to agonist-induced conformational changes that start at the agonist binding site and propagate down to G protein binding site . The agonist-bound crystal structure of β2AR without a binding partner (PDB:3PDS) recently revealed by Rosenbaum et al. is found to be identical to the inactive state of the receptor (PDB:2RH1). This suggests that in the absence of a G-protein, the receptor prefers to adopt the inactive conformation whether or not it is bound to an agonist. In other words, the agonist molecule is not sufficient alone to shift the equilibrium to the active state. Dror et al.  also proposed an intermediate state for G-protein binding site, which exists as a part of the receptor’s intrinsic dynamics. Binding of a G-protein to this binding site simply promotes a transition to the active conformation, which is further stabilized by an agonist bound at the extracellular region. The most important feature about the dynamics of β2AR is the strong coupling that exists between the intracellular G-protein binding site and the extracellular ligand-binding site of the receptor [7, 19]. The receptor behaves like a pair of pincers where the intracellular part becomes narrower as the extracellular part becomes wider, and vice versa.
Due to its unstructured nature, ICL3 region is either unresolved in crystallographic experiments or completely removed and replaced by T4-lysozyme (T4L) to facilitate the crystallization. Thus, none of the experimental and simulation studies have discussed the possible effect of ICL3 on the intrinsic dynamics of the receptor. Its replacement by T4L to facilitate crystallization did not prevent agonist-induced conformational changes based on fluorescence spectroscopy measurements . Yet, it is well accepted that its direct interaction with G-protein probably have a significant role on the receptor’s dynamics and the activation/inactivation pathway [21, 22].
In this study, the effect of ICL3 on receptor’s conformational dynamics was investigated via two distinct models of the receptor. Both models were generated from the inactive state of the receptor (PDB:2RH1) after removal of T4L. Moreover, the partial inverse-agonist carazolol was removed from the binding site of both models, since the goal of this work was to provide data about the intrinsic dynamics of the receptor, i.e., the ensemble of conformations accessible to its apo form. According to the current view on ligand binding, the equilibrium distribution of conformational states may be shifted upon ligand binding.
In the so-called loop model, the ICL3 region was modeled as an unstructured loop of 32-residues length and inserted between two open ends of TM5 and TM6 (residues 230 and 263), whereas in the second model, these two open ends were “clipped” or simply covalently attached to each other. The “clipped” model of the receptor, serving as a reference in our study, has been commonly used in recent simulation studies as well [18, 23]. Both models were subjected to 1 μs MD simulation in a POPC membrane bilayer at 310 K. The essential dynamics analysis was carried out to reveal important allosteric coupling within the receptor in the presence of ICL3. Two distinct snapshots taken from the loop model’s trajectory were further used as docking targets for an agonist and an antagonist molecule in order to investigate the effect of ICL3 on binding site conformations. Finally, elastic network analysis was performed on different conformations and loop models to reveal a consistent picture on receptor intrinsic dynamics.
Results and discussion
Presence of ICL3 affects RMSDs and loop mobility
The most significant difference between two time ranges for the loop model is in the mobility of the ICL3 region, which fell down to 2–3 Å in the second time range from 16 Å (out of the range of Figure 3). In addition, a relatively lower decrease in RMSF is observed in all parts of the protein including the loops and the more stable helices in the second time range. However, such a difference cannot be observed in the mobility of the clipped model based on the two time frames (not shown). These results indicate that the fluctuation of ICL3 region in the loop model is directly reflected on every part of the protein structure, including the transmembrane regions. Once ICL3 becomes closely packed under the receptor at around 700 ns (see next section), the mobility of the transmembrane region decreases slightly and becomes more similar to that of the clipped model.
Another important observation about the RMSF profiles in Figure 3 is the mobility of the ECL2 loop region, which is at the extracellular side of the membrane and plays an important role as an access point to the binding site. For the loop model, during the last most stable 300 ns, the RMSF of ECL2 decreases to 5 Å from 6.5 Å, as a consequence of the decrease in the mobility of ICL3. But still, the mobility of ECL2 in the loop model is higher than that in the clipped model irrespective of the time ranges considered. The higher ECL2 mobility allows a wider range of conformational sampling, which would include the open/closed forms of the gateway to the binding site, making the loop model’s binding site more accessible and accommodating for diffusing ligands than the clipped model.
The conformational change of ICL3 gives rise to a “very inactive” state of the receptor
This close packing of ICL3 is quantitatively represented in Figure 4C, which shows the sudden change in the x-offset with respect to the y-offset. The value of x-offset is the difference in the x coordinates between the center of masses for the core of the receptor and the ICL3 region. Similarly, the value of y-offset is calculated from the difference in y coordinates of these center of masses. Three distinct clusters are observed in time ranges of [0–470] ns, [470–670] ns and [670–1000] ns. The first (red) and third clusters (blue) correspond to the open and packed states of ICL3, respectively. The second one (green) represents a transition between the two states. Interestingly, the second cluster’s starting time at 470 ns corresponds approximately to the time at which small helical formations appear inside ICL3 (see Additional file 1: Figure S1). It is an open question whether these helical formations might trigger the transition to the packed state. Additionally, the onset of the third cluster around 670 ns corresponds to the time, at which TM6 starts to deviate from the reference crystal structures (see Figure 4A).
The variation in the distance between Ser207 and Asp113 is illustrated in Figure 5b (top view, looking down from the extracellular side), which shows two conformations of the receptor, with minimum and maximum distance values of 8.3 Å and 18.6 Å (shown as light blue and purple, respectively). Clearly, the intracellular part of TM5 is slightly moving into the core region (see Figure 5b), while its extracellular part is moving away from the core region of the receptor. As a result, Ser207, which is located at the extracellular part of TM5, drifts away from Asp113 on TM3, position of which does not change notably. Similar motions are observed for TM4 and TM6 as well. As a result, the enlarged binding site becomes unfavorable for agonist binding due to lack of some key interactions.
Figure 5c illustrates the correlation between the TM6 shift at the intracellular part of the receptor and the Ser207-Asp113 distance at the extracellular part of the receptor. In the loop model, the change in the RMSD value of TM6 with respect to the active state happens at around the same time as the increase in the Ser207-Asp113 distance, whereas in the clipped model, no such correlation is observed (see Figure 5d). Three different states of the structure are observed in the loop model at around [0–640], [640–770] and [770–1000] ns intervals. In contrast, the conformational variations of the clipped model remain in a restricted area, which corresponds to the first conformational state ([0–660] ns) of the loop model.
The loop model shows a conformational variation in its second intracellular loop (ICL2), which is correlated with the motion of ICL3. Additional file 2: Figure S2 shows four different stages of conformational variation of ICL2 changes with respect to intracellular part of TM6. Both RMSD values are calculated with reference to the active state (PDB:3SN6). However, no such major structural change is observed in the clipped model. The motion of ICL2 in the loop model is also illustrated, in which the ICL2 between TM3 and TM4 steps aside as the ICL3 comes closer to the middle region (initial stage: blue, second stage (at 700 ns): cyan, final stage: red, and active crystal: green). No such conformational rearrangement in ICL2 is observed in the clipped model. Additional file 3: Figure S3 illustrates the change in the RMSD value of ECL2 with the change in the Ser207-Asp113 distance. As the distance increases in the loop model, there is a conformational variation in ECL2 with respect to the active state. However, no such correlation is observed in the clipped model. Clearly, the structural variation in ECL2 is directly affected by the change in the distance as a result of a shift of TM5 away from the binding site, which is in turn a consequence of the ICL3 motion and TM6 shift at the intracellular part of the receptor.
Ionic lock (Arg131-Glu268) is not a molecular switch
During 1 μs long MD simulation of the loop and clipped models, which represent the inactive state of the receptor, the ionic lock profiles are monitored as shown in Additional file 4: Figure S4. Consistent with previous work , the ionic lock seems to be on and off during the simulations of both models (in upper panels). Thus, this ionic lock cannot distinguish between active and inactive states. The cause behind the breakage/formation of this ionic lock is found to be the result of a change in the rotational state of the Χ angle of Glu268, which coordinates perfectly well with the ionic distance profiles (in lower panels).
Furthermore, the increase in the distance between two side chains that form the ionic lock, namely Arg131-N and Glu268-O, coincides properly with the increase in the distance between their alpha-carbons. In the profiles of the inactive state, the backbone distance fluctuates at around 9.5 Å and reaches 12.4 Å at most. However, in the known crystal structure of the active receptor (PDB:3SN6), the distance between alpha-carbons is around 16 Å as a result of a significant outward shift in the intracellular part of TM6. Thus this backbone distance could be one possible measurement for detection of activation.
Essential dynamics analysis reveals the transition to the “very inactive” state in the first principal mode
Additionally, the distance profiles between Ser207 and Asp113 in the loop model was recalculated using Cα atoms only as shown in Figure 7b. A high correspondence between the original and the two reconstructed profiles is observed as in Figure 7a. Also, a plot of RMSD value versus the Ser207-Asp113 distance clearly shows that the essential modes (first and cumulative five) describe the distribution in the original trajectory (Figure 7c). Thus, the closure of the ICL3 driven by the first mode (Figure 7d) is strongly coupled with the opening of the binding site indicated by the Ser207-Asp113 distance. For the clipped model, same profiles are plotted in Additional file 5: Figure S5. The profiles obtained from the projection of the first and cumulative five modes do not explain satisfactorily the dynamics of the extracellular and intracellular regions of the receptor.
Elastic network modeling reveals coupling between global modes and ICL3 conformational transition
MD simulations were performed on a receptor model including a specific unstructured conformation of ICL3 obtained from MODELLER. In order to show the independence of the observed phenomena from the initial MD structure, ANM was performed on four distinct conformations of the loop model. These were selected as the initial, average and final structures of the 1 μs-long loop trajectory, and a receptor model containing an alternative unstructured conformation of ICL3, also provided by MODELLER. The RMSD between the alternative loop model and the one used in our MD simulations is ~ 20 Å for the ICL3 residues after an alignment of transmembrane regions.
where p i and u j represent the ith and jth normalized eigenvectors from PCA and ANM, respectively. The squared inner dot products are generally summed over the first k =10 modes, which describe the collective subspaces of each method. The average overlap values are 0.64, 0.72, 0.62 and 0.66 between the first 10 modes of PCA and ANM performed with the initial, average and final structures of MD run and the alternative loop model, respectively. These values are quite high (relatively closer to 1), representing satisfactory overlap between ANM and PCA subspaces.
The overlap matrices calculated based on residue displacements of ICL3 and TM6 region only (see Additional file 6: Figure S6) indicate several ANM modes exhibiting high overlap with the first mode of MD. Thus, slow modes of ANM clearly drive the significant conformational change of ICL3 and TM6 towards the receptor core, independent of the ICL3 conformation/model used. In summary, our ANM analysis justifies that the ICL3 dynamics observed in MD run can be attributed to be a feature of intrinsic receptor dynamics in conformity with a recent study carried on catalytic loop motions for different enzymes .
The clustering of MD snapshots reveals more conformational variations in the loop model
The clustering profile of the transmembrane region shows four distinct clusters for 1 μs simulation of the loop model as illustrated in Figure 9a. Two of those clusters dominate over the other two, since they contain 54% and 37% of the total snapshots, which are observed at around [0–600] ns and [600–1000] ns, respectively. On the other hand, two distinct clusters are obtained for the clipped model, and only one of them dominates for 92% of the time. For each of the three short MD simulations, there is only one single cluster that dominates during 100 ns. These results indicate that the transmembrane region of the loop model alternates between two distinct conformations, while the clipped model’s transmembrane region prefers to adopt only one. Interestingly, the second conformation in the loop model observed between 600 and 1000 ns coincides with the time at which the ICL3 changes its conformation and the receptor adopts a “very inactive” state with an expanded binding site. A highly similar clustering profile is obtained for the binding site region as shown in Additional file 7: Figure S7. This is an expected outcome considering that the binding site is embedded in the transmembrane region.
Figure 9b shows the cluster profile of the intracellular part of the receptor, which consists of residues interacting with the G-protein based on the active crystal structure (PDB:3SN6) . Clearly, the loop model’s intracellular part samples three distinct states while the clipped model’s intracellular part only samples one conformation. In the three short simulations, there is also one single conformation dominating the others. The structural flexibility of the intracellular part is critical in making contact with the G protein. For the loop model, three snapshots were selected from each cluster shown in Figure 9b as representatives and illustrated in Additional file 8: Figure S8 with a bottom view to show the contact of the receptor with the helical segment of gamma subunit of G protein. In the active crystal structure (PDB:3SN6) taken as a reference and placed on top of the figure, G protein’s helical segment nicely fits the binding cavity. At the initial stages of the simulation, the binding cavity is almost preserved. Towards the end of the simulation, the motion of ICL3 closes down the G-protein binding site almost completely, as shown in the last frame leaving no contact point for the G protein.
The clustering profile in Figure 9c, shows two dominant clusters for the ICL2 region in both model. The two distinct states in the loop model are sampled for about 29% and 58% of the time, while in the clipped model two major clusters are sampled for 39% and 59% of the time. Each of the short simulations of the loop model still does not show structural variation during 100 ns, similar to other two cases above. These results indicate that ICL3 has no significant effect on the conformational sampling of ICL2.
Finally, the cluster profile in Figure 9d shows four distinct clusters for ECL2 region of the loop model, with each consisting of a considerable amount of snapshots and sampled consecutively during the simulation. On the other hand, the clipped model’s one μs simulation as well as three short simulations of the loop model, impart no conformational variation to the ECL2 region. The ECL2 loop region is the second extracellular loop covering the top of the receptor and plays a critical role of providing a passage to the binding site region. Therefore, the ability of ECL2 to sample various conformations, being a functionally important feature for the receptor, is clearly enhanced in the presence of ICL3.
Docking results of epinephrine and ICI to an open and a closed form of the binding site
Two frames of the loop model are selected from the 1 μs trajectory to represent the two extreme cases of Ser207-Asp113 distance value (see Figure 5a, b). One of the conformers is a closed form with a distance value of 8.31 Å, which is in the range of active states (8–10 Å) , while the other conformer represents an open form with the maximum distance value of 18.63 Å. The docked ligands are a natural agonist epinephrine, and ICI, which is an antagonist with a known crystal conformation (PDB:3NY8). The epinephrine is chosen due to its relatively small size, and ICI is selected because it is a large antagonist with an experimentally determined conformation.
The docking results of ICI indicate a complimentary situation. Due to its large size, the antagonist ICI cannot fit into a narrow binding site in the closed form, but it can be favorably docked into a wider binding site, such as the open form. Figure 10C shows the highest score conformations docked in the open (magenta) and closed (light blue) form of the receptor (see Figure 10d for a comparison to the bound state of ICI in the crystal structure shown as black sticks). The RMSD values of docked ICI to open and closed forms with respect to its native state are determined as 3.95 Å and 8.16 Å, respectively. Clearly, ICI when bound into a wider binding site is able to interact with experimentally reported key residues for antagonist binding. ICI interacts with eight residues within a radius of 3.5 Å, Asp113, Tyr199, Ser203, Ser204, Phe208, Trp286, Tyr308 and Asn312, of which five (Asp113, Ser203, Ser204, Trp286, Asn312) are reported as key residues  (see Additional file 9: Figure S9 for the specific interactions of ICI). Furthermore, ICI is correctly oriented in the binding pocket with the hydrophobic catechol ring interacting with Ser203 and Ser204, and its polar end interacting with Asp113. However, in the closed form, ICI is improperly positioned in the binding pocket, interacting within a radius of 3.5 Å with Thr110, Asp113, Val114, Phe193, Tyr199, Phe289, Asn293, Lys305 and Tyr308 of which only two (Asp113 and Val114) are known key residues.
To sum up, the agonist is favorably bound to the closed form, which coincides with the ICL3-open conformation during the first half of the simulation. The ICL3-open conformation may correspond to an intermediate state that promotes G-protein binding, which seems to be stabilized by the presence of the agonist at the binding site. An MD simulation, where the binding site would be constrained to the agonist bound geometry, would give more insight in this perspective. In the second half of the simulation, the open geometry of the binding site is observed to which the antagonist favorably docks. This “very inactive” state corresponds to the closely packed ICL3 that completely blocks the G-protein binding site.
The crystal structure of β2AR has been resolved for the first time in 2007, and since then in silico studies have been conducted to unravel structure-dynamics-function relationship of this G-protein coupled receptor. However, the missing intracellular loop ICL3, which is known to interact with the G protein at the cytoplasmic side, has not been considered or elaborated so far in these studies. In this respect, our MD study exposed the marked effect of ICL3 on collective dynamics and justified the correlated motion between the intracellular G-protein binding site and the extracellular ligand-binding site of the receptor [7, 19]. The receptor behaves like a pair of pincers where the intracellular part becomes narrower as the extracellular part becomes wider, and vice versa.
In accordance with this coupling behavior, our μs long MD simulation of β2AR, which included the modeled intracellular loop ICL3, revealed a so-called “very inactive” state of the receptor, which has not been reported before. In the second half of the simulation, the ICL3 moved toward the core of the receptor and completely blocked the G-protein binding site. Consequently, the intracellular part of TM6, adjacent to ICL3, also shifted toward the core of the receptor. This conformational change in TM6 is in the opposite direction with respect to experimentally observed deformation during activation , which is observed as an expansion or outward movement towards the lipid membrane. This new inactive state of the receptor may provide insight into the design of novel therapeutic drugs.
Another important observation is the correlated motion between the binding site and the G-protein binding site regions. At around the same time when the ICL3 blocked the G-protein binding site, the extracellular binding site of the receptor expanded. The expansion was detected based on the distance profile between two anchor residues, Ser207 and Asp113, critical for agonist binding. In line with these findings, our docking studies indicated favorable antagonist binding to the expanded binding site (closed ICL3) and agonist binding to the closed binding site (open ICL3).
This coupled conformational change seems to be transmitted from the intracellular part to the extracellular part of the receptor via TM5 and TM6. As Sara Linse stated in her review , “a 7TM receptor is like a bundle of rods immersed in the membrane and if a ligand grips the bundle at one end, the bundle opens up like a bouquet of roses at the other end”. In our study, the intracellular end is gripped or held tight by ICL3, and consequently the extracellular part opened up. In other words, ICL3 played the dominant role in inducing the change in the intracellular part, which induced in turn the extracellular part. This dominancy of ICL3 is expected due to its high mobility, which is also a desirable quality for initiation of interactions with intracellular proteins [20, 21]. Thus, we suggest that when left without ICL3, the receptor would not be able to sample that inactive state at all.
This transition to the “very inactive” state took place within a time frame of about 0.1 μs (starting at ~ 0.6 μs). In the last 0.3 μs of the simulation, ICL3, which was observed to be the most mobile region of the receptor during the whole simulation, preserved its close state. Principal component analysis of 1 μs long MD trajectory showed that the first principal mode, which explains 69% of the overall motion governs the transition from the initial inactive state to the “very inactive” state.
At this point it may be argued that our simulation conditions, such as the absence of any ligand at the binding site and/or intracellular proteins that may interact with ICL3, depict a non-physiological environment. Even though, the receptor function is clearly linked with its interacting partners, our aim was to elucidate the intrinsic conformational dynamics of the intact receptor. Based on the widely accepted population shift mechanism [35–37], we tried to uncover the pre-existing conformational states of the apo receptor, which may be shifted and/or modified by the presence of binding partner(s). In fact, our ANM analysis using different conformers/models of ICL3 strengthened our MD results on receptor dynamics. Either the first or second collective mode in ANM was found to be coupled with the specific motion of ICL3, independent of the model used. In contrast, none of these conformational transitions, nor any allosteric coupling between intra- and extracellular parts, were observed in the clipped model simulation lacking ICL3 region. Thus, we stress that the presence of ICL3 provides a more realistic constriction than those of clipped and non-clipped (loose ends) models so far used in β2AR simulations.
Future works will be focusing on the loop model, which will consist of constraining the binding site region to observe the reverse transition (or release) from the “very inactive” state to the inactive/intermediate state and possibly the active state with an opening of the G-protein binding site.
Preparation of the receptor models
The X-ray crystallographic structure of human β2AR in complex with T4 lysozyme (T4L) (PDB:2RH1) at 2.40 Å resolution  was used as the initial conformation. After removal of T4L, the missing intracellular loop region (ICL3) between residues 231 and 262 was modeled as an unstructured loop of 32-residues length via MODELLER . Even though ICL3 possibly exists as an unstructured loop , our modeled loop in Figure 2a interestingly resembles a conformer generated from a fully extended chain by Dror et al. . Moreover, it is expected that ICL3 can sample various conformations during our 1 μs long run.
In the second model, Leu230 and Lys263, which are the two ends of the missing region, were covalently attached to each other via a peptide bond to form the “clipped” model (Figure 2b), which has been commonly used in simulation studies.
Preparation of the environment with the receptor
Each model was later embedded into a palmitoyloleoyl-phosphatidylcholine (POPC) membrane bilayer along the z-axis using VMD’s Membrane Plug-in Tool, v1.1 . The receptor was positioned with an oblique angle of 8° between its main principal component along the membrane and the z-axis . A total of fifteen internal water molecules detected experimentally in the crystal structure were retained because they make hydrogen bonds with the most conserved residues of the receptor and thus possibly contribute to its structural stability. Using VMD’s solvate module, the protein-lipid system was solvated in both intracellular and extracellular sides with a thickness of 15 Å and 13 Å for the “loop” and “clipped” model, respectively. Finally, the protein-lipid-water system was ionized with Na+ and Cl- ions to make the total charge of the system to be equal to zero, which is necessary for Particle-Mesh Ewald summation method used in electrostatic energy calculations. The resulting periodic box dimensions were (86×86×100) and (77×69×90) in Angstrom for the “loop” and “clipped” models, respectively.
Molecular dynamics simulations
Using the suggested procedure for membrane protein system preparation , both models were subjected to three preparation stages. The system consisted of three components of different types, each having a different response time to outside forces. Thus, to reach the equilibrium fast, it was practical to keep some components fixed, while other components were free to move. The first equilibrium stage consisted of melting the lipid tails where only the lipids were free to relax while the protein and waters were held fixed. At the end of the simulation, the unrealistically aligned lipid molecules, transformed into a more disordered, liquid-like structure. In the second preparation step, the protein’s motion was constrained while lipid and water molecules were free to move. Finally in the third stage, the protein was released, and all components were allowed to relax. In every preparation stage, the system was subjected to 1000 steps of energy minimization followed by 0.5 ns MD simulation. At the end of the third stage, the area per lipid was stabilized at 0.635 nm2/lipid in agreement with the experimentally measured value of 0.65 nm2/lipid . Also, the surface area of the membrane in -xy directions decreased due to close packing of lipid molecules with the protein.
Each model was later subjected to 1 μs MD simulation with NAMD v2.7 software tool . CHARMM22 [43, 44] and CHARMM27 [45, 46] forcefields were used to describe the interaction potential of the protein and the lipid respectively, and waters were treated explicitly using TIP3P model . The system was composed of a total of 68,001 and 42,701 atoms for the “loop” and “clipped” model, respectively. MD simulation was performed at constant NPT at 310 K using Langevin dynamics for all non-hydrogen atoms, with a Langevin damping coefficient of 5 ps-1. The system was kept at a constant pressure of 1 atm by using a Nose–Hoover Langevin piston  with a period of 100 fs and damping timescale of 50 ps. 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  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.
For the “loop model”, three additional 100 ns long MD simulations starting with different initial velocities were performed alongside 1 μs MD simulation. The aim was to possibly explore different conformational subspaces than those visited during the long trajectory of the loop model.
Docking was performed using the software tool AutoDock v4.0 . The docking site was selected based on the location of the partial inverse agonist carazolol in the complex structure (PDB:2RH1). Two distinct snapshots taken from MD trajectory were used as target structures. Lamarckian genetic algorithm was used to explore the conformational space. A total of 100 runs were performed for each structure with each run consisting of 1.0×106 and 1.5×106 energy evaluations for epinephrine and ICI ligands, respectively. Grid box constructed with a spacing of 0.375 Å had dimensions of 24 Å × 24 Å × 24 Å for all dockings. For each docking experiment, the pose with the highest score (lowest binding energy of AutoDock) was used as the most probable solution for that complex.
Elastic network analysis
The collective/global modes of the protein were extracted via the anisotropic network model (ANM) [51, 52], which describes the protein as a coarse-grained elastic network of harmonic springs based on a minimum-energy folded conformation. The network is formed simply, by connecting the close-neighboring alpha-carbon atom (called nodes) pairs in the folded structure. The slow or the low-frequency modes extracted from normal mode analysis of the elastic network are known to successfully describe the functional conformational changes.
In our current work, the receptor’s loop model was embedded into a coarse-grained membrane environment according to the methodology provided in Lezon et al.  The membrane consisting of spheres arranged in an FCC lattice had a diameter of 80 Å and a thickness of 33 Å. The cutoff value for pairwise interactions between nodes was taken as 11 Å. The force constants of harmonic springs were selected as 1.0, 2.0 and 4.0 for protein-protein, protein-membrane and membrane-membrane type of pairwise interactions, respectively. In this model, the membrane environment serves as a constriction and thereby inhibits the unrealistically large fluctuations of the transmembrane helices that would be observed if ANM were applied to the protein alone .
OO have carried out the molecular dynamics simulations of loop and clipped models. AU performed the anisotropic network analysis. EDA performed the docking calculations. All authors have participated in the analysis and interpretation of MD trajectory, and writing of the manuscript. All authors read and approved the final manuscript.
G protein coupled receptors
Root mean square deviation
Root mean square fluctuation
Principal component analysis
Anisotropic network model.
This work has been partially supported by The Scientific and Technological Research Council of Turkey (TÜBİTAK, Project # 109 M281), State Planning Organization of Turkey (DPT) (Project # 2009 K120520) and Kadir Has University BAP (Project # 2010-BAP-04). PD acknowledges partial support by Bogazici University BAP Project (No:5714).
- Palczewski K, Kumasaka T, Hori T, Behnke CA, Motoshima H, Fox BA, Trong IL, Teller DC, Okada T, Stenkamp RE, Yamamoto M, Miyano M: Crystal structure of rhodopsin: a G protein-coupled receptor. Science 2000, 289: 739–745. 10.1126/science.289.5480.739View ArticlePubMedGoogle Scholar
- Teller DC, Okada T, Behnke CA, Palczewski K, Stenkamp RE: Advances in determination of a high-resolution three-dimensional structure of rhodopsin, a model of G-protein-coupled receptors (GPCRs). Biochemistry 2001, 40: 7761–7772. 10.1021/bi0155091PubMed CentralView ArticlePubMedGoogle Scholar
- Cherezov V, Rosenbaum DM, Hanson MA, Rasmussen SG, Thian FS, Kobilka TS, Choi HJ, Kuhn P, Weis WI, Kobilka BK, Steven RC: High-resolution crystal structure of an engineered human beta2-adrenergic G protein-coupled receptor. Science 2007, 318: 1258–1265. 10.1126/science.1150577PubMed CentralView ArticlePubMedGoogle Scholar
- Rasmussen SG, Choi HJ, Rosenbaum DM, Kobilka TS, Thian FS, Edwards PC, Burghammer M, Ratnala VR, Sanishvili R, Fischetti RF, Schertler GF, Weis WI, Kobilka BK: Crystal structure of the human beta2 adrenergic G-protein-coupled receptor. Nature 2007, 450: 383–387. 10.1038/nature06325View ArticlePubMedGoogle Scholar
- Hanson MA, Cherezov V, Griffith MT, Roth CB, Jaakola VP, Chien EY, Velasquez J, Kuhn P, Stevens RC: A specific cholesterol binding site is established by the 2.8 A structure of the human beta2-adrenergic receptor. Structure 2008, 16: 897–905. 10.1016/j.str.2008.05.001PubMed CentralView ArticlePubMedGoogle Scholar
- Warne T, Serrano-Vega MJ, Baker JG, Moukhametzianov R, Edwards PC, Henderson R, Leslie AG, Tate CG, Schertler GF: Structure of a beta1-adrenergic G-protein-coupled receptor. Nature 2008, 454: 486–491. 10.1038/nature07101PubMed CentralView ArticlePubMedGoogle Scholar
- Bokoch MP, Zou Y, Rasmussen SG, Liu CW, Nygaard R, Rosenbaum DM, Fung JJ, Choi H, Thian FS, Kobilka TS, Puglisi JD, Weis WI, Pardo L, Prosser RS, Mueller L, Kobilka BK: Ligand-specific regulation of the extracellular surface of a G-protein-coupled receptor. Nature 2010, 463: 108–112. 10.1038/nature08650PubMed CentralView ArticlePubMedGoogle Scholar
- Wacker D, Fenalti G, Brown MA, Katritch V, Abagyan R, Cherezov V, Stevens RC: Conserved binding mode of human beta2 adrenergic receptor inverse agonists and antagonist revealed by X-ray crystallography. J Am Chem Soc 2010, 132: 11443–11445. 10.1021/ja105108qPubMed CentralView ArticlePubMedGoogle Scholar
- Rosenbaum DM, Zhang C, Lyons JA, Holl R, Aragao D, Arlow DH, Rasmussen SG, Choi HJ, Devree BT, Sunahara RK, Chae PS, Gellman SH, Dror RO, Shaw DE, Weis WI, Caffrey M, Gmeiner P, Kobilka BK: Structure and function of an irreversible agonist-beta(2) adrenoceptor complex. Nature 2011, 469: 236–240. 10.1038/nature09665PubMed CentralView ArticlePubMedGoogle Scholar
- Rasmussen SG, DeVree BT, Zou Y, Kruse AC, Chung KY, Kobilka TS, Thian FS, Chae PS, Pardon E, Calinski D, Mathiesen JM, Shah ST, Lyons JA, Caffrey M, Gellman SH, Steyaert J, Skiniotis G, Weis WI, Sunahara RK, Kobilka BK: Crystal structure of the beta2 adrenergic receptor-Gs protein complex. Nature 2011, 477: 549–555. 10.1038/nature10361PubMed CentralView ArticlePubMedGoogle Scholar
- Rasmussen SG, Choi HJ, Fung JJ, Pardon E, Casarosa P, Chae PS, Devree BT, Rosenbaum DM, Thian FS, Kobilka TS, Schnapp A, Konetzki I, Sunahara RK, Gellman SH, Pautsch A, Steyaert J, Weis WI, Kobilka BK: Structure of a nanobody-stabilized active state of the beta(2) adrenoceptor. Nature 2011, 469: 175–180. 10.1038/nature09648PubMed CentralView ArticlePubMedGoogle Scholar
- Bahar I, Chennubhotla C, Tobi D: Intrinsic dynamics of enzymes in the unbound state and relation to allosteric regulation. Curr Opin Struct Biol 2007, 17: 633–640. 10.1016/j.sbi.2007.09.011PubMed CentralView ArticlePubMedGoogle Scholar
- De Lean A, Stadel JM, Lefkowitz R: A ternary complex model explains the agonist-specific binding properties of the adenylate cyclase-coupled beta-adrenergic receptor. J Biol Chem 1980, 255: 7108–7117.PubMedGoogle Scholar
- Ghanouni P, Gryczynski Z, Steenhuis JJ, Lee TW, Farrens DL, Lakowicz JR, Kobilka BK: Functionally different agonists induce distinct conformations in the G protein coupling domain of the beta 2 adrenergic receptor. J Biol Chem 2001, 276: 24433–24436. 10.1074/jbc.C100162200View ArticlePubMedGoogle Scholar
- Swaminath G, Xiang Y, Lee TW, Steenhuis J, Parnot C, Kobilka BK: Sequential binding of agonists to the beta2 adrenoceptor. Kinetic evidence for intermediate conformational states. J Biol Chem 2004, 279: 686–691.View ArticlePubMedGoogle Scholar
- Swaminath G, Deupi X, Lee TW, Zhu W, Thian FS, Kobilka TS, Kobilka B: Probing the beta2 adrenoceptor binding site with catechol reveals differences in binding and activation by agonists and partial agonists. J Biol Chem 2005, 280: 22165–22171. 10.1074/jbc.M502352200View ArticlePubMedGoogle Scholar
- Dror RO, Arlow DH, Borhani DW, Jensen MO, Piana S, Shaw DE: Identification of two distinct inactive conformations of the beta2-adrenergic receptor reconciles structural and biochemical observations. Proc Natl Acad Sci USA 2009, 106: 4689–4694. 10.1073/pnas.0811065106PubMed CentralView ArticlePubMedGoogle Scholar
- Dror RO, Arlow DH, Maragakis P, Mildorf TJ, Pan AC, Xu H, Borhani DW, Shaw DE: Activation mechanism of the beta2-adrenergic receptor. Proc Natl Acad Sci USA 2011, 108: 18684–18689. 10.1073/pnas.1110499108PubMed CentralView ArticlePubMedGoogle Scholar
- Nygaard R, Zou Y, Dror RO, Mildorf TJ, Arlow DH, Manglik A, Pan AC, Liu CW, Fung JJ, Bokoch MP, Sun TT, Shaw DE, Mueller L, Prosser RS, Kobilka BK: The dynamic process of β 2 -adrenergic receptor activation. Cell 2013, 152: 532–542. 10.1016/j.cell.2013.01.008PubMed CentralView ArticlePubMedGoogle Scholar
- Rosenbaum DM, Cherezov V, Hanson MA, Rasmussen SG, Thian FS, Kobilka TS, Choi HJ, Yao XJ, Weis WI, Stevens RC, Kobilka BK: GPCR engineering yields high-resolution structural insights into beta2-adrenergic receptor function. Science 2007, 318: 1266–1273. 10.1126/science.1150609View ArticlePubMedGoogle Scholar
- O’Dowd BF, Hnatowich M, Regan JW, Leader WM, Caron MG, Lefkowitz RJ: Site-directed mutagenesis of the cytoplasmic domains of the human beta 2-adrenergic receptor. Localization of regions involved in G protein-receptor coupling. J Biol Chem 1988, 263: 15985–15992.PubMedGoogle Scholar
- Liggett SB, Caron MG, Lefkowitz RJ, Hnatowich M: Coupling of a mutated form of the human beta 2-adrenergic receptor to Gi and Gs. Requirement for multiple cytoplasmic domains in the coupling process. J Biol Chem 1991, 266: 4816–4821.PubMedGoogle Scholar
- Isin B, Estiu G, Wiest O, Oltvai ZN: Identifying ligand binding conformations of the β2-adrenergic receptor by using its agonists as computational probes. PLoS One 2012, 7(12):e50186. 10.1371/journal.pone.0050186PubMed CentralView ArticlePubMedGoogle Scholar
- Katritch V, Reynolds KA, Cherezov V, Hanson MA, Roth CB, Yeager M, Abagyan R: Analysis of full and partial agonists binding to beta(2)-adrenergic receptor suggests a role of transmembrane helix V in agonist-specific conformational changes. J Mol Recognit 2009, 22: 307–318. 10.1002/jmr.949PubMed CentralView ArticlePubMedGoogle Scholar
- Simpson LM, Wall ID, Blaney FE, Reynolds CA: Modeling GPCR active state conformations: the beta(2)-adrenergic receptor. Proteins-Structure Function and Bioinformatics 2011, 79: 1441–1457. 10.1002/prot.22974View ArticleGoogle Scholar
- Amadei A, Linssen ABM, Berendsen HJC: Essential dynamics of proteins. Proteins-Structure Function and Genetics 1993, 17: 412–425. 10.1002/prot.340170408View ArticleGoogle Scholar
- Kurkcuoglu Z, Bakan A, Kocaman D, Bahar I, Doruker P: Coupling between catalytic loop motions and enzyme global dynamics. Plos Comput Biol 2012, 8(9):e1002705. 10.1371/journal.pcbi.1002705PubMed CentralView ArticlePubMedGoogle Scholar
- MATLAB 220.127.116.119(R2010a). Natick, Massachusetts: The MathWorks Inc; 2010.
- Laskowski RA: PDBsum new things. Nucleic Acids Res 2009, 37: D355-D359. 10.1093/nar/gkn860PubMed CentralView ArticlePubMedGoogle Scholar
- Swaminath G, Lee TW, Kobilka B: Identification of an allosteric binding site for ZN(2+) on the beta(2) adrenergic receptor. J Biol Chem 2003, 278: 352–356.View ArticlePubMedGoogle Scholar
- MOE 2011.10. 1010 Sherbooke St. West, Suite #910, Montreal, QC, Canada, H3A2R7: Chemical Computing Group Inc; 2011.
- The PyMOL Molecular Graphics System P. 0.99. Schrödinger: LLC;
- Schlessinger A, Punta M, Yachdav G, Kajan L, Rost B: Improved disorder prediction by combination of orthogonal approaches. PLoS ONE 2009, 4(2):e4433. 10.1371/journal.pone.0004433PubMed CentralView ArticlePubMedGoogle Scholar
- Linse SS: The nobel prize in chemistry 2012 - advanced information. mm 2013. http://www.nobelprize.org/nobel_prizes/chemistry/laureates/2012/advanced-chemistryprize2012.pdf Google Scholar
- Monod J, Wyman J, Changeux JP: On the nature of allosteric transitions: a plausible model. J Mol Biol 1965, 12: 88–118. 10.1016/S0022-2836(65)80285-6View ArticlePubMedGoogle Scholar
- Weber G: Ligand binding and internal equilibriums in protein. Biochemistry 1972, 11: 864–878. 10.1021/bi00755a028View ArticlePubMedGoogle Scholar
- Ma B, Kumar S, Nussinov R: Folding and binding cascades: shifts in energy landscapes. Proc Natl Acad Sci U S A 1999, 96(18):9970–9972. 10.1073/pnas.96.18.9970PubMed CentralView ArticlePubMedGoogle Scholar
- Narayanan E, John B, Mirkovic N, Fiser A, Ilyin VA, Pieper U, Stuart AC, Marti-Renom MA, Madhusudhan MS, Yerkovich B, Sali A: Tools for comparative protein structure modeling and analysis. Nucleic Acids Res 2003, 31: 3375–3380. 10.1093/nar/gkg543View ArticleGoogle Scholar
- Humphrey W, Dalke A, Schulten K: VMD: visual molecular dynamics. J Mol Graph 1996, 14: 33–38. 10.1016/0263-7855(96)00018-5View ArticlePubMedGoogle Scholar
- Lomize MA, Lomize AL, Pogozheva ID, Mosberg HI: OPM: orientations of proteins in membranes database. Bioinformatics 2006, 22: 623–625. 10.1093/bioinformatics/btk023View ArticlePubMedGoogle Scholar
- Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kalé L, Schulten K: Scalable molecular dynamics with NAMD. J Comput Chem 2005, 26: 1781–1802. 10.1002/jcc.20289PubMed CentralView ArticlePubMedGoogle Scholar
- Petrache HI, Dodd SW, Brown MF: Area per lipid and acyl length distributions in fluid phosphatidylcholines determined by (2)H NMR spectroscopy. Biophys J 2000, 79: 3172–3192. 10.1016/S0006-3495(00)76551-9PubMed CentralView ArticlePubMedGoogle Scholar
- Mackerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Roux B, Schlenkrich M, Smith J, Stote R, Straub J, Watanabe M, Wiorkiewicz-Kuczera J, Yin D, Karplus M: Self-consistent parameterization of biomolecules for molecular modeling and condensed phase simulations. Faseb Journal 1992, 6: A143-A143.Google Scholar
- Mackerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Roux B, Schlenkrich M, Smith J, Stote R, Straub J, Watanabe M, Wiorkiewicz-Kuczera J, Yin D, Karplus M: All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B 1998, 102: 3586–3616. 10.1021/jp973084fView ArticlePubMedGoogle Scholar
- Schlenkrich M, Brickmann J, MacKerell AD Jr, Karplus M: An empirical potential energy function for phospholipids: criteria for parameter optimization and applications. In Biological Membranes: A Molecular Perspective from Computation and Experiment. 1st edition. Edited by: Merz KMJr, Roux B. Birkhauser, Boston; 1996.Google Scholar
- Feller SE, Yin D, Pastor RW, MacKerell AD Jr: Molecular dynamics simulation of unsaturated lipids at Low hydration: parametrization and comparison with diffraction studies. Biophys J 1997, 73: 2269–2279. 10.1016/S0006-3495(97)78259-6PubMed CentralView ArticlePubMedGoogle Scholar
- Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML: Comparison of simple potential functions for simulating liquid water. J Chem Phys 1983, 79: 926–935. 10.1063/1.445869View ArticleGoogle Scholar
- Feller SE, Zhang YH, Pastor RW: Computer-simulation of liquid/liquid interfaces 2. Surface-tension area dependence of a bilayer and monolayer. J Chem Phys 1995, 103: 10267–10276. 10.1063/1.469928View ArticleGoogle Scholar
- Ryckaert JP, Ciccotti G, Berendsen HJC: Numerical-integration of cartesian equations of motion of a system with constraints - molecular-dynamics of N-alkanes. J Comput Phys 1977, 23: 327–341. 10.1016/0021-9991(77)90098-5View ArticleGoogle Scholar
- Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, Olson AJ: AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility. J Comput Chem 2009, 30: 2785–2791. 10.1002/jcc.21256PubMed CentralView ArticlePubMedGoogle Scholar
- Doruker P, Atilgan AR, Bahar I: Dynamics of proteins predicted by molecular dynamics simulations and analytical approaches: application to alpha-amylase inhibitor. Proteins 2000, 40: 512–524. 10.1002/1097-0134(20000815)40:3<512::AID-PROT180>3.0.CO;2-MView ArticlePubMedGoogle Scholar
- Atilgan AR, Durell SR, Jernigan RL, Demirel MC, Keskin O, Bahar I: Anisotropy of fluctuation dynamics of proteins with an elastic network model. Biophys J 2001, 80: 505–515. 10.1016/S0006-3495(01)76033-XPubMed CentralView ArticlePubMedGoogle Scholar
- Lezon TR, Bahar I: Constraints imposed by the membrane selectively guide the alternating access dynamics of the glutamate transporter GltPh. Biophys J 2012, 102: 1331–1340. 10.1016/j.bpj.2012.02.028PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.