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

Unraveling the Structural Basis of Biased Agonism in the β2-Adrenergic Receptor Through Molecular Dynamics Simulations

Biased agonism in G protein-coupled receptors is a phenomenon resulting in the selective activation of distinct intracellular signaling pathways by different agonists, which may exhibit bias toward either Gs, Gi, or arrestin-mediated pathways. This study investigates the structural basis of ligand-induced biased agonism within the context of the β2-adrenergic receptor (β2-AR). Atomistic molecular dynamics simulations were conducted for β2-AR complexes with two stereoisomers of methoxynaphtyl fenoterol (MNFen), that is, compounds eliciting qualitatively different cellular responses. The simulations reveal distinct interaction patterns within the binding cavity, dependent on the stereoisomer. These changes propagate to the intracellular parts of the receptor, triggering various structural responses: the dynamic structure of the intracellular regions of the (R,R)-MNFen complex more closely resembles the "Gs-compatible" and "β-arrestin-compatible" conformation of β2-AR, while both stereoisomers maintain structural responses equidistant from the inactive conformation. These findings are confirmed by independent coarse-grained simulations. In the context of deciphered molecular mechanisms, Trp313 plays a pivotal role, altering its orientation upon interactions with (R,R)-MNFen, along with the Lys305-Asp192 ionic bridge. This effect, accompanied by ligand interactions with residues on TM2, increases the strength of interactions within the extracellular region and the binding cavity, resulting in a slightly more open conformation and a minor (by ca. 0.2 nm) increase in the distance between the TM5-TM7, TM1-TM6, TM6-TM7, and TM1-TM5 pairs. On the other hand, an even slighter decrease in the distance between the TM1-TM4 and TM2-TM4 pairs is observed.

Ligand-Induced Biased Activation of GPCRs: Recent Advances and New Directions from In Silico Approaches

G-protein coupled receptors (GPCRs) are the largest family of membrane proteins engaged in transducing signals from the extracellular environment into the cell. GPCR-biased signaling occurs when two different ligands, sharing the same binding site, induce distinct signaling pathways. This selective signaling offers significant potential for the design of safer and more effective drugs. Although its molecular mechanism remains elusive, big efforts are made to try to explain this mechanism using a wide range of methods. Recent advances in computational techniques and AI technology have introduced a variety of simulations and machine learning tools that facilitate the modeling of GPCR signal transmission and the analysis of ligand-induced biased signaling. In this review, we present the current state of in silico approaches to elucidate the structural mechanism of GPCR-biased signaling. This includes molecular dynamics simulations that capture the main interactions causing the bias. We also highlight the major contributions and impacts of transmembrane domains, loops, and mutations in mediating biased signaling. Moreover, we discuss the impact of machine learning models on bias prediction and diffusion-based generative AI to design biased ligands. Ultimately, this review addresses the future directions for studying the biased signaling problem through AI approaches.

Early Events in β2AR Dimer Dynamics Mediated by Activation-Related Microswitches

G-Protein-Coupled Receptors (GPCRs) make up around 3-4% of the human genome and are the targets of one-third of FDA-approved drugs. GPCRs typically exist as monomers but also aggregate to form higher-order oligomers, including dimers. β2AR, a pharmacologically relevant GPCR, is known to be targeted for the treatment of asthma and cardiovascular diseases. The activation of β2AR at the dimer level remains under-explored. In the current study, molecular dynamics (MD) simulations have been performed to understand activation-related structural changes in β2AR at the dimer level. The transition from inactive to active and vice versa has been studied by starting the simulations in the apo, agonist-bound, and inverse agonist-bound β2AR dimers for PDB ID: 2RH1 and PDB ID: 3P0G, respectively. A cumulative total of around 21-μs simulations were performed. Residue-based distances, RMSD, and PCA calculations suggested that either of the one monomer attained activation-related features for the apo and agonist-bound β2AR dimers. The TM5 and TM6 helices within the two monomers were observed to be in significant variation in all the simulations. TM5 bulge and proximity of TM2 and TM7 helices may be contributing to one of the early events in activation. The dimeric interface between TM1 and helix 8 were observed to be well maintained in the apo and agonist-bound simulations. The presence of inverse agonists favored inactive features in both the monomers. These key features of activation known for monomers were observed to have an impact on β2AR dimers, thereby providing an insight into the oligomerization mechanism of GPCRs.

Determinants of Neutral Antagonism and Inverse Agonism in the β2-Adrenergic Receptor

Free-energy profiles for the activation/deactivation of the β2-adrenergic receptor (ADRB2) with neutral antagonist and inverse agonist ligands have been determined with well-tempered multiple-walker (MW) metadynamics simulations. The inverse agonists carazolol and ICI118551 clearly favor single inactive conformational minima in both the binary and ternary ligand-receptor-G-protein complexes, in accord with the inverse-agonist activity of the ligands. The behavior of neutral antagonists is more complex, as they seem also to affect the recruitment of the G-protein. The results are analyzed in terms of the conformational states of the well-known microswitches that have been proposed as indicators of receptor activity.

Biased Signaling in Mutated Variants of β2-Adrenergic Receptor: Insights from Molecular Dynamics Simulations

The molecular basis of receptor bias in G protein-coupled receptors (GPCRs) caused by mutations that preferentially activate specific intracellular transducers over others remains poorly understood. Two experimentally identified biased variants of β2-adrenergic receptors (β2AR), a prototypical GPCR, are a triple mutant (T68F, Y132A, and Y219A) and a single mutant (Y219A); the former bias the receptor toward the β-arrestin pathway by disfavoring G protein engagement, while the latter induces G protein signaling explicitly due to selection against GPCR kinases (GRKs) that phosphorylate the receptor as a prerequisite of β-arrestin binding. Though rigorous characterizations have revealed functional implications of these mutations, the atomistic origin of the observed transducer selectivity is not clear. In this study, we investigated the allosteric mechanism of receptor bias in β2AR using microseconds of all-atom Gaussian accelerated molecular dynamics (GaMD) simulations. Our observations reveal distinct rearrangements in transmembrane helices, intracellular loop 3, and critical residues R1313.50 and Y3267.53 in the conserved motifs D(E)RY and NPxxY for the mutant receptors, leading to their specific transducer interactions. Moreover, partial dissociation of G protein from the receptor core is observed in the simulations of the triple mutant in contrast to the single mutant and wild-type receptor. The reorganization of allosteric communications from the extracellular agonist BI-167107 to the intracellular receptor-transducer interfaces drives the conformational rearrangements responsible for receptor bias in the single and triple mutants. The molecular insights into receptor bias of β2AR presented here could improve the understanding of biased signaling in GPCRs, potentially opening new avenues for designing novel therapeutics with fewer side-effects and superior efficacy.

