Advances in enhanced sampling simulations of G-protein-coupled receptors

1 Introduction

G-protein-coupled receptors (GPCRs) constitute the largest and most pharmacologically rich family of human membrane proteins, with more than 800 members and roughly one-third of approved drugs acting upon them (Fredriksson et al., 2003; Hauser et al., 2017). Their function rests on the ability to interconvert between inactive, intermediate, and fully active conformations in response to chemically diverse ligands, a plasticity that underlies biased agonism, allostery, and basal activity (Hilger et al., 2018; Weis and Kobilka, 2018). Although cryo-EM and X-ray crystallography have populated the Protein Data Bank with hundreds of GPCR-effector complexes, these static endpoints alone cannot reveal the dynamic energetic and/or kinetic bottlenecks that govern signal transduction (Latorraca et al., 2017).

Molecular Dynamics (MD) simulation has become an essential complement, bridging structural data and functional behavior at atomistic resolution (Johnston and Filizola, 2011; Miao and McCammon, 2016). Yet GPCR activation typically occurs over milliseconds and longer timescales, beyond the reach of brute-force MD on commodity hardware. Enhanced-sampling methods, including Metadynamics (Laio and Parrinello, 2002; Barducci et al., 2008), Gaussian-accelerated MD (GaMD) (Miao et al., 2015), supervised MD (Cuzzolin et al., 2016), string-of-swarms approaches (Pan et al., 2008), and the recently introduced On-the-fly Probability Enhanced Sampling (OPES) (Invernizzi and Parrinello, 2020; Invernizzi et al., 2020) have helped narrow this timescale gap by facilitating barrier crossing and improving the characterization of GPCR free-energy landscapes. However, challenges remain in achieving convergence, deriving quantitatively reliable kinetics, and validating computational predictions experimentally. OPES reformulates adaptive biasing as the on-the-fly estimation of a target probability distribution rather than the accumulation of a history-dependent bias, affording better control over convergence and a more robust treatment of suboptimal or multidimensional sets of CVs (Invernizzi and Parrinello, 2020; 2022). Its multi-replica OneOPES extension (Rizzi et al., 2023) combines a convergence-focused replica with exploratory replicas biased along auxiliary CVs and has been applied extensively to GPCR activation free energies.

Beyond receptor activation and deactivation, enhanced sampling simulations have also significantly advanced our understanding of GPCR allostery (Vuckovic et al., 2023), ligand recognition (Han et al., 2025b), and protein–protein interactions (Wang and Miao, 2019; Wang et al., 2025). Recent studies have revealed how allosteric modulators reshape GPCR’s conformational ensembles and communication networks. Enhanced sampling approaches have further improved characterization of ligand binding and dissociation pathways in diverse GPCR systems. In addition, simulations of peptide and G-protein binding or dissociation have provided atomistic insights into the dynamic mechanisms governing GPCR kinetics and signaling.

Here, we will survey recent advances in enhanced-sampling methods for GPCR research, with emphasis on transferable and reproducible simulation protocols. We also highlight how these approaches have revealed allosteric coupling, hidden intermediate states along activation and deactivation pathways, and ligand-dependent modulation of conserved GPCR microswitches (Table 1).

TABLE 1. Comparison of enhanced sampling approaches discussed in this perspective.

Several reviews have already surveyed GPCR conformational dynamics and the simulation methods used to probe them, including general accounts of GPCR dynamics (Latorraca et al., 2017), the development of enhanced conformational sampling for the GPCR activation landscape (Abrol et al., 2022), computational studies of GPCR-G protein interactions (Wang and Miao, 2019) and, more recently, how functional dynamics can be exploited for drug discovery (Conflitti et al., 2025). Our focus is complementary to these earlier efforts. Rather than broadly cataloguing enhanced-sampling algorithms or GPCR conformational dynamics, we concentrate largely on studies from the past 5 years in which enhanced sampling has yielded converged, experimentally testable free-energy landscapes and kinetic rates. We organize these studies along connected signaling trajectory spanning receptor activation, deactivation, allosteric modulation, ligand and peptide binding, and ultimately G-protein dissociation, a process that has received comparatively limited attention in previous reviews.

2 GPCR activation and deactivation

The conformational cycle of GPCRs, encompassing inactive, intermediate, and fully active states coupled to G-proteins, β-arrestins or other effectors, is governed by free energy landscapes shaped by ligand chemistry, lipid environment, conserved microswitches, and structured allosteric water network (Hilger et al., 2018; Weis and Kobilka, 2018). Quantifying these landscapes computationally has long been an aspiration; the studies below show that it is now within reach.

