Contributed by Arieh Warshel, September 3, 2020 (sent for review June 29, 2020; reviewed by Timothy Clark and D. Thirumalai)
A detailed molecular-level understanding of the mechanism of activation and coupling selectivity of the Î¼-opioid receptor (Î¼-OR) is largely missing in the literature, despite remarkable advances in structural and biochemical studies. Here we have performed extensive free energy calculations using our electrostatic-based coarse-grained model and obtained a complete thermodynamic cycle and two-dimensional free energy landscapes for the activation process. This allowed us to propose a plausible sequence of events that drives the activation process. Additionally, we have been able to explain the selectivity of the coupling between Î¼-OR and G proteins, by using free energy calculations. This advance can serve as a basis for further progress in understanding the biased signaling of Î¼-OR.
Understanding the activation mechanism of the Î¼-opioid receptor (Î¼-OR) and its selective coupling to the inhibitory G protein (Gi) is vital for pharmaceutical research aimed at finding treatments for the opioid overdose crisis. Many attempts have been made to understand the mechanism of the Î¼-OR activation, following the elucidation of new crystal structures such as the antagonist- and agonist-bound Î¼-OR. However, the focus has not been placed on the underlying energetics and specificity of the activation process. An energy-based picture would not only help to explain this coupling but also help to explore why other possible options are not common. For example, one would like to understand why Î¼-OR is more selective to Gi than a stimulatory G protein (Gs). Our study used homology modeling and a coarse-grained model to generate all of the possible âend statesâ of the thermodynamic cycle of the activation of Î¼-OR. The end points were further used to generate reasonable intermediate structures of the receptor and the Gi to calculate two-dimensional free energy landscapes. The results of the landscape calculations helped to propose a plausible sequence of conformational changes in the Î¼-OR and Gi system and for exploring the path that leads to its activation. Furthermore, in silico alanine scanning calculations of the last 21 residues of the C terminals of Gi and Gs were performed to shed light on the selective binding of Gi to Î¼-OR. Overall, the present work appears to demonstrate the potential of multiscale modeling in exploring the action of G protein-coupled receptors.
The opioid overdose crisis is one of the leading causes of accidental death in the United States. Opioids, which are clinically used as analgesic drugs, bind to various opioid receptors to activate the inhibitory G protein (Gi) signaling pathway and thus initiate its analgesic functions. The Î¼-opioid receptor (Î¼-OR) is the main opioid target for pain management (1). It can also interact with Î²-arrestin through a different activation pathway and may induce adverse effects like respiratory depression (2). It is found that some opioid ligands can preferentially favor Gi over Î²-arrestin signaling and thus reduces the adverse effect (3). Therefore, a detailed understanding of the mechanism of activation and identification of the sequential conformational changes, that lead to the activation, can have immense pharmaceutical implications.
Î¼-OR is one of only a few G protein-coupled receptors (GPCRs) for which high-resolution structures are available, in both the active (4, 5) and inactive forms (6). These structural discoveries have laid the foundations for numerous studies that strived to explain the thermodynamic and kinetic aspects of the activation process (7ââââ11) as well as to develop new drugs (3, 12, 13). Recent spectroscopic studies demonstrated that Î¼-OR exists in a conformational equilibrium, where different external stimuli (agonist/antagonist) stabilize different intermediate states to induce the downstream signaling (14, 15). In general, the complete free energy landscape of the receptor activation is essential for connecting the structural plasticity of GPCR to its function (16). The biased transduction of the signaling process by GPCRs in the presence of different external stimuli can be connected to the changes in the free energy landscape of activation in the given physiological conditions. Thus, by free energy landscape calculations we can potentially explain the bias signaling of Î¼-OR.
The sequence of motions of the domains in the GPCRâG protein complex, that propagate the signal from the extracellular side of the receptor to the intracellular side, can be probed by experimental studies (17). However, a more detailed energy-based molecular picture requires extensive analysis of the free energy surface, using computer simulations. The identification of the relevant intermediate structures that connect the energy minima (for the active and inactive states) on the energy landscape should be a subject of major interest. Such an identification is important for understanding the biased functionality of GPCRs (8) and for future drug development.
Apparently, proper understanding of the activation process of the receptor is incomplete without considering the whole ternary complex of the ligandâreceptorâG proteins. Although, there are some excellent recent computational studies where the ternary complexes were considered (18), the major challenge for incorporating the whole ternary complex in simulations is in obtaining adequate sampling of such an exceptionally large system (the ternary complex) (19). A multimicrosecond all-atom molecular dynamics (MD) simulation has suggested a mechanism of nucleotide release from G proteins but did not determine how the activation of the receptor affects the sequence of events for the nucleotide release. The abovementioned process can be simulated consistently only when the effect of the complete receptorâG protein association is considered. For example, in another study, all three components of the ternary complexes were considered in simulating the biased coupling of the Î²2-AR receptor to different intracellular binding partners (IBPs) (in the presence of different ligands), using enhanced-sampling simulations (18). This computational work has paved an avenue for understanding the effect of cooperativity of different binding components in the formation of the complexes, but it is not clear whether the opening (an integral part of G protein activation) of the GsÎ± protein has been simulated consistently (18). Some computational (20) and experimental (21, 22) works support the idea that the G protein opening is a necessary but not sufficient condition for the overall activation process. Therefore, a consistent calculation of the energetics of the activation process not only needs all of the components of the trimeric system (the receptor, the IBPs, and the G proteins) but also requires consideration of all of the relevant end-point states and the transition among those states during the process of activation (23). Additionally, identifying the sequence of events that leads to the activation of GPCRâG protein (for example, Î¼-ORâGi) signaling pathway should enhance our ability to identify relevant intermediates along the activation process and to understand the important interactions between the extracellular ligands and Î¼-OR that stabilizes such intermediates. This should help in designing new drugs that preferentially activate G protein signaling pathways or reduce adverse effects of opioids. Of course, such advances require a tight combination with the corresponding experimental information.
Another important aspect that enhances our understanding of Î¼-OR transduced signaling pathways is the selectivity of Î¼-OR to a specific G protein. The downstream beneficiary signal transduction occurs mainly via the Î¼-ORâGi. Different GPCRs can selectively interact with one or more subtypes of the Î±-subunit of the G protein families and produce distinct cellular responses. The mechanisms of the selective GPCRâG protein coupling are still not clearly understood. The analysis of the primary sequences and structures of different GPCRâG protein complexes reveals the presence of a distinct selectivity barcode residues for different GÎ± and without a simple conserved motif of GPCR for the selective coupling (24). Interestingly, the promiscuity of many GPCRs in the coupling to their G proteins varies not only between different GPCR families but also within members of the same families (25). Thus, receptor-specific analysis of this selective coupling is very important. The analysis of one-dimensional sequences of the complexes or static three-dimensional structures are very informative for the coupling selectivity, but experiments like alanine scanning (26) or computer simulations (27) should be very important in providing microscopic understanding of this phenomena.
This study took Î¼-OR as a model system, due to its importance in the context of fighting the ongoing opioid overdose crisis. A complete thermodynamic cycle of the Î¼-ORâGi activation was generated, by evaluating the free energy change of each step using our electrostatic-based coarse-grained (CG) model (28). Next, the end points of the thermodynamic cycle were used to generate intermediate states along the activation pathways, by using âtargeted molecular dynamicsâ (TMD) or, alternatively, by our newly developed protocol for exploring the path between the endpoint conformations. The intermediates were used to simulate detailed two-dimensional (2D) landscapes of the activation process to provide a plausible pathway of the activation. Finally, the comparison of the free energy of stabilization of different Î¼-OR-G protein complexes and computational alanine scanning of the last 21 C-terminal residues of the G proteins (Gi and Gs) helped us to shed light on the selective binding of Gi to Î¼-OR.
Results and Discussion
The Thermodynamic Cycle for the Activation of the Î¼-OR Complex System.
In order to understand the activation mechanism of a process, it is very useful to start by computing the energies of all the âend statesâ of the relevant thermodynamic cycle. This asymptotic treatment can guide the detailed exploration of possible mechanistic options. Additionally, the relative stability of the âend statesâ can be employed in estimating the accuracy of the modeling protocol.
This work considered two âend statesâ of Î¼-OR (active and inactive) and four states of GiÎ± (shown in Fig. 1) in defining the relevant thermodynamic cycle. The Î±-subunits of the G proteins consists of two domains, namely âRas-like domainâ (RD) and âalpha-helical domainâ (AHD). These domains become separated when the G protein is bound to a GPCR, whereas in the absence of the receptor there is generally no separation (22). The domain-separated conformation of the G protein is termed here âopen,â while the other is named âclosed.â The C-terminal (Î±5-helix) of the Î±-subunit of the G protein (circled in Fig. 1) also changes the structure from disordered (dynamic) to ordered (rigid) upon receptor binding (29). In comparison to the receptor-free state of the G protein, its Î±5-helix in the receptor-bound state moves upward (with respect to the nucleotide binding site) toward the receptor (Fig. 1). Such upward and ordered conformation of the G protein is termed here as âinâ and the disordered conformation as âout.â Thus, in the entire activation process of any G protein, two drastically different conformational changes are observed: 1) the separation of the RD and AHD (âclosedâ to âopenâ) and 2) the upward movement and disordered-to-ordered structural change of the Î±5-helix (âoutâ to âinâ). Subsequently, the four âend statesâ of G protein are named 1) âclosed-out,â 2) âopen-out,â 3) âclosed-in,â and 4) âopen-in.â Ideally, for each âend stateâ of Î¼-OR there are four Gi protein âend states,â and thus we should consider a thermodynamic cycle with a total of 16 âend statesâ of the complex. Four out of those 16 states, for example the bound and unbound complexes of the inactive Î¼-OR with the âclosed-inâ and âopen-inâ conformation of GiÎ±, are not included in the thermodynamic cycle. The âinâ state of GiÎ± forms a very unstable complex, with the âinactiveâ form of Î¼-OR, and thus those four states should not be parts of any feasible activation pathway. For each âend stateâ in the thermodynamic cycle shown in Fig. 2, the free energy of a state was estimated by calculating the folding free energy of the state by means of our CG model. The corresponding procedure for generating all âend stateâ complexes is provided in SI Appendix, section S2.1.1. The thermodynamic cycle calculations were performed without considering the agonist as part of the complex. However, the contribution of agonists to the thermodynamic cycle of the activation is discussed in SI Appendix, section S3.
The most interesting findings from the calculated thermodynamic cycle are summarized below.
The binding of GiÎ± to Î¼-OR is always favorable (AâG, BâH, CâI, DâJ, EâK, and FâL), although the probability is higher when Î¼-OR is in an active-like conformation (BâH, CâI, DâJ, and EâK). This suggests that Î¼-OR can form a complex with GiÎ± in the absence of an agonist. In fact, a bioluminescence resonance energy transfer study supports this, by reporting that Î´-OR forms a constitutive complex with a heterotrimeric Gi in the absence of agonist (30). Fig. 2 also shows that the binding affinity of GiÎ± to Î¼-OR increases when the GiÎ± is in the âopen-inâ conformation (DâJ), in accordance with the observation that in the crystal structures of the active complexes the G protein is always in the âopen-inâ conformation. The calculated binding free energy for (DâJ) seems to be overestimated with respect to the experimental binding affinity of the GÎ±i1-hÎ´-OR (unliganded), calculated using plasmon-waveguide resonance spectroscopy (31). However, the entropy loss accompanying the formation of a complex of two large proteins could be quite significant. We note in this respect that this phenomenon is not properly addressed by our CG-based free energy calculations and thus the loss of entropy during the binding process is expected to be underestimated.
The calculated thermodynamic cycle also implies that the opening of GiÎ± (outward movement of the AHD) is not favorable without the âinâ conformation of Gi (IâJ vs. HâK) and that the Gi should be bound to Î¼-OR to favor the opening (CâD vs. IâJ). Contrary to these results, the authors of ref. 20 suggested that the separation of RD and AHD should be spontaneous in the absence of the nucleotide, irrespective of the presence or absence of the receptor. In this case it is important to understand that the AHD part of the âopenâ conformation of the GiÎ± was modeled with respect to the AHD of GsÎ± and consequently the extent of opening is around 60 Ã in the modeled structures (like in the âopenâ conformation of GsÎ±). Unfortunately, no activated G protein structure, other than Protein Data Bank (PDB) ID code 3SN6, is available as a reference for generating GiÎ± with different extent of opening. It is possible that such a large opening of the Gi is unreasonable, and thus possibly we obtained an overestimation in the calculated energy of the âend statesâ with the âopenâ conformation of Gi. It is also possible that the extent of separation in Gi is different in the receptor-bound and -unbound states of the Gi, which would be hard to guess without more structural information.
Fig. 2 also shows that the stability of the active form of Î¼-OR (without the agonist) can be increased upon binding of GiÎ±, which subsequently allows the opening of the G protein (AâGâHâIâJ). This phenomenon is consistent with the experimentally found basal activity of Î¼-OR (active in the absence of an agonist) (32). It is worth mentioning that a transition of an inactive form of receptor to an active form can be simulated by long MD simulations and the activation state of a given receptor conformation can be detected from those simulations by a distance-based activation index model (33), but free energy calculations and consideration of a complete thermodynamic cycle of activation provides clearer energy-based explanations of the processes, like rationalization of the basal activity.
Overall, our thermodynamic cycle is very informative and quite successful in explaining several important experimental observations, especially the basal activity of Î¼-OR. The cycle can be helpful for comparing the relative stability of different states, and this information can be used in generating a qualitative energy landscape, as was done in ref. 23. However, a more comprehensive analysis of the sequence of conformational changes in the Î¼-ORâGi complex during the activation, can only be obtained by performing extensive 2D free energy landscape calculations.
Two-Dimensional Plots for Exploring Plausible Activation Pathways.
Our calculated 2D free energy landscapes were generated in order to describe the probable sequence of conformational changes leading to the activation of the Î¼-ORâGi system. The questions that we would like to answer by this analysis are as follows: 1) What is the most important conformational change in Gi (âinâ-âoutâ or âopenâ-âclosedâ) for the activation process; 2) how can one conformational change in Gi facilitate the other; and 3) what is the sequence of conformational changes that leads to the activation of Gi in the presence of Î¼-OR? To address these questions we performed four different free energy landscape calculations using our CG model. In the 2D free energy landscapes, one reaction coordinate corresponds to the change from the âactiveâ to âinactiveâ state of Î¼-OR and the other one indicates the changes in the Gi (for example, the change from the âclosed-outâ to âopen-inâ conformation and its variations). A detailed description of the methods used to generate the 2D landscape is provided in SI Appendix, section S2.1.2.
The 2D free energy landscape in Fig. 3 indicates that the transition from the âoutâ to âinâ conformation of the G protein is facilitated when the Î¼-OR is in the âactiveâ like conformation. In the most inactive form of Î¼-OR, the folding free energy of the complex is rather high, regardless the conformation of the G protein. Our calculations showed that the stability of the complex increases only when the TM6 of Î¼-OR is sufficiently open (one of the indicators of the âactiveâ conformation of the receptor). Note that the AHD of the GiÎ± was not included in the calculation for generating the plot in Fig. 3. Thus, Fig. 3 can be used to obtain a picture of the activation pathway where the effect of AHD on the activation process is missing. On the other hand, Fig. 3 can explain the formation of a stable complex of an active Î¼-OR with mini-Gi (GiÎ± with no AHD) and proves that the stability of the active conformation is not necessarily dependent on the opening of AHD (34ââ36). A recent study also shows that an opioid receptor can bind mGsi (a mini-Gs in which the last nine residues of the C terminal are of Gi) (37) and form a stable activated receptorâG protein complex. The analysis of the plot also hints that in general our 2D free energy landscape calculations could be used successfully (by also considering the motion of the AHD part) for explaining the sequence of changes that leads to the activation of the Î¼-OR-(full) Gi system.
The addition of the AHD to Fig. 3 would make the transformation coordinate along the x axis of the landscape complex, since two distinct conformational changes of the G protein have to be simulated. Thus, to understand the effect of the AHD part on the activation process, these changes were simulated separately and generated three different 2D surface by the following approach: 1) The AHD part of the G protein was included in the calculation, but it was assumed that no domain separation occurred (no opening of G protein); 2) the AHD domain separation was simulated in the presence of the âoutâ conformation of the Î±5-helix of Gi; and 3) the AHD domain separation was simulated in the presence of the âinâ conformation of the Î±5-helix of Gi.
Fig. 4 describes the effect of the âclosedâ form of the AHD on Fig. 3. It can be seen from Fig. 4 that the parts of the complexes, where the receptor is in the âinactiveâ state, became marginally less stable after adding the AHD part to the RD part of the GiÎ±, relative to similar complexes in Fig. 3 (compare the top right parts of Figs. 3 and 4). This result is consistent with the experimental conclusion that the close proximity of the two domains of GiÎ± is less favorable when the Î±5-helix of GiÎ± moves toward the activated receptor (20, 22). On the other hand, the overall relative stabilities of the protein complexes are not changed significantly by adding the AHD in the âclosedâ form with respect to similar complexes without the AHD. For example, the relative difference in folding free energy between the most stable section of the landscape and the major of part of Fig. 3 (orange coloration in the plot) are â20 kcal/mol. The same is true for Fig. 4 (light orange coloration). Interestingly, it can also be concluded from Fig. 4 that the upward movement of the Î±5-helix precedes the domain separation and stabilizes the active conformation of the receptor.
Fig. 5 may lead to the conclusion that the opening of the G protein is relatively difficult (due to the high free energy barrier) when the Î±5-helix is in the âoutâ conformation irrespective of the conformation of the receptor. Another conclusion that can be reached based on the plot is that even when the Î¼-OR is in the âactiveâ conformation there is a population of the very unstable complexes. On the other hand, the bottom right corner of Fig. 5 shows that there is a region on the energy surface corresponding to a stable receptorâG protein complex formation. Therefore, the results in Fig. 5 indicate that the domain separation (with Î±5-helix in the âoutâ conformation) requires crossing a barrier. Note that to simulate the AHD and RD domain separation we had to use TMD. It is possible that the effect of the bias (0.5 kcal/mol) used to generate the intermediates between âclosed-outâ and âopen-outâ had not been removed completely by the relaxation (0.5 ns) performed after the execution of the TMD for each intermediate structure. As a result, we might have obtained high-energy intermediate conformations between the âclosed-outâ and the âopen-outâ conformation. Finally, our results in Fig. 6 are useful in suggesting how important the TMD simulations are in providing useful qualitative results about the pathway of the activation of the Î¼-ORâGiÎ± complex. This is true despite the possibility that the TMD might not always generate low-energy intermediates. This problem could have been resolved by our renormalization approach (38).
Interestingly, by comparing Figs. 5 and 6 it can be concluded that the energy barrier for the opening of the G protein is reduced when the G protein is in its âinâ conformation, and the G protein formed a more stable complex with the active-like receptor conformation. Thus, the âoutâ-to-âinâ conformation change should start before the domain separation (see Fig. 4), or in other words the upward movement of the Î±5-helix reduces the contacts between the Ras domain and the AHD domain and facilitates the domain separation. This is also in accordance with the conclusion of ref. 20.
Overall, our thermodynamic cycle and the semiquantitative 2D free energy landscapes can provide a tentative sequence of conformational changes that rationalize the pathway for the basal activity of the Âµ-ORâGi complex. From the thermodynamic cycle we know that the binding of the Gi protein to Âµ-OR is favorable for all six end-point states (G, H, I, J, K, and L states in Fig. 2). It is also known from experiment that the G proteins are most stable in the âclosed-outâ conformation, when they are not bound to the GPCRs. Consequently, the plausible sequence of changes during the activation is the following. The process starts with an initial weakly interacting complex formation between the âclosed-outâ conformation of Gi and inactive form of Âµ-OR. Next, the Gi tends to extend the Î±5-helix toward the intracellular cavity of the receptor to form a âclosed-inâ conformation. Here the âclosed-outââtoââclosed-inâ transition also favors the shifting of the equilibrium from the âinactiveâ to âactiveâ conformation of Âµ-OR. From Figs. 3 and 4 we can conclude that the transition between the âinactiveâ and âactiveâ conformations of Âµ-OR slightly precedes the âoutâ-to-âinâ conversion of Gi. Alternatively, we may say that the cooperative effect of both Âµ-OR and Gi stabilizes the âactiveâ Âµ-OR âclosed-inââGi complex. The Âµ-OR âclosed-inââGi complex formation then favors the domain separation of the Gi and thus leads to the formation of the active Âµ-OR âopen-inââGi complex. Thus, our extensive 2D free energy landscape free energy calculation can provide a qualitative understanding of the basal activation pathway and supports many of the experimental findings. However, accurate activation barrier calculations for the transitions should be simulated to infer the kinetics of the transition processes and to obtain a more concrete understanding of the sequence of events in the activation process. At this point, we would like to mention that the kinetics of the transitions are very important for understanding the fate of a transformation process, and in particular the transition between parallel pathways. However, a consistent calculation of the kinetics of transition requires the finding of distinguishable pathways. Unfortunately, while our âshortest path findingâ method produces sufficiently smooth transitions it does not guarantee that it follows the best least-energy paths. Thus, we cannot at this stage infer the kinetics of the transition process with high confidence, although the 2D landscape by itself can be used to get a qualitative idea of the kinetics of the processes.
Selectivity of the Coupling in the Î¼-ORâGi Complex.
Î¼-OR selectively forms complexes with Gis to transduce the downstream signaling, but the molecular mechanism of this specificity of Î¼-OR toward Gi is elusive. As pointed out in the Introduction, a receptor-specific analysis is necessary to discover the selectivity determinant for Î¼-OR but is largely missing in the literature. On the other hand, consistent free energy calculations that compare hypothetical complexes with canonical conformation of Î¼-OR can be used to explore the influence of the canonical conformation of Î¼-OR on its selectivity. It is mentioned in refs. 5, 39, and 40 that the probable reason for the observed selectivity is the relatively small opening of the transmembrane helix 6 (TM6) in Î¼-OR. According to this idea, the smaller opening reduces the sterically accessible space for the accommodation of the bulky residues at the C-terminal of the GÎ± that are more common in other G proteins (like Gs), compared to those in Gi. An example that supports the above-mentioned deduction is the finding that the receptor Î²2-AR couples preferentially to Gs and the active form of the receptor has much wider opening at the cytoplasmic side of the receptor compared to that of active Î¼-OR (SI Appendix, Fig. S3). Thus, we may wonder whether an active Î¼-OR can take a more extended conformation than the one already found in the crystal structure and can accommodate Gs. It is known that the structural plasticity of the GPCRsâ cytoplasmic pocket is high in the presence of an agonist (41). Therefore, in an ideal scenario Î¼-OR in its agonist-bound form may assume a conformation with the cytosolic pocket as wide as in receptor Î²2-AR to accommodate Gs. However, it is already known that Î¼-OR can be very selective to Gi. In this case the probable reasons for the observed Gi selectivity are 1) Î¼-OR cannot adopt any extended conformation to accommodate the bulky C-terminal of G proteins (like Gs) and 2) even if an extended conformation is energetically accessible, the complex formation with Gs is not stable. We were interested in checking what would be the most probable reason for the G protein selectivity shown by Î¼-OR.
By assuming that an extended conformation of the active form of Î¼-OR (ext. Î¼-OR) is possible, we may try to see if that conformation can stabilize the Î¼-ORâGs complex more than that with Gi. A comparison of the stability of different complexes where G proteins are bound to a crystal structure like conformation of Î¼-OR or ext. Î¼-OR can be used to evaluate the above-mentioned hypothesis.
In order to generate a model for the ext. Î¼-OR, we have considered the structure of the active form of Î²2-AR (a Gs binding GPCR) as a template. Since the sequence identity between Î²2-AR and Î¼-OR is small, the directly modeled structure could not be used for reliable comparison. Hence, we propagated a relaxation run and processed the last few frames of the 80-ns-long production trajectory. An inward movement of the TM6 was observed after the MD simulation, while it remains sufficiently open to be able to accommodate bulkier Gs. The relaxed structures of the active Î¼-OR and the ext. Î¼-OR were combined with the mini-Gi and Gs proteins to produce four types of the structures of the complex (SI Appendix). The structures generated were used to support the calculation of the folding free energy (related to the stability of the complexes). The reason for involving mini-G proteins in these calculations is discussed in SI Appendix, section S2.1.3. The results of the free energy calculations (Table 1) show that the stability of ext. Î¼-OR-mini-Gi does not change considerably with respect to the Î¼-ORâmini-Gi complex. The decrease in the stability in ext. Î¼-ORâmini-Gi is supported by the fact the ext. active Î¼-OR conformation is less feasible for Î¼-OR in the presence of Gi. The substantial destabilization for the Î¼-ORâmini-Gs complex, compared to Î¼-ORâmini-Gi, is also expected, since Î¼-OR in the active like conformation has small cytosolic opening to accommodate the bulky C-terminal of Gs (35). Surprisingly, Table 1 reveals that the complex of the ext. Î¼-OR with mini-Gs is even less stable compared to the Î¼-ORâmini-Gs complex. This provides a proof that the highly stabilizing intermolecular interaction between mini-Gi and Î¼-OR is responsible for the selectivity. In fact, even if the intracellular side of the Î¼-OR is open enough to accommodate Gs, the resulting complex would not be stable. It is worth mentioning that we tried to generate two different Î¼-ORâmini-Gs complex structures, since the position of Î±5-helix of Gi and Gs in their respected GPCR complexes (PDB ID codes 6DDF and 3SN6) are aligned differently (5). The calculated stabilities of the complexes formed by those two different reference structures are similar (Table 1). That result suggests that the conclusion about the stability of the Î¼-ORâmini-Gs complex remains the same, regardless of the references used to model the structures of Î¼-ORâmini-Gs. Hence, our energy-based calculations reveal that even if the bulky residues can fit inside some other active like conformations of the Î¼-OR, the specificity is related to more residue-specific interactions occurring between the Î¼-OR and the Gi.
To check in what way the residue-specific interactions between the Î¼-OR and G protein act as a selectivity determinant, we estimated the residue-based interaction energy between the G proteins and the Î¼-OR. The corresponding computational alanine scan on the last 21 residues of the C-terminal of the G proteins can help us find out how the C-terminal residues act as a selectivity determinant for the coupling specificity. The scan on C-terminal residues is promising since the C terminal of a G protein is widely believed to be a primary determinant for the coupling selectivity. Moreover, it is also found in the literature that the last five residues at the C terminal of the G proteins are most important for selectivity (42), and even a single mutation at those positions can alter the selective binding of the G protein to wild-type GPCRs (43). Additionally, a recent study reports that the active conformation of the Î¼-OR is stabilized by the last 12 residues of GiÎ± (44).
Fig. 7 depicts the relative difference in solvation free energy of changing the residues at each position of the last 21 residues of the GsÎ± and GiÎ± to alanine at that position. These free energies are calculated based on the semimicroscopic protein dipole Langevin dipole (PDLD/S-LRA) method (45), provided by the MOLRIS-XG software (46). The difference in generalized solvation free energy shows how stable a given amino acid is at its own position. A positive difference in the solvation free energy is an indicator that the residue (X) is having stabilizing interactions with the rest of the complex, and an alanine mutation at that position would contribute to the destabilization of the overall complex. As seen in Fig. 7, for most of the positions of Gi the stability of the complex decreases by alanine mutations [positive ]. Only a few residues are showing negative values of , but the extent of stabilization (by alanine mutation) at those sequence positions is also very small. Contrary to Gi, only for a few residues in Gs is the value positive. Thus, by comparing the last 21 residues of GsÎ± and GiÎ± it can be concluded that these residues are integral parts of the selectivity determinant region. A selectivity analysis of this kind can be also used to explain the biased signaling mechanism. The difference in molecular interactions at the interface region of GPCR and its intracellular partners (Î²-arrestin and G protein) can be recognized by mutating different residues on Î¼-OR at the interface region. However, this latter analysis would require structural information about the Î¼-ORâÎ²-arrestin complex which is currently not available.
Here we have explored two aspects of the Âµ-ORâGi protein system. First, we tried to analyze the activation mechanism of the receptorâG protein complex that helps in transducing signaling from the extracellular part of the cell to its intracellular part. Second, we suggested a molecular basis for a specific coupling of the active Âµ-OR and Gi. The free energy landscape calculations showed that the upward movement and the disordered-to-ordered transition of the Î±5-helix of Gi is very important for the shift of equilibrium between the âinactiveâ and âactiveâ forms of the receptor. The change in the Î±5-helix of Gi also favors the domain separation of GiÎ±. Our calculated thermodynamic cycle of the activation process also explains the experimentally observed basal activity of Âµ-OR. The computational alanine scanning and folding free energy calculations of the modeled GPCRâG protein complexes provide comprehensive free energy-based understanding of the molecular basis of specific coupling between the active Âµ-OR and Gi. The folding free energy calculations revealed that any large extended Âµ-OR conformation can accommodate the C terminal of GsÎ± but cannot form more stable complex than Âµ-ORâGiÎ±. Furthermore, the all-atomistic relative solvation free energy calculations have shown that the alanine mutations at different positions of the last 21 residues largely destabilize the Âµ-ORâGi complex compared to that of Âµ-ORâGs. This last finding allowed us to conclude that the last 21 residues of G proteins can act as a selectivity barcode that Âµ-OR searches to find the correct binding partner to transduce downstream signaling. These semiquantitative understandings of the activation mechanism and selectivity to Gi can be crucial in rational design of drugs to reduce the opioid overdose crisis.
In the future, a consistent implementation of our approaches to other GPCRâG protein complexes can provide an energy-based explanation of any universal mechanism of activation in a class of GPCRs, as it has already been supported by other computational-, structural-, or evolutionary-based comparisons (47ââââ51). Additionally, our adopted approach can be extended for revealing the mechanism of the biased signaling by studying the influence of different extracellular binding partners on the free energy landscape. While the effect of an agonist, [D-Ala2, N-MePhe4, Gly-ol]-enkephalin (DAMGO), on the thermodynamic cycle of activation of Âµ-OR is discussed in SI Appendix, Fig. S4, it would be interesting to investigate the effect of agonist binding on 2D free energy landscapes (like Figs. 3â6). That study requires huge number of calculations. Those calculations include 1) docking studies to find proper binding poses of the agonist for each apo Î¼-ORâGÎ± complexes (or each data point of the 2D free energy landscape) and 2) all-atom binding simulations for each of those complexes. Such simulations would be very insightful but tedious, and as a result separate investigations would be performed in the future to understand the effect of an agonist on the sequence of conformational changes of the receptor and G proteins to form the activated complex. Finally, the mechanistic study of activation with other intracellular binding partners, such as Î²-arrestin, would also provide a clearer picture on biased signaling.
Materials and Methods
Modeling of the End-Point States and Calculation of the Thermodynamic Cycle of Activation.
To generate the end-point states of the thermodynamic cycle of activation of Î¼-OR, Modeller (52) was used. Modeller was adopted as a tool for multitemplate-based modeling. The detailed modeling protocol is described in SI Appendix, section S2.1.1.
The modeled structures were extensively relaxed using GROMACS 2018.4 (53). CHARMM-GUI (54) was employed for the generation of the input files that were needed by GROMACS. The Î¼-OR and G proteins were equilibrated separately. In the case of the Î¼-OR simulations the membranes (composed of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) were included. Ions (Na+ and Clâ) were added to the system for all GROMACS simulations. Details of the MD simulation protocol is provided in SI Appendix, section S2.1. The last few atomic coordinate frames of each production run of the MD simulations were extracted from the trajectory and processed by the Molaris-XG software (46) for computing the free energy using its CG model. CG energy minimization runs were performed, using the Molaris-XG software, on the CG structures, before performing any free energy calculation. In the CG representations the membrane was emulated with a grid of effective atoms. The inclusion of ions might be important in all-atom free energy calculations, but it is not relevant to our CG free energy calculations, which include effective dielectric that accounts well for the energetics. We can also use ionic strength corrections, but this does not seem to be needed in the present case where we deal with differences in free energy and the effect of the ionic strength is cancelled in some respect. A concise description of the theory of the CG model implemented in Molaris-XG software is provided in SI Appendix, section S1.1.
Two-Dimensional Free Energy Landscape Calculations.
The intermediate structures that are required for the 2D free energy landscape calculation were generated using two methods, namely TMD and âshortest path exploration first.â The intermediates of the Î¼-OR and the G proteins were generated separately. Intermediates of Î¼-OR were generated using âshortest path explorationâ and 11 out of 3,000 generated conformations were selected for calculations.
Out of two distinct conformational changes of GiÎ± the intermediates for âoutâ-to-âinâ conformational change were generated using âshortest path exploration,â whereas for the other change (âclosedâ to âopenâ) the TMD method was applied. Twenty-one distinct intermediate configurations of Gi were selected for the free energy calculations.
The structures of all complexes were obtained by combining the intermediate configurations of Î¼-OR and the G proteins using the MatchMaker module (55) [part of the UCSF Chimera software (56)] taking PDB ID code 6DDF (5) as a reference. In order to generate the 2D energy landscape we used CG-based free energy calculation. Detailed descriptions of the TMD and âshortest path exploration firstâ methods are provided in SI Appendix, section S1.3.
Modeling of the Extended Î¼-OR- and Active Î¼-OR-G Protein Complexes.
The structure of Î²2-adrenergic receptor PDB ID code 3SN6 (21) and the sequence of Î¼-OR were used to generate the hypothetical extended active form of Î¼-OR with the help of the SWISS-MODEL (57) software. GROMACS MD simulations were run to equilibrate the modeled structure (details are mentioned in SI Appendix, section S2.1). The structures of all complexes were built by overlaying the Î¼-OR and the mini-G protein structures with respect to PDB ID code 6DDF using the MatchMaker module. The stability of the complexes was estimated by the folding free energy calculations using the CG model.
Solvation and Binding Free Energy Calculation.
We used the PDLD/S-LRA method (58) for the solvation and binding free energy calculations. A relaxation run (0.1 ns) was performed before all solvation and binding free energy calculations. The Molaris-XG software was used for the relaxation runs and for PDLD/S-LRA calculations.
All study data are included in the paper and SI Appendix.
This work was supported by NSF grant MCB 1707167 and NIH grant R35 GM122472. We thank the University of Southern California High Performance Computing and Communication Center, as well as the Extreme Science and Engineering Discovery Environmentâs Comet facility at the San Diego Supercomputing Center, for computational resources. D.M. thanks Dr. Raphael Alhadeff for his guidance and helpful discussion.
- âµ1To whom correspondence may be addressed. Email: .
Author contributions: D.M. and A.W. designed research; D.M. and V.K. performed research; D.M. and A.W. analyzed data; and D.M., V.K., and A.W. wrote the paper.
Reviewers: T.C., Friedrich-Alexander-UniversitÃ¤t Erlangen-NÃ¼rnberg; and D.T., University of Texas at Austin.
The authors declare no competing interest.
This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2013364117/-/DCSupplemental.
Published under the PNAS license.