Coevolution-Driven Method for Efficiently Simulating Conformational Changes in Proteins Reveals Molecular Details of Ligand Effects in the β2AR Receptor

With the advent of AI-powered structure prediction, the scientific community is inching closer to solving protein folding. An unresolved enigma, however, is to accurately, reliably, and deterministically predict alternative conformational states that are crucial for the function of, e.g., transporters, receptors, or ion channels where conformational cycling is innately coupled to protein function. Accurately discovering and exploring all conformational states of membrane proteins has been challenging due to the need to retain atomistic detail while enhancing the sampling along interesting degrees of freedom. The challenges include but are not limited to finding which degrees of freedom are relevant, how to accelerate the sampling along them, and then quantifying the populations of each micro- and macrostate. In this work, we present a methodology that finds relevant degrees of freedom by combining evolution and physics through machine learning and apply it to the conformational sampling of the β2 adrenergic receptor. In addition to predicting new conformations that are beyond the training set, we have computed free energy surfaces associated with the protein's conformational landscape. We then show that the methodology is able to quantitatively predict the effect of an array of ligands on the β2 adrenergic receptor activation through the discovery of new metastable states not present in the training set. Lastly, we also stake out the structural determinants of activation and inactivation pathway signaling through different ligands and compare them to functional experiments to validate our methodology and potentially gain further insights into the activation mechanism of the β2 adrenergic receptor.

Activation/Deactivation Free-Energy Profiles for the β2-Adrenergic Receptor: Ligand Modes of Action

We use enhanced-sampling simulations with an effective collective variable to study the activation of the β2-adrenergic receptor in the presence of ligands with different efficacy. The free-energy profiles are computed for the ligand-free (apo) receptor and binary (apo-receptor + G-protein α-subunit and receptor + ligand) and ternary complexes. The results are not only compatible with available experiments but also allow unprecedented structural insight into the nature of GPCR conformations along the activation pathway and their role in the activation mechanism. In particular, the simulations reveal an unexpected mode of action of partial agonists such as salmeterol and salbutamol that arises already in the binary complex without the G-protein. Specific differences in the polar interactions with residues in TM5, which are required to stabilize an optimal TM6 conformation that facilitates G-protein binding and receptor activation, play a major role in differentiating them from full agonists.

Identification of a β-arrestin-biased negative allosteric modulator for the β2-adrenergic receptor

Catecholamine-stimulated β2-adrenergic receptor (β2AR) signaling via the canonical Gs-adenylyl cyclase-cAMP-PKA pathway regulates numerous physiological functions, including the therapeutic effects of exogenous β-agonists in the treatment of airway disease. β2AR signaling is tightly regulated by GRKs and β-arrestins, which together promote β2AR desensitization and internalization as well as downstream signaling, often antithetical to the canonical pathway. Thus, the ability to bias β2AR signaling toward the Gs pathway while avoiding β-arrestin-mediated effects may provide a strategy to improve the functional consequences of β2AR activation. Since attempts to develop Gs-biased agonists and allosteric modulators for the β2AR have been largely unsuccessful, here we screened small molecule libraries for allosteric modulators that selectively inhibit β-arrestin recruitment to the receptor. This screen identified several compounds that met this profile, and, of these, a difluorophenyl quinazoline (DFPQ) derivative was found to be a selective negative allosteric modulator of β-arrestin recruitment to the β2AR while having no effect on β2AR coupling to Gs. DFPQ effectively inhibits agonist-promoted phosphorylation and internalization of the β2AR and protects against the functional desensitization of β-agonist mediated regulation in cell and tissue models. The effects of DFPQ were also specific to the β2AR with minimal effects on the β1AR. Modeling, mutagenesis, and medicinal chemistry studies support DFPQ derivatives binding to an intracellular membrane-facing region of the β2AR, including residues within transmembrane domains 3 and 4 and intracellular loop 2. DFPQ thus represents a class of biased allosteric modulators that targets an allosteric site of the β2AR.

Elucidation of a dynamic interplay between a beta-2 adrenergic receptor, its agonist, and stimulatory G protein

G protein-coupled receptors (GPCRs) represent the largest group of membrane receptors for transmembrane signal transduction. Ligand-induced activation of GPCRs triggers G protein activation followed by various signaling cascades. Understanding the structural and energetic determinants of ligand binding to GPCRs and GPCRs to G proteins is crucial to the design of pharmacological treatments targeting specific conformations of these proteins to precisely control their signaling properties. In this study, we focused on interactions of a prototypical GPCR, beta-2 adrenergic receptor (β2AR), with its endogenous agonist, norepinephrine (NE), and the stimulatory G protein (Gs). Using molecular dynamics (MD) simulations, we demonstrated the stabilization of cationic NE, NE(+), binding to β2AR by Gs protein recruitment, in line with experimental observations. We also captured the partial dissociation of the ligand from β2AR and the conformational interconversions of Gs between closed and open conformations in the NE(+)-β2AR-Gs ternary complex while it is still bound to the receptor. The variation of NE(+) binding poses was found to alter Gs α subunit (Gsα) conformational transitions. Our simulations showed that the interdomain movement and the stacking of Gsα α1 and α5 helices are significant for increasing the distance between the Gsα and β2AR, which may indicate a partial dissociation of Gsα The distance increase commences when Gsα is predominantly in an open state and can be triggered by the intracellular loop 3 (ICL3) of β2AR interacting with Gsα, causing conformational changes of the α5 helix. Our results help explain molecular mechanisms of ligand and GPCR-mediated modulation of G protein activation.