For CV-based enhanced sampling methods, a primary obstacle has been the laborious, error-prone construction of path collective variables (PCVs). In the original Root Mean Square Deviation (RMSD)-based formulation (Branduardi et al., 2007), hereafter referred to as reference PATH (RPATH), the path is defined by the RMSD of the instantaneous structure from the ordered sequences of the reference conformations, so that every intermediate milestone must be tied to a curated real structure and refined recursively. This is the construction used to describe β1-adrenergic receptor (ADRB1) activation in an earlier study by the same group (Aureli et al., 2025a). This limitation was recently overcome by the Euclidean PATH (EPATH) formulation (Aureli et al., 2026). In EPATH, milestones are placed within a two-dimensional Root Mean Square of Deviation (RMSD) space spanned by the inactive and active endpoints, eliminating the need to curate intermediate reference conformations. When combined with OneOPES framework (Rizzi et al., 2023), EPATH reproduces the β1-adrenergic receptor (ADRB1) activation free energy previously obtained with the tailored RPATH CV approach (Aureli et al., 2025a). On the µ-opioid receptor (MOR), both the RPATH and EPATH agree on a ΔG of activation of 18.0 ± 1.0 kcal/mol, but EPATH reached this value in roughly 200 ns versus the ∼800 ns required by RPATH. This is a fourfold reduction in convergence time with no manual milestone curation. Importantly, EPATH naturally widens the conformational ensemble in the transition region, as quantified by per-milestone Shannon Entropy, making the protocol genuinely transferable across pharmacologically diverse class A GPCRs (Aureli et al., 2026).

In a companion study, OneOPES was deployed with biologically motivated CVs (distances within the PIF, DRY, and NPxxY microswitches, together with a YY-cavity hydration descriptor) to map multidimensional free-energy landscapes of apo- and adrenaline-bound ADRB1 in a 1-pamitoyl-2-oleoyl-sn-3-glycero-phosphocholine (POPC)/cholesterol bilayer (Aureli et al., 2025a). Three independent multi-replica simulations converged on ΔG = 13 ± 1.0 kcal/mol for the apo receptor and −3 ± 1.0 kcal/mol for the holo form, with adrenaline preferentially stabilizing and intermediate “pre-active” basin in striking agreement with high-resolution backbone NMR (Grahl et al., 2020). The same framework showed that D^2.50 protonation lowers the activation barrier to 3 ± 1.0 kcal/mol, whereas Na⁺ binding raises it to 15 ± 1.0 kcal/mol, quantitatively rationalizing the inverse-agonist effect of sodium and the activation-coupled protonation switch.

The same problem was approached from a different angle by extending supervised MD into a multiple-walker (mwSuMD) variant that requires no energy bias (Deganutti et al., 2025). With 100 ps walker windows and the single metric (SM) score acceptance metric, mwSuMD reconstructed the vasopressin V2AR:AVP complex to a 3.85 Å RMSD and dissociated AVP about 10 times faster than parent SuMD. Crucially, mwSuMD captured the full sequence from inactive GLP-1R, through PF06882961 binding and TM6 outward rearrangement, to G_s engagement and GDP release. Although no explicit biasing potential was applied, progression was guided by supervised selection of short MD segments and thus, the resulting trajectory characterizes a plausible GPCR-G-protein activation and nucleotide-release cycle rather than unbiased ensemble statistics or physical kinetics.

A more granular dissection of intracellular hydration came from (Aureli et al., 2025b) which mapped the YY-motif (Y^5.58/Y^7.53) water network in ADRB1. In the active basin of the apo receptor the cavity hosts 1 to 2 waters bridging the two tyrosines; adrenalin binding loosens the rigid N^7.45/N^7.49/Y^7.53 scaffold and favors longer water chains, most commonly three to four molecules between the two tyrosines, recasting the intracellular hydration landscape. A collective variable specifically constructed to monitor sodium ion (Na⁺) further revealed that Na⁺ entry is kinetically gated by a hydrophobic cluster (M^2.53/V^3.36/W^6.48) and confirmed the opposing roles of sodium and D^2.50 protonation in shifting the active/inactive equilibrium (Aureli et al., 2025b).

Bridging simulation and structural biology, Bi et al. trapped a near fully-activated intermediate of the A2A adenosine receptor (denoted S4) in complex with mini-Gαsβγ at 2.6 Å resolution (Bi et al., 2025), using an R291A mutation to block the S4 to fully active conformation (S5). ¹⁹F-quantitative Nuclear Magnetic Resonance uncovered the complex conformational states (S4), and GTP-hydrolysis and BODIPY-nucleotide assays showed that S4 initiates rate-limited nucleotide exchange at a non-canonical site 2 while the Gα helical domain remains only partially open. GaMD simulations along with a newly resolved cryo-EM structure validated the S4 conformation and revealed that the S4 to S5 transition requires a 2 Å Gα α5 “lift-up”, an 8° clockwise receptor TM6 rotation, and reorganization of the PIF, CWxP, and NPxxY microswitches, providing the missing dynamic bridge between pre-coupled and fully active complexes.

Recently, metadynamics was used to characterize the interconversion of β-arrestin 1 between distinct orientations within the ghrelin-β-arrestin 1 complex (Gomes et al., 2025). Because 2-µs all-atom MD simulations did not capture transitions between the β1AR-like and neurotensin-like (NTSR1) arrangements, metadynamic was applied to overcome the free-energy barriers separating these states. The enhanced-sampling simulations promoted transitions in both directions and revealed two major free-energy minima corresponding to the alternative β-arrestin orientations. They further showed that PI(4,5)P₂ modulates the relative stability and exchange of these states by anchoring the β-arrestin C-edge loop to the membrane. Thus, enhanced sampling provided mechanistic and energetic description of lipid-dependent β-arrestin orientation that was inaccessible within the timescale of the conventional MD simulations.

Beyond class A, GaMD was applied to dissect activation and deactivation of the adhesion GPCR ADGRG2 in complex with the IP15 Stachel-mimicking agonist and the 4PH-E and 4PH-D antagonists (Adediwura and Miao, 2023). Their simulations sampled active, intermediate and inactive conformations with TM3-TM6 distances of approximately 18.7, 14.7, and 11.1–11.8 Å, respectively. Antagonist binding induced 180°–250° rotations of the W^6.53 toggle switch and triggered an Extracellular loop 2 (ECL2) CWI/SYG residue interplay that selectively stabilized the inactive state, establishing that adhesion GPCR deactivation follows mechanistically distinct, peptide-tuneable pathways.

To make activation and deactivation of GPCR simulations broadly applicable, a well-tempered multiple-walker (32 walker) metadynamics protocol (Calderón et al., 2023) was proposed using the universal A¹⁰⁰ activation index (Ibrahim et al., 2019) (a linear combination of five inter-helix distances) as a single CV. Applied to 5-HT₁A in apo, serotonin-bound, aripiprazole-bound, and G-protein bound, the protocol recovered activation/deactivation profiles compatible with cryo-EM and pharmacology experiments in 3–7 µs of total sampling, and supported a conformational-selection mechanism in which ligands redistribute populations among pre-existing microswitch sub-states.

In another study, enhanced-sampling simulations were combined with free energy calculations and machine learning to uncover how ligands with different efficacy profiles modulate activation of the β2AR (Fleetwood et al., 2021). By simulating β2AR bound to full agonists, partial agonists, antagonists, and inverse agonists, the authors observed that ligands reshape the receptor’s conformational ensembles rather than switching between inactive and active states, supporting a multi-state model for GPCR activation. Data-driven dimensionality reduction and classification approaches revealed that each ligand stabilizes distinct receptor conformations, with signaling hotspots localized around TM5 and the conserved NPxxY motif in TM7 (Fleetwood et al., 2021). Overall, this study provided an important computational framework focused on ligand-induced alteration of the activation free-energy landscapes.

Moreover, an optimized string-of-swarms method was used to expose how the agonist BI-167107 reshapes the β2AR landscape (Fleetwood et al., 2020). The TM5 bulge is strongly coupled to the Ile^3.40/Phe^6.44 connector, whose state probabilistically gates desolvation of the D^2.50 cavity, the NPxxY twist, and the Y^5.58/Y^7.53 interaction, ultimately controlling TM6 outward motion. The data-driven CVs separated active from inactive structures of ten further Class A receptors, demonstrating direct cross-family transferability of the framework.

3 GPCR allostery

For decades, the primary “orthosteric” site has been targeted for drug design of GPCRs. However, the orthosteric site is highly conserved among different GPCRs subtypes and often activates/deactivates multiple GPCRs simultaneously with poor selectivity, thereby causing toxic side effects (Allen and Roth, 2011). Alternatively, allosteric modulators including positive and negative allosteric modulators (PAMs and NAMs) have been discovered to bind to distant “allosteric” sites in GPCRs (Keov et al., 2011). Because the allosteric site is more divergent in residue sequences, allosteric modulators offer higher receptor selectivity and therefore, they serve as promising selective therapeutics for GPCRs (Vuckovic et al., 2023). However, despite decades of research, molecular-level understanding of the general principles that govern the pharmacological effects exerted by GPCR allosteric modulators remains limited.