Mechanistic insight into the impact of a bivalent ligand on the structure and dynamics of a GPCR oligomer

Development of effective bivalent ligands has become the focus of intensive research toward modulation of G protein-coupled receptor (GPCR) oligomers, particularly in the field of GPCR pharmacology. Experimental studies have shown that they increased binding affinity and signaling potency compared to their monovalent counterparts, yet underlying molecular mechanism remains elusive. To address this, we performed accelerated molecular dynamics simulations on bivalent-ligand bound Adenosine 2A receptor (A2AR) dimer in the context of a modeled tetramer, which consists of A2AR and dopamine 2 receptor (D2R) homodimers and their cognate G proteins. Our results demonstrate that bivalent ligand impacted interactions between pharmacophore groups and ligand binding residues, thus modulating allosteric communication network and water channel formed within the receptor. Moreover, it also strengthens contacts between receptor and G protein, by modulating the volume of ligand binding pocket and intracellular domain of the receptor. Importantly, we showed that impact evoked by the bivalent ligand on A2AR dimer was also transmitted to apo D2R, which is part of the neighboring D2R dimer. To the best of our knowledge, this is the first study that provides a mechanistic insight into the impact of a bivalent ligand on dynamics of a GPCR oligomer. Consequently, this will pave the way for development of effective ligands for modulation of GPCR oligomers and hence treatment of crucial diseases such as Parkinson's disease and cancer.

Molecular Mechanisms of Diverse Activation Stimulated by Different Biased Agonists for the β2-Adrenergic Receptor

β2AR is an important drug target protein involving many diseases. Biased drugs induce specific signaling and provide additional clinical utility to optimize β2AR-based therapies. However, the biased signaling mechanism has not been elucidated. Motivated by the issue, we chose four agonists with divergent bias (balanced agonist, G-protein-biased agonist, and β-arrestin-biased agonists) and utilized Gaussian accelerated molecular dynamics simulation coupled with a dynamic network to probe the molecular mechanisms of distinct biased activation induced by the structural differences between the four agonists. Our simulations reveal that the G-protein-biased agonist induces an open conformation with the outward shifts of TM6 and TM7 for the intracellular domain, which will be beneficial to couple G protein. In contrast, the β-arrestin-biased agonists regulate an occluded conformation with a slightly outward movement of TM6 and an inward shift of TM7, which should favor β-arrestin signaling. The balanced agonist does not induce an observable outward shift for TM6 but, along with a slight tilt for TM7, leads to an inactive-like conformation. In addition, our results reveal the first time that ICL3 presents specific conformations with different agonists. The G-protein-biased agonist drives ICL3 to open so that the G protein-binding pocket can be available, while the β-arrestin-biased agonists induce ICL3 to form a closed conformation with a stable local α-helix. MM/PBSA analysis further reveals that the hydroxyl groups in the resorcinol of the G-protein-biased agonist form strong interactions with Y5.38 and S5.42, thus preventing tilting of the TM5 extracellular end. The catechol of the balanced agonist and the β-arrestin-biased ones induces the rearrangement of two hydrophobic residues F6.52 and W6.48. However, different from the balanced agonist, the ethyl substituent of β-arrestin-biased agonists forms additional hydrophobic interactions with W6.48 and F6.51 after the rearrangement, which should contribute to the β-arrestin bias. The shortest pathway analysis further reveals that the three residues Y7.43, N7.45, and N7.49 are crucial for allosterically regulating G-protein-biased signaling, while the two residues W6.48 and F6.44 make an important contribution to regulate β-arrestin-biased signaling. For the balanced agonist NE, the allosteric regulation pathway simultaneously involves the residue associated with G-protein-biased signaling like S5.46 and the residues related to β-arrestin-biased signaling like W6.48 and F6.44, thus producing unbiased signaling. The observations could advance our understanding of the biased activation mechanism on class A GPCRs and provide a useful guideline for the design of biased drugs.

Identification of ligand-specific G protein-coupled receptor states and prediction of downstream efficacy via data-driven modeling

Ligand binding stabilizes different G protein-coupled receptor states via a complex allosteric process that is not completely understood. Here, we have derived free energy landscapes describing activation of the β2 adrenergic receptor bound to ligands with different efficacy profiles using enhanced sampling molecular dynamics simulations. These reveal shifts toward active-like states at the Gprotein-binding site for receptors bound to partial and full agonists, and that the ligands modulate the conformational ensemble of the receptor by tuning protein microswitches. We indeed find an excellent correlation between the conformation of the microswitches close to the ligand binding site and in the transmembrane region and experimentally reported cyclic adenosine monophosphate signaling responses. Dimensionality reduction further reveals the similarity between the unique conformational states induced by different ligands, and examining the output of classifiers highlights two distant hotspots governing agonism on transmembrane helices 5 and 7.

Deciphering collaborative sidechain motions in proteins during molecular dynamics simulations

The dynamic structure of proteins is essential for their functions and may include large conformational transitions which can be studied by molecular dynamics (MD) simulations. However, details of these transitions are difficult to automatically track. To facilitate their analysis, we developed two scores of correlation between sidechain dihedral angles. The CIRCULAR and OMES scores are computed from, respectively, dihedral angle values and rotamer distributions. As a case study, we applied our methods to an activation-like transition of the chemokine receptor CXCR4, observed during accelerated MD simulations. The principal component analysis of the correlation matrices was consistent with the networking structure of the top ranking pairs. Both scores identify a set of residues whose “collaborative” sidechain rotamerization immediately preceded or accompanied the conformational transition of CXCR4. Detailed analysis of the sequential order of these rotamerizations suggests that an allosteric mechanism, involving the outward motion of an asparagine residue in transmembrane helix 3, might be a prerequisite to the large scale conformational transition of CXCR4. This case study provides the proof-of-concept that the correlation methods developed here are valuable exploratory techniques to help decipher complex reactional pathways.

Energy Landscapes Reveal Agonist Control of G Protein-Coupled Receptor Activation via Microswitches