Burger et al. investigated GPCR allostery in the muscarinic acetylcholine receptors (mAChRs) and determined the first active-state cryo-EM structure of the M5 mAChR and identified a previously unknown extrahelical allosteric binding site for selective PAMs (Burger et al., 2025). The authors showed that the earlier M5-selective PAM ML380 does not bind to previously known extracellular allosteric site, leading to the discovery of a new intracellular extrahelical pocket located between transmembrane helices TM3 and TM4 above intracellular loop 2 (ICL2). GaMD simulations of the ACh–VU6007678 complex confirmed a stable extrahelical TM3–TM4/ICL2 binding mode, revealing interactions with F130^3.52, R146^4.41, and M150^4.45 that stabilize ICL2 and the active M5 receptor conformation, consistent with the cryo-EM structure and mutagenesis data (Burger et al., 2025).

A recent computational study (Kong et al., 2025) further demonstrated the utility of enhanced sampling simulations for elucidating allosteric cooperativity in drug discovery. Using microsecond GaMD simulations combined with MM/GBSA free energy calculations, principal component analysis, dynamic cross-correlation analysis, and elastic network model-based communication pathway analysis, the authors investigated the synergistic effects of orthosteric ligand-binding pocket (LBP) antagonists and allosteric binding function-3 (BF3) site antagonists on the androgen receptor (AR). The simulations revealed that simultaneous binding of orthosteric and allosteric antagonists induced greater conformational rearrangements than either ligand alone, including enhanced flexibility and anticorrelated motions of helix H12 relative to helices H3 and H4. Furthermore, residue interaction networks and signaling pathway analyses identified both short- and long-range allosteric communication pathways linking the BF3 and LBP binding sites to the activation function-2 (AF2) region, providing an explanation for experimentally observed positive cooperativity (Kong et al., 2025). This work highlights how GaMD simulations can capture dynamic allosteric coupling between spatially distinct ligand-binding sites.

A combination of conventional molecular dynamics (cMD) simulations with enhanced sampling metadynamics, protein energy network (PEN) analysis, and transient pocket detection were used to investigate the activation mechanism of the adenosine A1 receptor (A1R) (Maria-Solano and Choi, 2023). Free-energy landscape of receptor activation identified inactive, intermediate, pre-active, and fully active conformational states associated with the TM6 inward-to-outward transition. PEN analysis was then applied to construct residue interaction networks by calculating pairwise residue interaction energies from the MD trajectories. In parallel, geometry-based transient pocket detection (MDpocket) was used to identify cavities that were formed dynamically during the receptor activation. By integrating these analyses, the authors demonstrated that allosteric communication networks evolved during receptor activation, with stronger coupling between extracellular and intracellular regions in the active state and showed that transient allosteric pockets were dynamically connected to these communication pathways (Maria-Solano and Choi, 2023).

A new GaMD, Deep Learning, and free-energy prOfiling Workflow (GLOW) was developed to investigate dynamic mechanism of GPCR allostery (Do et al., 2022). The authors ran a total of 66μs GaMD simulations on 44 GPCR systems with and without PAMs/NAMs and observed that allosteric modulation primarily acts through “conformational selection” in receptor dynamics. Deep learning analysis of residue contact maps and free-energy landscapes further revealed that modulator binding significantly reduces GPCR conformational flexibility and restricts the receptors to fewer low-energy conformational states (Do et al., 2022).

Coarse-grained MD (CG-MD) simulations provided important insights into lipid-mediated allosteric regulation of GPCRs (Song et al., 2019). studied the adenosine A2A receptor in the inactive, active, and active–mini-G_s states embedded in an in vivo-mimetic membrane containing multiple lipid species. Simulations identified nine specific interaction sites for different membrane lipids such as GM3 (monosialotetrahexosylganglioside), cholesterol, and PIP2 (phosphatidylinositol 4,5-bisphosphate), with lipid binding strengths and interaction lifetimes showing clear dependence on the receptor conformational state. In particular, PIP2 preferentially stabilized the active state of receptor by binding at the intracellular TM6/TM7 site, thereby supporting the outward movement of TM6, while simultaneously enhancing A2A receptor–mini-G_s association through lipid-mediated interactions (Song et al., 2019). These findings demonstrate how CG-MD could capture membrane-mediated allosteric mechanisms and revealed how specific lipid–GPCR interactions regulate receptor activation and G-protein coupling.

CG-MD simulations were further combined with native mass spectrometry, mutagenesis, and functional assays to investigate phosphatidylinositol 4,5-bisphosphate (PIP2) interactions with several class A GPCRs, including β1AR, A2AR, and NTSR1 (Yen et al., 2018). The simulations identified conserved PIP2-binding sites at the intracellular ends of the transmembrane helices, where interactions with receptor residues facilitated lipid association. Importantly, experimental and simulation-based analyses showed that PIP2 preferentially stabilized the active receptor states and enhanced the G-protein coupling. These findings established PIP2 as a membrane-associated allosteric modulator capable of influencing GPCR conformational states and G-protein coupling selectivity, highlighting the importance of specific lipid–receptor interactions in GPCR allostery (Yen et al., 2018).

4 Ligand binding and unbinding in GPCRs

GPCRs remain the primary conduit by which cells detect environmental stimuli (Wang and Miao, 2019). Extracellular ligands (endogenous agonists or drug molecules) initiate signaling by binding to GPCRs, triggering conformational changes that enable interaction with heterotrimeric G-proteins and the regulation of downstream signaling. Computational modeling and simulations have become a powerful tool for investigating the mechanisms and kinetics of ligand binding and unbinding to GPCRs at an atomic resolution (Huang and Caflisch, 2011; Pan et al., 2013).

Weighted ensemble (WE) simulations combined with machine learning–based dimensionality reduction and clustering approaches (PCA, t-SNE, and DBSCAN) were used to characterize the ligand dissociation pathways and kinetics of norepinephrine (NE) and epinephrine (EP) from β1AR/β2AR receptors (Han et al., 2025b). The authors identified receptor-specific unbinding mechanisms, showing that NE and EP dissociate from β1AR primarily through the canonical front-side ECL2 pathway, whereas both ligands preferentially exit β2AR through a previously unreported back-side ECL2/TM3–TM4 pathway. Importantly, the simulations could recover the dissociation kinetic rates, with the predicted β1AR–NE k _*off* matching experimental values (Han et al., 2025b).

In another study, multi-microsecond cMD, GaMD (Miao et al., 2015) and WE simulations (Ahn et al., 2021) were combined along with molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) free-energy calculations, and machine learning–based conformational clustering to elucidate the atomistic mechanisms underlying NE binding, dissociation, and G-protein coupling in β1AR/β2AR (Han et al., 2025a). The authors identified a new secondary binding site for NE in β2AR and WE simulations further revealed distinct subtype-specific dissociation pathways for NE, where NE prominently interacts with ECL3 in β2AR, whereas in β1AR, NE dissociation proceeds mainly through TM6/TM7 pathway (Han et al., 2025a).

Adiabatic-bias MD (ABMD) simulations were performed using path CVs to understand the unbinding mechanism of tiotropium from the muscarinic M3 receptor (M3R) and characterize the effects of orthosteric-site mutations on ligand residence time (Coricello et al., 2025). Rather than simulating the complete dissociation process, the simulations focused on the rate-limiting translocation of tiotropium from the orthosteric site to the extracellular vestibule, capturing metastable states along the unbinding pathway. The calculated translocation times showed a good correlation with experimental dissociation rates, highlighting PCV-ABMD as an efficient approach for investigating mutation-dependent kinetic selectivity in GPCRs (Coricello et al., 2025).

The cryo-EM structures of the endogenous agonist epinephrine and the selective agonist dexmedetomidine bound to the G_i-coupled α_2A-adrenergic receptor (α_2A-AR), were resolved revealing important mechanistic differences in ligand recognition and G-protein coupling between α- and β-adrenergic receptors (Lou et al., 2024). Using cryo-EM, GaMD simulations, mutagenesis, and functional signaling assays, the authors demonstrated that epinephrine adopts a distinct binding conformation in α_2A-AR compared with β-adrenergic receptors, due to subtype-specific interactions involving the conserved TM7 residue F427^7.38 in α_2A-AR versus N^7.38 in β-ARs, which act as key determinants of ligand selectivity. Importantly, comparison of α_2A-AR–Gi and β₁-AR–Gi complexes uncovered distinct receptor-specific G_i-coupling interfaces, providing new insights into GPCR-G-protein selectivity and signaling bias (Lou et al., 2024). Overall, this work provides an important structural framework for understanding adrenergic receptor’s subtype specificity, G_i-coupling, and agonist recognition.

Steered MD (SMD) was applied to estimate ligand dissociation rates and residence times in the adenosine A_2A receptor (Akhunzada et al., 2022). Using multiple SMD replicas at different pulling velocities and the Bell–Evans model to relate the force required to induce ligand dissociation at different pulling rates to the underlying dissociation kinetics, the authors predicted the dissociation rate constants and residence times of the antagonist ZMA241385 and agonist NECA. Although the predicted absolute residence times were shorter than the experimental values, the approach successfully showed that ZMA241385 exhibited longer residence time than NECA. The study also highlighted the potential of SMD to integrate ligand unbinding kinetics with binding thermodynamics (Akhunzada et al., 2022).