Agonist binding to G protein-coupled receptors (GPCRs) leads to conformational changes in the transmembrane region that activate cytosolic signaling pathways. Although high-resolution structures of different receptor states are available, atomistic details of allosteric signaling across the membrane remain elusive. We calculated free energy landscapes of β₂ adrenergic receptor activation using atomistic molecular dynamics simulations in an optimized string of swarms framework, which shed new light on how microswitches govern the equilibrium between conformational states. Contraction of the extracellular binding site in the presence of the agonist BI-167107 is obligatorily coupled to conformational changes in a connector motif located in the core of the transmembrane region. The connector is probabilistically coupled to the conformation of the intracellular region. An active connector promotes desolvation of a buried cavity, a twist of the conserved NPxxY motif, and an interaction between two conserved tyrosines in transmembrane helices 5 and 7 (Y–Y motif), which lead to a larger population of active-like states at the G protein binding site. This coupling is augmented by protonation of the strongly conserved Asp79^2.50. The agonist binding site hence communicates with the intracellular region via a cascade of locally connected microswitches. Characterization of these can be used to understand how ligands stabilize distinct receptor states and contribute to development drugs with specific signaling properties. The developed simulation protocol can likely be transferred to other class A GPCRs.

Quantum Mechanical and Molecular Mechanics Modeling of Membrane-Embedded Rhodopsins

Computational chemistry provides versatile methods for studying the properties and functioning of biological systems at different levels of precision and at different time scales. The aim of this article is to review the computational methodologies that are applicable to rhodopsins as archetypes for photoactive membrane proteins that are of great importance both in nature and in modern technologies. For each class of computational techniques, from methods that use quantum mechanics for simulating rhodopsin photophysics to less-accurate coarse-grained methodologies used for long-scale protein dynamics, we consider possible applications and the main directions for improvement.

Cell-based assays and molecular dynamics analysis of a boron-containing agonist with different profiles of binding to human and guinea pig beta2 adrenoceptors

The design of beta2 adrenoceptor (β2AR) agonists is attractive because of their wide-ranging applications in medicine, and the details of agonist interactions with β2AR are interesting because it is considered a prototype for G-protein coupled receptors. Preclinical studies for agonist development have involved biological assays with guinea pigs due to a similar physiology to humans. Boron-containing Albuterol derivatives (BCADs) designed as bronchodilators have improved potency and efficacy compared with their boron-free precursor on guinea pig β2ARs (gpβ2ARs), and two of the BCADs (BR-AEA and boronterol) conserve these features on cells expressing human β2ARs (hβ2ARs). The aim of this study was to test the BCAD Politerol on gpβ2ARs and hβ2ARs in vitro and in silico. Politerol displayed higher potency and efficacy on gpβ2AR than on hβ2AR in experimental assays, possible explanations are provided based on molecular modeling, and molecular dynamics simulations of about 0.25 µs were performed for the free and bound states adding up to 2 µs in total. There were slight differences, particularly in the role of the boron atom, in the interactions of Politerol with gpβ2ARs and hβ2ARs, affecting movements of transmembrane domains 5-7, known to be pivotal in receptor activation. These findings could be instrumental in the design of compounds selective for hβ2ARs.

Structural insights into positive and negative allosteric regulation of a G protein-coupled receptor through protein-lipid interactions

Lipids are becoming known as essential allosteric modulators of G protein-coupled receptor (GPCRs). However, how they exert their effects on GPCR conformation at the atomic level is still unclear. In light of recent experimental data, we have performed several long-timescale molecular dynamics (MD) simulations, totalling 24 μs, to rigorously map allosteric modulation and conformational changes in the β₂ adrenergic receptor (β2AR) that occur as a result of interactions with three different phospholipids. In particular, we identify different sequential mechanisms behind receptor activation and deactivation, respectively, mediated by specific lipid interactions with key receptor regions. We show that net negatively charged lipids stabilize an active-like state of β2AR that is able to dock G_sα protein. Clustering of anionic lipids around the receptor with local distortion of membrane thickness is also apparent. On the other hand, net-neutral zwitterionic lipids inactivate the receptor, generating either fully inactive or intermediate states, with kinetics depending on lipid headgroup charge distribution and hydrophobicity. These chemical differences alter membrane thickness and density, which differentially destabilize the β2AR active state through lateral compression effects.

Exploring the mechanism of F282L mutation-caused constitutive activity of GPCR by a computational study

G-protein-coupled receptors (GPCRs) are important drug targets and generally activated by ligands. However, some experiments found that GPCRs also give rise to constitutive activity through some mutations (viz., CAM), which are usually associated with different kinds of diseases. However, the mechanisms of CAMs and their roles in interactions with drug-ligands are unclear in experiments. Herein, we used microsecond molecular dynamics simulations to study the effect of one important F282L mutation on β2AR in order to address the questions above. With the aid of principle component and correlation analysis, our results revealed that the F282L mutation could increase the instability of the overall structure, increase the dramatic fluctuations of NPxxY and extracellular loops, and decrease restraint of the helices through weakening interhelical H-bonding and correlations between residues, which could partly contribute to the constitutive activity reported by the experiments. The observations from the protein structure network (PSN) analysis indicate that the mutant exhibits less information flow than the wild β2AR and weakens the role of TM5 and TM6 in the signal transmission, but it enhances the impact of TM3 on the orthosteric pathway and TM4 on the allosteric one. In addition, the results from the virtual screening reveal that the mutant prefers to select agonists rather than antagonists, similar to the active state but opposite of the inactive state, further confirming that the F282L mutation advances the activation of β2AR. Our observations provide valuable information for understanding the mechanism of the mutation-caused constitutive activity of GPCR and related drug-design.

Evaluation of variability in high-resolution protein structures by global distance scoring

Systematic analysis of the statistical and dynamical properties of proteins is critical to understanding cellular events. Extraction of biologically relevant information from a set of high-resolution structures is important because it can provide mechanistic details behind the functional properties of protein families, enabling rational comparison between families. Most of the current structural comparisons are pairwise-based, which hampers the global analysis of increasing contents in the Protein Data Bank. Additionally, pairing of protein structures introduces uncertainty with respect to reproducibility because it frequently accompanies other settings for superimposition. This study introduces intramolecular distance scoring for the global analysis of proteins, for each of which at least several high-resolution structures are available. As a pilot study, we have tested 300 human proteins and showed that the method is comprehensively used to overview advances in each protein and protein family at the atomic level. This method, together with the interpretation of the model calculations, provide new criteria for understanding specific structural variation in a protein, enabling global comparison of the variability in proteins from different species.