GaMD enhanced sampling simulations were used to explore the mechanisms of a nonselective antagonist (Caffeine) binding to adenosine A2A receptor (Do et al., 2021). The simulations captured spontaneous ligand binding and unbinding events and importantly, this study identified previously unknown low-energy intermediate binding state located at the extracellular mouth between ECL2 and TM1 (Do et al., 2021). The study further demonstrated that ECL2, which exhibits high sequence divergence among adenosine receptor subtypes, plays a central role in guiding ligand recognition, highlighting the ECL2 region as a promising target site for drug design.

Tau-Random Acceleration MD (τRAMD) was applied to investigate ligand dissociation kinetics and pathways in the β2AR and muscarinic M2 receptor (Kokh and Wade, 2021). By applying an adaptively oriented external force to accelerate ligand unbinding, τRAMD enabled estimation of relative residence times while simultaneously identifying metastable intermediates associated with ligand dissociation. Importantly, the calculated relative residence times showed strong correlation with experimental measurements across different GPCRs, while simulations of the M2 receptor further demonstrated how the positive allosteric modulator LY2119620 altered the iperoxo agonist dissociation pathways and substantially increased its residence time. These findings highlight τRAMD as an efficient approach to linking GPCR ligand-binding kinetics to the receptor structural and dynamic features governing ligand dissociation (Kokh and Wade, 2021).

Metadynamics has been combined with machine learning approaches to investigate GPCR ligand dissociation kinetics (Lamim Ribeiro et al., 2020). The authors integrated Automatic Mutual Information Noise Omission (AMINO) and Reweighted Autoencoded Variational Bayes for Enhanced Sampling (RAVE) with infrequent metadynamics to learn optimized reaction coordinates for the dissociation of morphine and buprenorphine from the μ-opioid receptor. The approach enabled estimation of ligand residence times and dissociation rate constants from microsecond-timescale simulations, with the predicted rates agreeing well with the experimental values. The simulations further revealed that morphine dissociation proceeded through a vestibular intermediate followed by a rate-limiting transition step to the solvent-exposed state, whereas buprenorphine exhibited a distinct kinetic bottleneck associated with the transition from the orthosteric site to the vestibular region (Lamim Ribeiro et al., 2020).

Supervised metadynamics (suMetaD), an enhanced-sampling approach combining supervised MD with metadynamics, was developed to investigate the binding and unbinding kinetics of six antagonists targeting the adenosine A2A receptor (Deganutti et al., 2017). By directing ligand motion along pathways compatible with association or dissociation while simultaneously estimating the energetic profile, suMetaD identified representative transition-state conformations and kinetic bottlenecks along the ligand binding/unbinding pathways. Comparison with surface plasmon resonance (SPR) measurements revealed that the transition-state energetics were strongly influenced by protein–ligand and protein–water interactions. The study further demonstrated how enhanced sampling simulations can complement experimental kinetic measurements by providing structural interpretations of transition states (Deganutti et al., 2017).

5 Peptide binding and unbinding to GPCRs

Peptides can bind to GPCRs and initiate various biological processes, including signal transduction, cellular regulation, immune response, and so on (Sriram and Insel, 2018). Peptides exhibit a unique capability to bind proteins with shallow pockets, which are often difficult to modulate by small molecules. Peptides have therefore emerged as promising drug candidates (Ahrens et al., 2012; Fosgerau and Hoffmann, 2015) and accurate characterization of peptide–protein interactions becomes important for both biological research and drug development.

Peptide GaMD (Pep-GaMD) (Wang and Miao, 2020) simulations were performed on the adhesion GPCR ADGRD1 in both the presence and absence of the G_s protein, with unbound and bound tethered agonist (TA) conformations (Joshi and Miao, 2026). Free energy landscapes were reconstructed through energetic reweighting (Miao et al., 2014), and Density Peak clustering (Chen et al., 2020) was used to identify representative low-energy conformational states along the binding pathway. The simulations captured spontaneous binding of the TA from an extended β-strand conformation into the receptor’s orthosteric pocket (Figure 1A), where it underwent a large conformational transition into the U-shaped α-helical active conformation (Figure 1B) as observed in the cryo-EM structures (Qu et al., 2022; Xiao et al., 2022). Importantly, the study identified two previously unknown intermediate states during binding: an extracellular “I1″ state in which the TA transiently interacted with ECL3, TM5, TM6, and TM7 at the receptor surface, and an “I2″ state in which the TA partially entered the orthosteric pocket before adopting its final active conformation. In the absence of the G_s-protein, Pep-GaMD captured only partial TA binding together with receptor deactivation, demonstrating that G-protein coupling is critical for stabilizing the fully active ADGRD1 conformation and enabling complete tethered agonist binding (Joshi and Miao, 2026).