Structural Features and Ligand Selectivity for 10 Intermediates in the Activation Process of β2-Adrenergic Receptor

It has already been suggested by researchers that there should be multiple intermediate states in the activation process for G-protein-coupled receptors (GPCRs). However, the intermediate states are very short-lived and hardly captured by the experiments, leading to very limited understanding of their structural features and drug efficacies. In this work, a novel joint strategy of targeted molecular dynamics simulation, conventional molecular dynamics simulation, and virtual screening is developed to address the problems. The results from 10 intermediate conformations obtained from the work reveal that the ligand pocket is very unstable and fluctuates between the inactive state and the active one in the case of ligand-free, in particular for ECL2 as a gate-keeper of the ligand-binding. The ligand-binding site could be stable in the active state with a small volume and a completely closed ECL2, only when the G-protein-binding region is fully activated. In addition, the activations of the ligand-binding pocket and G-protein-binding site are relatively independent and exhibit a loose allosteric coupling, which contributes to the existence of multiple intermediate conformations. Interestingly, the screening performance of the agonists does not increase on increasing the overall activity of the intermediate state, but is dependent on the activated extent of the ligand pocket. The receptor is prone to bind the agonist when closing ECL2 and reducing the ligand-binding pocket volume, whereas it is more favorable for binding the antagonist when opening ECL2 and increasing the pocket volume. These observations added to previous studies could help us better understand the activation mechanism of GPCRs and provide valuable information for drug design.

Improving virtual screening of G protein-coupled receptors via ligand-directed modeling

G protein-coupled receptors (GPCRs) play crucial roles in cell physiology and pathophysiology. There is increasing interest in using structural information for virtual screening (VS) of libraries and for structure-based drug design to identify novel agonist or antagonist leads. However, the sparse availability of experimentally determined GPCR/ligand complex structures with diverse ligands impedes the application of structure-based drug design (SBDD) programs directed to identifying new molecules with a select pharmacology. In this study, we apply ligand-directed modeling (LDM) to available GPCR X-ray structures to improve VS performance and selectivity towards molecules of specific pharmacological profile. The described method refines a GPCR binding pocket conformation using a single known ligand for that GPCR. The LDM method is a computationally efficient, iterative workflow consisting of protein sampling and ligand docking. We developed an extensive benchmark comparing LDM-refined binding pockets to GPCR X-ray crystal structures across seven different GPCRs bound to a range of ligands of different chemotypes and pharmacological profiles. LDM-refined models showed improvement in VS performance over origin X-ray crystal structures in 21 out of 24 cases. In all cases, the LDM-refined models had superior performance in enriching for the chemotype of the refinement ligand. This likely contributes to the LDM success in all cases of inhibitor-bound to agonist-bound binding pocket refinement, a key task for GPCR SBDD programs. Indeed, agonist ligands are required for a plethora of GPCRs for therapeutic intervention, however GPCR X-ray structures are mostly restricted to their inactive inhibitor-bound state.

G protein-coupled receptors (GPCRs) are a major target for drug discovery. These receptors are highly dynamic membrane proteins, and have had limited tractability using with biophysical screens that are widely adopted for globular protein targets. Thus, structure-based virtual screening (SBVS) holds great promise as a complement to physical screening for rational design of novel drugs. Indeed, the increasing number of atomic-detail GPCR X-ray crystal structures has coincided with an increase in prospective SBVS studies that have identified novel compounds. However, experimentally solved GPCR structures do not meet the full demand for SBVS, as the GPCR structural landscape is incomplete, lacking both in coverage of available GPCRs, and diversity in both receptor conformations and the chemistry of co-crystalised ligands. Here we present a novel computational GPCR binding pocket refinement method that can generate predictive GPCR/ligand complexes with improved SBVS performance. This ligand-directed modeling workflow uses parallel processing and efficient algorithms to search the GPCR/ligand conformational space faster and more efficiently than the widely used protein refinement method molecular dynamics. In this study, the resulting models are evaluated both structurally, and in retrospective SBVS. We demonstrate improved performance of refined models over their starting structures in the majority of our test cases.

Identification of a Different Agonist-Binding Site and Activation Mechanism of the Human P2Y1 Receptor

The human P2Y1 receptor (P2Y1R) is a purinergic G-protein-coupled receptor (GPCR) that functions as a receptor for adenosine 5'-diphosphate (ADP). An antagonist of P2Y1R might potentially have antithrombotic effects, whereas agonists might serve as antidiabetic agents. On the basis of the antagonist-bound MRS2500-P2Y1R crystal structure, we constructed computational models of apo-P2Y1R and the agonist-receptor complex 2MeSADP-P2Y1R. We then performed conventional molecular dynamics (cMD) and accelerated molecular dynamics (aMD) simulations to study the conformational dynamics after binding with agonist/antagonist as well as the P2Y1R activation mechanism. We identified a new agonist-binding site of P2Y1R that is consistent with previous mutagenesis data. This new site is deeper than those of the agonist ADP in the recently simulated ADP-P2Y1R structure and the antagonist MRS2500 in the MRS2500-P2Y1R crystal structure. During P2Y1R activation, the cytoplasmic end of helix VI shifts outward 9.1 Å, the Ser1463.47-Tyr2375.58 hydrogen bond breaks, a Tyr2375.58-Val2626.37 hydrogen bond forms, and the conformation of the χ1 rotamer of Phe2696.44 changes from parallel to perpendicular to helix VI. The apo-P2Y1R system and the MRS2500-P2Y1R system remain inactive. The newly identified agonist binding site and activation mechanism revealed in this study may aid in the design of P2Y1R antagonists/agonists as antithrombotic/antidiabetic agents, respectively.

Activation-Induced Conformational Changes of Dopamine D3 Receptor Promote the Formation of the Internal Water Channel

The atomic-level dopamine activation mechanism for transmitting extracellular ligand binding events through transmembrane helices to the cytoplasmic G protein remains unclear. In the present study, the complete dopamine D3 receptor (D3R), with a homology-modeled N-terminus, was constructed to dock different ligands to simulate conformational alterations in the receptor’s active and inactive forms during microsecond-timescale molecular dynamic simulations. In agonist-bound systems, the D3R N-terminus formed a “lid-like” structure and lay flat on the binding site opening, whereas in antagonist and inverse agonist-bound systems, the N-terminus exposed the binding cavity. Receptor activation was characterized using the different molecular switch residue distances, and G protein-binding site volumes. A continuous water pathway was observed only in the dopamine-G_αi-bound system. In the inactive D3Rs, water entry was hindered by the hydrophobic layers. Finally, a complete activation mechanism of D3R was proposed. Upon agonist binding, the “lid-like” conformation of the N-terminus induces a series of molecular switches to increase the volume of the D3R cytoplasmic binding part for G protein association. Meanwhile, water enters the transmembrane region inducing molecular switches to assist in opening the hydrophobic layers to form a continuous water channel, which is crucial for maintaining a fully active conformation for signal transduction.

Study of the mechanism of protonated histidine-induced conformational changes in the Zika virus dimeric envelope protein using accelerated molecular dynamic simulations

The Zika virus has drawn worldwide attention because of the epidemic diseases it causes. It is a flavivirus that has an icosahedral protein shell constituted by an envelope glycoprotein (E-protein) and membrane protein (M-protein) in the mature virion. The multistep process of membrane fusion to infect the host cell is pH-induced. To understand the mechanism of the conformational changes in the (E-M)₂ protein homodimer embedded in the membrane, two 200-ns accelerated dynamic simulations were performed under different pH conditions. The low pH condition weakens the interactions and correlations in both E-protein monomers and in the E-M heterodimer. The highly conserved residues, His249, His288, His323 and His446, are protonated under low pH conditions and play key roles in driving the fusion process. The analysis and discussion in this study may provide some insight into the molecular mechanism of Zika virus infection.

Functional Characterization of a Novel Series of Biased Signaling Dopamine D3 Receptor Agonists

Dopamine receptors play an integral role in controlling brain physiology. Importantly, subtype selective agonists and antagonists of dopamine receptors with biased signaling properties have been successful in treating psychiatric disorders with a low incidence of side effects. To this end, we recently designed and developed SK609, a dopamine D3 receptor (D3R) selective agonist that has atypical signaling properties. SK609 has shown efficacy in reversing akinesia and reducing L-dopa-induced dyskinesia in a hemiparkinsonian rats. In the current study, we demonstrate that SK609 has high selectivity for D3R with no binding affinity on D2R high- or low-affinity state when tested at a concentration of 10 μM. In addition, SK609 and its analogues do not induce desensitization of D3R as determined by repeated agonist treatment response in phosphorylation of ERK1/2 functional assay. Most significantly, SK609 and its analogues preferentially signal through the G-protein-dependent pathway and do not recruit β-arrestin-2, suggesting a functional bias toward the G-protein-dependent pathway. Structure-activity relationship (SAR) studies using analogues of SK609 demonstrate that the molecules bind at the orthosteric site by maintaining the conserved salt bridge interactions with aspartate 110 on transmembrane 3 and aryl interactions with histidine 349 on transmembrane 6, in addition to several hydrophobic interactions with residues from transmembranes 5 and 6. The compounds follow a strict SAR with reference to the three pharmacophore elements: substituted phenyl ring, length of the linker connecting phenyl ring and amine group, and orientation and hydrophobic branching groups at the amine among SK609 analogues for efficacy and functional selectivity. These features of SK609 and the analogues suggest that biased signaling is an inherent property of this series of molecules.

Coarse-Grained Prediction of Peptide Binding to G-Protein Coupled Receptors

In this study, we used the Martini Coarse-Grained model with no applied restraints to predict the binding mode of some peptides to G-Protein Coupled Receptors (GPCRs). Both the Neurotensin-1 and the chemokine CXCR4 receptors were used as test cases. Their ligands, NTS_8-13 and CVX15 peptides, respectively, were initially positioned in the surrounding water box. Using a protocol based on Replica Exchange Molecular Dynamics (REMD), both opening of the receptors and entry of the peptides into their dedicated pockets were observed on the μs time-scale. After clustering, the most statistically representative orientations were closely related to the X-ray structures of reference, sharing both RMSD lower than 3 Å and most of the native contacts. These results demonstrate that such a model, that does not require access to tremendous computational facilities, can be helpful in predicting peptide binding to GPCRs as well as some of the receptor's conformational changes required for this key step. We also discuss how such an approach can now help to predict, de novo, the interactions of GPCRs with other intra- or extracellular peptide/protein partners.

Molecular Mechanism of Biased Ligand Conformational Changes in CC Chemokine Receptor 7

Biased ligand binding to G protein-coupled receptors enables functional selectivity of intracellular effectors to mediate cellular function. Despite the significant advances made in characterizing the conformational states (transmembrane helical arrangements) capable of discriminating between G protein and arrestin binding, the role of the ligand in stabilizing such conformations remains unclear. To address this issue, we simulate microsecond dynamics of CC chemokine receptor 7 (CCR7) bound to its native biased ligands, CCL19 and CCL21, and detect a series of molecular switches that are mediated by various ligand-induced allosteric events. These molecular switches involve three tyrosine residues (Y112(3.32), Y255(6.51), and Y288(7.39)), three phenylalanine residues (F116(3.36), F208(5.47), and F248(6.44)), and a polar interaction between Q252(6.48) and R294(7.45) in the transmembrane domain of CCR7. Conformational changes within these switches, particularly hydrogen bond formation between Y112(3.32) and Y255(6.51), lead to global helical movements in the receptor's transmembrane helices and contribute to the transitioning of the receptor to distinct states. Ligand-induced helical movements in the receptor highlight the ability of biased ligands to stabilize the receptor in different states through a dynamic network of allosteric events.