All-atom Protein-Protein Interaction GaMD (PPI-GaMD) simulations (Wang and Miao, 2022) were applied on the peptide agonist–bound adhesion GPCR complexes (ADGRG2-IP15 and ADGRG1-P7) after removal of the G-protein (Joshi and Miao, 2025). The simulations captured spontaneous dissociation of the peptide agonists from the orthosteric pocket through a conserved extracellular pathway formed by ECL2, ECL3, TM7, and TM1 in both receptors. In addition, the simulations identified distinct low-energy intermediate states for peptides and the receptor in the active, intermediate, and fully inactive conformational states during deactivation. In ADGRG2, two partially dissociated states (PD1 and PD2) of the IP15 peptide agonist were characterized, whereas only a single partially dissociated state was observed for the P7 peptide agonist in ADGRG1, highlighting subtype-specific differences in peptide dissociation. The simulations further revealed structural hallmarks of receptor deactivation, including relaxation of the TM6 kink, inward movement of the intracellular end of TM6, and closure of the orthosteric binding pocket following peptide dissociation (Joshi and Miao, 2025). Overall, these studies provide a detailed computational framework for understanding TA binding and unbinding in adhesion GPCRs, facilitating rational design of peptide modulators targeting ADGRs.

6 Dissociation kinetics of G-proteins from GPCRs

GPCRs couples with different subtypes of intracellular heterotrimeric G-proteins, which subsequently dissociate from their receptors to regulate various downstream signalling pathways (Lefkowitz, 2013; Sungkaworn et al., 2017; Alegre et al., 2021). Structural biology techniques such as X-ray crystallography and cryo-EM provides important insights into structural changes of GPCRs bound to different G-proteins (Piper et al., 2022). However, these techniques primarily offer static snapshots of protein complexes, leaving the dynamic aspects of protein-protein interactions largely unexplored. This further highlights the need to study GPCR–G protein interactions.

The first detailed characterization of G-protein dissociation kinetics from GPCRs was provided by combining PPI-GaMD simulations with fluorescence-and BRET-based kinetic assays (Wang et al., 2025). Using five experimentally characterized GPCR–G-protein complexes, the authors demonstrated that PPI-GaMD simulations (Wang and Miao, 2022) could successfully capture spontaneous G-protein dissociation events and predict dissociation kinetic rate constants that agreed well with the experimental data. The combined free energy profile obtained from PPI-GaMD production simulations, identified three low-energy states, i.e., the dissociated “D,” intermediate “I1”and “I2” states (Figure 2A). The simulations identified two dissociation pathways involving distinct intermediate states: one dominated by interactions between the G-protein α5 helix and receptor ICL2 (Figure 2B), and another mediated primarily through TM5/TM6 contacts (Figure 2E) (Wang et al., 2025). In the intermediate I1 state, receptor residues Y141^ICL2 and Q142^ICL2 formed hydrophobic interactions with the C-terminus of the Gα5 helix (Figure 2D) while in the intermediate I2 state, most of the hydrophobic interactions were formed between receptor TM5/6 and the C terminus of the Gα5 helix (Figure 2G). The cryo-EM structure is shown as reference in Figures 2C,F.

The authors also extended their work to study the effects of allosteric modulators on G_i-protein dissociation from the adenosine A₁ receptor and observed that PAMs of the A₁R, including T62 and MIPS521, significantly reduced G_i-protein dissociation rates and stabilized receptor–G-protein interactions (Wang et al., 2025). Overall, the study provides an important framework for quantitatively characterizing GPCR-G protein signaling kinetics with broader implications for advancing allosteric drug discovery for GPCRs.

7 Discussion

Enhanced sampling studies of GPCRs surveyed here demonstrate the fields’ progression from proof-of-concept applications toward quantitative and predictive characterization of receptor conformational dynamics, activation mechanisms, allosteric modulation, free energy landscapes and signaling kinetics. Three principles emerge with growing clarity. First, GPCR activation is intrinsically multidimensional and intermediate-rich, demanding data-driven and/or biologically motivated collective variables rather than any single geometric proxy or CV-free enhanced sampling approaches. Second, structured water networks and conserved microswitches (PIF, DRY, NPxxY, YY, CWxP) act as allosteric relays whose reorganization is rate-limiting and ligand-tunable. Third, enhanced sampling simulations now make it possible to characterize complete GPCR signaling cycles, from ligand binding and receptor activation to G-protein coupling, GDP release, etc. Major understanding from these studies is that GPCR signaling is fundamentally governed by highly dynamic and interconnected allosteric networks. Enhanced sampling simulations combined with cryo-EM, NMR, machine learning, and free-energy profiling consistently revealed that ligands and allosteric modulators reshape receptor conformational ensembles by stabilizing specific low-energy states and altering communication pathways between extracellular ligand-binding regions and intracellular protein signaling interfaces.