G-Protein/β-Arrestin-Linked Fluctuating Network of G-Protein-Coupled Receptors for Predicting Drug Efficacy and Bias Using Short-Term Molecular Dynamics Simulation

The efficacy and bias of signal transduction induced by a drug at a target protein are closely associated with the benefits and side effects of the drug. In particular, partial agonist activity and G-protein/β-arrestin-biased agonist activity for the G-protein-coupled receptor (GPCR) family, the family with the most target proteins of launched drugs, are key issues in drug discovery. However, designing GPCR drugs with appropriate efficacy and bias is challenging because the dynamic mechanism of signal transduction induced by ligand—receptor interactions is complicated. Here, we identified the G-protein/β-arrestin-linked fluctuating network, which initiates large-scale conformational changes, using sub-microsecond molecular dynamics (MD) simulations of the β₂-adrenergic receptor (β₂AR) with a diverse collection of ligands and correlation analysis of their G protein/β-arrestin efficacy. The G-protein-linked fluctuating network extends from the ligand-binding site to the G-protein-binding site through the connector region, and the β-arrestin-linked fluctuating network consists of the NPxxY motif and adjacent regions. We confirmed that the averaged values of fluctuation in the fluctuating network detected are good quantitative indexes for explaining G protein/β-arrestin efficacy. These results indicate that short-term MD simulation is a practical method to predict the efficacy and bias of any compound for GPCRs.

Activation and conformational dynamics of a class B G-protein-coupled glucagon receptor

The human glucagon receptor (GCGR) is a class B G-protein-coupled receptor (GPCR). The GCGR can be activated by glucagon and regulates the release of glucose. The GCGR has been proposed to be an important drug target for type 2 diabetes. Based on the structural model of a full-length glucagon-bound GCGR (glu-GCGR), we performed accelerated molecular dynamics (aMD) simulations, potential of mean force (PMF) calculations, cross-correlation analysis and community network analysis to study the activation mechanism and the conformational dynamics during the activation process. The PMF map depicts three different conformational states of the GCGR: the inactive, intermediate and active states. The activation of the GCGR is characterized by the outward movement of the intracellular side of helix VI. In the active state of the GCGR, the Arg173(2.46)-Ser350(6.41) and Glu245(3.50)-Thr351(6.42) hydrogen bonds break, and the χ1 rotamer of Phe322(5.54) changes from perpendicular to parallel to helix VI. The binding of the agonist glucagon decreases the correlated motions of the extracellular loops (ELCs) and the helices around the glucagon-binding site. During the activation of the GCGR, the connections between the intracellular sides of helices become weaker, and the connections between glucagon and ECLs and the extracellular sides of helices become stronger. These facilitate G-protein coupling on the intracellular side and glucagon binding on the extracellular side, and stabilize the GCGR in the active state. We expect that this study can provide useful information on the activation mechanism of the GCGR and facilitate the future design of GCGR inhibitors.

Understanding the molecular basis of agonist/antagonist mechanism of GPER1/GPR30 through structural and energetic analyses

The G-protein coupled receptors (GPCRs) represent the largest superfamily of membrane proteins in charge to pass the cell signaling after binding with their cognate ligands to the cell interior. In breast cancer, a GPCR named GPER1 plays a key role in the process of growth and the proliferation of cancer cells. In a previous study, theoretical methods were applied to construct a model of GPER1, which later was submitted to molecular dynamics (MD) simulations to perform a docking calculation. Based on this preceding work, it is known that GPER1 is sensitive to structural differences in its binding site. However, due to the nature of that past study, conformational changes linked to the ligand binding were not observed. Therefore, in this study, in order to explore the conformational changes coupled to the agonist/antagonist binding, MD simulations of about 0.25μs were performed for the free and bound states, summarizing 0.75μs of MD simulation in total. For the bound states, one agonist (G-1) and antagonist (G-15) were chosen since is widely known that these two molecules cause an impact on GPER1 mobility. Based on the conformational ensemble generated through MD simulations, we found that despite G-1 and G-15 being stabilized by similar map of residues, the structural differences between both ligands impact the hydrogen bond pattern not only at the GPER1 binding site but also along the seven-helix bundle, causing significant differences in the conformational mobility along the extracellular and cytoplasmic domain, and to a lesser degree in the curvatures of helix 2, helix 3 and helix 7 between the free and bound states, which is in agreement with reported literature, and might be linked to microscopic characteristics of the activated-inactivated transition. Furthermore, binding free energy calculations using the MM/GBSA method for the bound states, followed by an alanine scanning analysis allowed us to identify some important residues for the complex stabilization.

Molecular Dynamics Simulations of the Human Glucose Transporter GLUT1

Glucose transporters (GLUTs) provide a pathway for glucose transport across membranes. Human GLUTs are implicated in devastating diseases such as heart disease, hyper- and hypo-glycemia, type 2 diabetes and cancer. The human GLUT1 has been recently crystalized in the inward-facing open conformation. However, there is no other structural information for other conformations. The X-ray structures of E. coli Xylose permease (XylE), a glucose transporter homolog, are available in multiple conformations with and without the substrates D-xylose and D-glucose. XylE has high sequence homology to human GLUT1 and key residues in the sugar-binding pocket are conserved. Here we construct a homology model for human GLUT1 based on the available XylE crystal structure in the partially occluded outward-facing conformation. A long unbiased all atom molecular dynamics simulation starting from the model can capture a new fully opened outward-facing conformation. Our investigation of molecular interactions at the interface between the transmembrane (TM) domains and the intracellular helices (ICH) domain in the outward- and inward-facing conformation supports that the ICH domain likely stabilizes the outward-facing conformation in GLUT1. Furthermore, inducing a conformational transition, our simulations manifest a global asymmetric rocker switch motion and detailed molecular interactions between the substrate and residues through the water-filled selective pore along a pathway from the extracellular to the intracellular side. The results presented here are consistent with previously published biochemical, mutagenesis and functional studies. Together, this study shed light on the structure and functional relationships of GLUT1 in multiple conformational states.

Functional loop dynamics of the streptavidin-biotin complex

Accelerated molecular dynamics (aMD) simulation is employed to study the functional dynamics of the flexible loop(3-4) in the strong-binding streptavidin-biotin complex system. Conventional molecular (cMD) simulation is also performed for comparison. The present study reveals the following important properties of the loop dynamics: (1) The transition of loop(3-4) from open to closed state is observed in 200 ns aMD simulation. (2) In the absence of biotin binding, the open-state streptavidin is more stable, which is consistent with experimental evidences. The free energy (ΔG) difference is about 5 kcal/mol between two states. But with biotin binding, the closed state is more stable due to electrostatic and hydrophobic interactions between the loop(3-4) and biotin. (3) The closure of loop(3-4) is concerted to the stable binding of biotin to streptavidin. When the loop(3-4) is in its open-state, biotin moves out of the binding pocket, indicating that the interactions between the loop(3-4) and biotin are essential in trapping biotin in the binding pocket. (4) In the tetrameric streptavidin system, the conformational change of the loop(3-4) in each monomer is independent of each other. That is, there is no cooperative binding for biotin bound to the four subunits of the tetramer.

Insights into a defined secondary binding region on β-adrenoceptors and putative roles in ligand binding and drug design

Putative roles of a secondary binding region shared among beta-adrenoceptors.

Computational structural biology of opioid receptors

The publication of high-resolution structures for all of the opioid receptor subfamilies has unveiled exciting opportunities for mechanistic insight into the molecular mechanisms underlying the biology of nociception, reward, and higher cognitive functions, as well as promises for progress in several clinical areas such as pain management, physiological dependence, addiction, and mood disorders. To turn this promise into novel and improved therapeutic entities, however, this information needs to be supplemented with research strategies that explore the dynamic behavior of the proteins and their interactions with other receptors and ligands in their physiological environment.Here we describe state-of-the-art molecular dynamics computational protocols, based on all-atom and coarse-grained modeling techniques, designed to estimate crucial thermodynamic and kinetic parameters describing the binding of small-molecule ligands and the formation of supramolecular complexes.

Structure-bias relationships for fenoterol stereoisomers in six molecular and cellular assays at the β2-adrenoceptor

Functional selectivity is well established as an underlying concept of ligand-specific signaling via G protein-coupled receptors (GPCRs). Functionally, selective drugs could show greater therapeutic efficacy and fewer adverse effects. Dual coupling of the β2-adrenoceptor (β2AR) triggers a signal transduction via Gsα and Giα proteins. Here, we examined 12 fenoterol stereoisomers in six molecular and cellular assays. Using β2AR-Gsα and β2AR-Giα fusion proteins, (R,S')- and (S,S')-isomers of 4'-methoxy-1-naphthyl-fenoterol were identified as biased ligands with preference for Gs. G protein-independent signaling via β-arrestin-2 was disfavored by these ligands. Isolated human neutrophils constituted an ex vivo model of β2AR signaling and demonstrated functional selectivity through the dissociation of cAMP accumulation and the inhibition of formyl peptide-stimulated production of reactive oxygen species. Ligand bias was calculated using an operational model of agonism and revealed that the fenoterol scaffold constitutes a promising lead structure for the development of Gs-biased β2AR agonists.

Modeling G protein-coupled receptors in complex with biased agonists

The biological response to the activation of G protein-coupled receptors (GPCRs) typically originates from the simultaneous modulation of various signaling pathways that lead to distinct biological consequences. Hence, 'biased agonists' (i.e., compounds that selectively activate one of the pathways while blocking the others) are highly sought-after molecules to provide fine-tuned pharmacological interventions. This review describes strategies that can be deployed to model the conformation of GPCRs in complex with ligands endowed with specific signaling profiles useful for the generation of hypotheses on the structural requirements for the activation of different signaling pathways or for rational computer-aided ligand discovery campaigns. In particular, it focuses on strategies potentially applicable to model the global or local conformational states of GPCRs stabilized by specific ligands.

NbIT--a new information theory-based analysis of allosteric mechanisms reveals residues that underlie function in the leucine transporter LeuT

Complex networks of interacting residues and microdomains in the structures of biomolecular systems underlie the reliable propagation of information from an input signal, such as the concentration of a ligand, to sites that generate the appropriate output signal, such as enzymatic activity. This information transduction often carries the signal across relatively large distances at the molecular scale in a form of allostery that is essential for the physiological functions performed by biomolecules. While allosteric behaviors have been documented from experiments and computation, the mechanism of this form of allostery proved difficult to identify at the molecular level. Here, we introduce a novel analysis framework, called N-body Information Theory (NbIT) analysis, which is based on information theory and uses measures of configurational entropy in a biomolecular system to identify microdomains and individual residues that act as (i)-channels for long-distance information sharing between functional sites, and (ii)-coordinators that organize dynamics within functional sites. Application of the new method to molecular dynamics (MD) trajectories of the occluded state of the bacterial leucine transporter LeuT identifies a channel of allosteric coupling between the functionally important intracellular gate and the substrate binding sites known to modulate it. NbIT analysis is shown also to differentiate residues involved primarily in stabilizing the functional sites, from those that contribute to allosteric couplings between sites. NbIT analysis of MD data thus reveals rigorous mechanistic elements of allostery underlying the dynamics of biomolecular systems.

We developed the new information theory-based analysis framework presented here, NbIT analysis, for the study of allosteric mechanisms in biomolecular systems from Molecular Dynamics trajectories. The illustrative application of NbIT to the analysis of the occluded state in the bacterial transporter LeuT, produced a quantitative representation of the allosteric behavior, and identified intramolecular channels that enable the long-distance information transmission. Our findings, identifying the roles of specific residues in the communication of the allosteric information, were validated by the recognition of residues that have been previously shown to play functional roles in this very well studied system. In addition, we show that application of NbIT analysis leads to the discrimination of functional roles by differentiating between residues that are essential to the dynamics within functional sites (e.g., the substrate binding sites), and residues whose role is to communicate between such functional sites. These results demonstrate that the information theoretical analysis presented here is a powerful tool for quantifying complex allosteric behavior in biomolecular systems and for identifying the crucial components underlying those behaviors.