The reviewed studies collectively demonstrate that GPCR function cannot be fully understood through static structures alone, but rather through highly dynamic conformational and interaction networks. Enhanced sampling simulations have revealed how allosteric modulators alter communication pathways and stabilize distinct receptor conformations (Vuckovic et al., 2023; Burger et al., 2025), while ligand binding and dissociation studies provide us with receptor subtype-specific pathways, novel intermediate states, and kinetic rates (Han et al., 2025b). Furthermore, studies of peptide binding (Joshi and Miao, 2026) and G-protein dissociation (Wang and Miao, 2022; Wang et al., 2025) provide atomistic insights into transient intermediates and kinetic processes that regulate GPCR activation and deactivation.

The continued integration of enhanced sampling simulations with deep learning and high-resolution structural biology will likely drive the next generation of GPCR research. Deep learning is already reshaping how GPCR simulation data is analyzed. Data-driven classification and dimensionality reduction of simulated ensembles have been used to identify ligand-specific receptor states and to predict downstream efficacy from the conformations that a ligand stabilizes (Fleetwood et al., 2021). Related approaches have recast ligand-specific differences in GPCR trajectories into representations amenable to deep learning object-recognition architectures, resolving ligand-dependent functional mechanisms in the 5-HT2A and D2 receptors (Plante et al., 2019). In the same vein, the GLOW workflow combined GaMD with deep learning of residue contact maps across 44 GPCR systems to show that allosteric modulators confine receptors to fewer low-energy states (Do et al., 2022). Machine learning models trained directly on large sets of GPCR trajectories have likewise been used to build residue-level activity maps and transition pathways, attributing to individual residues’ contribution to the GPCR activation (Mollaei and Barati Farimani, 2023). A complementary and rapidly growing direction is the use of deep learning to construct the CVs themselves, relacing hand-crafted geometric descriptors with variables learned from data through neural-network based discriminant analysis or by learning the slow modes of the dynamics (Bonati et al., 2020; Bonati et al., 2021). Path CVs can be learned from a deep locally non-linear embedding using a reactive trajectory without hand-picked landmarks (Fröhlking et al., 2024), an approach directly complementary to the RPATH and EPATH strategies discussed above.

Despite substantial progress, enhanced-sampling studies of GPCRs remain limited by several technical challenges. First, CV-based methods depend on the quality of the chosen CVs, and poor CVs can potentially distort free-energy barriers and intermediate states. Even optimized path CVs converge slowly; for example, the MOR RPATH CV required iterative refinement and converged about fourfold more slowly than its Euclidean counterpart (Aureli et al., 2026). CV-free approaches avoid this dependence but could introduce larger energetic noise for reweighting, limiting the accuracy of high-barrier energy landscapes (Miao et al., 2014). Convergence remains difficult to assess, as uncertainties of ∼1 kcal/mol are often estimated from only a few independent replicas, while no standardized criterion exists for determining whether a GPCR free-energy landscape has truly converged. Consequently, independent replicate simulations remain more reliable than block analysis of a single trajectory for assessing convergence (Aureli et al., 2025a). In addition, modelling choices can substantially alter the resulting landscapes. For ADRB1, protonation of D^2.50 reduced the activation barrier to 3 ± 1.0 kcal/mol, whereas Na + binding increased it to 15 kcal/mol, while lipid composition (especially cholesterol) and structured waters could exert important effects (Aureli et al., 2025a; Aureli et al., 2025b). Kinetic predictions remain more challenging than thermodynamic estimates because rate calculations depend on transition-state assumptions and the influence of simulation bias. Although predicted dissociation rates have agreed with experiments in certain systems, validation remains limited across GPCR families (Han et al., 2025b; Wang et al., 2025). Computational cost further restricts sampling, particularly for large GPCR-G protein complexes and multi-replica simulations, leaving class B and C GPCRs, adhesion receptors, oligomers and arrestin complexes relatively underexplored.

Addressing these limitations will require improved method development, reproducibility and experimental validation. Learned CVs are promising but can be circular because their standardized benchmark receptors and open sharing of complete simulations would enable more rigorous comparisons of competing methods (Calderón et al., 2023; Aureli et al., 2025a; Aureli et al., 2026). Finally, extending enhanced-sampling approaches beyond class A GPCRs to other classes such as adhesion GPCRs (Adediwura and Miao, 2023; Joshi and Miao, 2025; 2026), as well as studies of emerging GPCR functional paradigms such as lipid-dependent arrestin coupling (Gomes et al., 2025) and biased signaling (Tikhonova et al., 2013; Madhu et al., 2024) will be essential for establishing these methods as routine tools for GPCR pharmacology.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Myungjin Lee, University of Arizona, United States

Reviewed by: Wen Ma, University of Vermont, United States

Gian Marco Elisi, University of Urbino Carlo Bo, Italy