G Protein-Coupled Receptors (GPCRs) are important drug targets consisting of seven membrane-spanning helices connected by alternating intracellular and extracellular loops, and known for transducing extracellular signals across the cell membrane. Evidence accumulated over the past decade has suggested a pluridimensional functionality of these receptors induced by ligands with different efficacies. Specifically, ligands as diverse as photons, small molecules, and peptides can stabilize different receptor conformations which, in turn, can trigger the activation of different effectors, such as several subtypes of heterotrimeric G proteins, β-arrestins, and G protein-coupled receptor kinases (1). This differential engagement of transducer subtypes with varying signal magnitudes is a phenomenon known as “functional selectivity” or “biased agonism” (2–7), and is a prerequisite for triggering therapeutic or unwanted biological effects of GPCRs via activation of different downstream signaling pathways. However, it is not the only element at play. Temporal analyses of ligand binding to GPCRs, as well as GPCR-induced G protein activation/coupling and β-arrestin recruitment evoked by ligands with different efficacies, demonstrate the existence of another dimension of functional bias by GPCRs that directly impacts our understanding of GPCR signaling and therapeutic implications, thus suggesting the importance of incorporating quantifications of ligand binding and signaling kinetics in modern drug discovery efforts (8–10).

The mechanistic and kinetic bases of GPCR-mediated functional selectivity are poorly understood, notwithstanding the past decade’s technological advances in GPCR functional and structural biology. Among them are the availability of genetically-encoded biosensors for the optical detection of signals from a variety of transduction molecules in several cell types, tissues, and whole organisms (11), as well as revolutionary methodological developments in X-ray crystallography and cryogenic electron microscopy (cryo-EM). The latter have allowed to solve high-resolution experimental structures for 140 unique GPCRs, 520 unique ligand-GPCR complexes, 95 unique GPCRs in complex with heterotrimeric G proteins, 6 unique GPCRs in complex with arrestins, and 1 unique GPCR in complex with G protein-coupled receptor kinases (data retrieved from GPCRdb (12) and the Protein Data Bank (13) on 11/4/2022). Although these structures have provided important insights into ligand-GPCR and GPCR-transducer interactions, they are heavily engineered static snapshots. Thus, inferences of long-distance conformational changes propagating from the ligand binding site in the receptor to its cellular signaling partners, and eventually translating into specific physiological cell responses, remain highly speculative.

To probe the conformational heterogeneity of GPCRs and GPCR complexes and extend information from high-resolution structural methods, researchers have resorted to biophysical techniques such as nuclear magnetic resonance (14–21), double electron-electron resonance spectroscopy (22), hydrogen/deuterium exchange mass spectroscopy (23, 24), and single-molecule fluorescence resonance energy transfer (25–27). However, current technical challenges prevent these techniques from achieving atomic-level precision for the entire GPCR alone or in complex with their natural cellular signaling partners. Although molecular dynamics (MD) simulations can provide a critical bridge between the atomic-level insight from high-resolution structural methods and molecular motions (28), standard MD algorithms limit dynamic explorations to timescales that are shorter than most biological processes notwithstanding their use of massively-parallel high-performance computing platforms. Several enhanced conformational sampling methods have been put forward to overcome these limitations (29), and our group pioneered the use of one of them, i.e., metadynamics (MetaD) [see Box 1 and (30)], for more efficient studies of the conformational plasticity of liganded or unliganded GPCRs embedded in a lipid mimetic environment. The underlying principles of MetaD are that (a) the process to be investigated can be described by a small number of reaction coordinates (collective variables, CVs) and (b) the sum of destabilizing Gaussian potentials that are added to penalize sampled conformational states for faster simulation convergence is the mirror image of the free-energy profile of said process. Although MetaD’s main problem remains that of identifying the right CVs to efficiently study the dynamic process of interest, MetaD-based strategies are increasingly utilized nowadays to study GPCR ligand binding, conformational dynamics, and kinetics (see Figure 1 for illustrations of selected strategies), thus offering an outlet where one could exploit uniquely characterized metastable states of GPCRs as targets for the discovery of functionally selective ligands.

This minireview provides a general overview of MetaD strategies used alone or in combination with algorithms for statistical analyses of simulation data [e.g., Markov State Models (MSMs) (32, 33), information theory-based methods (34), and transfer entropy approaches (35)] and machine learning/artificial intelligence (AI)-based tools (36) to not only characterize ligand-specific conformational states of GPCRs that are difficult to resolve experimentally but also to derive information that can help expedite the GPCR drug discovery process. This includes structural determinants of allosteric communications from the causality of correlated motions, as well as ligand-specific kinetic elements of activation [e.g., (37, 38)]. The computational strategies and applications discussed herein are not exhaustive but include representative examples that have affirmed the power of MetaD in the study of various GPCR processes, including ligand binding, receptor activation, and related kinetics. The different simulation setups and protocols of the studies reported in this minireview are summarized in Table 1 .

To the best of our knowledge, the very first application of MetaD to GPCRs was published by our group more than a decade ago (31) and aimed at efficiently exploring the molecular mechanisms of GPCR ligand recognition (see Figure 1A for an example of output information). Specifically, we studied the free binding of the antagonist naloxone from the water environment to a homology model of the δ-opioid receptor (DOR) based on the X-ray crystal structure of the β₂-adrenergic receptor (β₂AR) as there were no available crystal structures of opioid receptors at the time (see Table 1 for details of system setup and simulations). Although one would have ideally used a single CV to contain the computational cost, which scales exponentially with the number of CVs (55), we could not achieve simulation convergence for this system with less than two CVs and a half microsecond of well-tempered (WT) MetaD simulations using multiple walkers (MW) (56). These simulations revealed, for the first time, unprecedented details of the free binding of a ligand to a GPCR, including the various intermediate energetic states visited by the ligand and their corresponding ligand-receptor interactions. Moreover, by restricting ligand sampling in the bulk region using an approach that allowed the ligand to only move in a conical region centered at the center of mass (COM) of the binding pocket (57), and applying the appropriate correction to the calculated free energy, we could derive ligand binding affinity estimates that were close to experimental values (58–61).

Subsequent applications by our group focused on predicting optimal ligand binding modes for atypical opioid ligands (42), as well as opioid allosteric modulators (41), which would be targeting regions on the receptor that are not conserved and highly flexible. Ligand binding in these studies was described by CVs accounting for the relative position and orientation of the ligands, as well as the number of interactions they established with the receptor. Although these simulations allowed for a thorough and efficient exploration of the dynamic process of ligand-GPCR binding at either orthosteric or allosteric sites (62, 63), they were computationally quite expensive as they required up to 7 μs of simulation time to reach convergence (42).

To ensure faster simulation convergence, a universal single CV that uses the highly conserved Trp^6.48 [superscript refers to the Ballesteros-Weinstein numbering scheme for GPCRs (64)] at the base of the orthosteric binding pocket and the orientation of the receptor in the membrane was proposed for ligand binding/unbinding to/from class A GPCRs using WT-MetaD and a funnel restraint (65) that limits the conformational sampling of the ligand in the bulk water (39) ( Table 1 ). Funnel MetaD (FM) has been successfully applied to a number of GPCR systems (62), and its use made more easily accessible by a recently reported graphical user interface (GUI)-based protocol termed FM Advanced Protocol (FMAP) (66). Also interesting is a protocol that has recently enabled the successful prediction of preferential binding modes of different GPCR systems using conformational ensembles derived from the clustering of MW-MetaD simulations (40). These strategies, however, are still not amenable for high-throughput given their computational cost. If estimates of protein-ligand binding free energies are not required but the goal is focused on assessing the relative stability of binding mode predictions, a reasonable alternative we and others have used in applications to GPCRs [e.g., (43) for a recent example from our lab] is a combination of induced-fit docking and MetaD (67). Based on these growing examples and the ready accessibility of these strategies via user-friendly graphical interfaces, we expect that MetaD-based approaches will become a standard tool for the prediction of GPCR ligand binding in the future, and will make a real impact on drug discovery.

We also were the first to study the activation pathway of a GPCR, specifically rhodopsin (RHO), using a combination of MetaD with adiabatic biased MD simulations (51). In particular, we used the MetaD variant known as path-sampling MetaD (PS-MetaD) to reconstruct the system’s free-energy landscape along predetermined transition trajectories between receptor inactive and active states as a function of the position along the path (s) and the distance (z) from it ( Figure 1B ). Our results suggested that at least four metastable macrostates containing receptor conformations with a different amplitude of the outward movement of transmembrane (TM) helix 6 are sampled by RHO in the transition from inactive to active conformations and these are connected by at least two different pathways (51). The conformations of two of these macrostates were very close to the available inactive and active experimental structures of RHO, whereas the other two macrostates contained conformations representing intermediate states. Notably, subsequent MetaD simulations we carried out on the β₂AR (46) (see section below for additional details and Table 1 ) confirmed the presence of one different inactive-like and one different active-like macrostates in addition to those whose conformations closely resembled the experimentally known inactive and active structures of the receptor. This finding was interesting since standard MD simulations of the β₂AR at the time had only been able to identify one intermediate state in terms of TM3-TM6 separation (68). Similar results were also obtained for the μ-opioid receptor (MOR) by carrying out adiabatic biased MD simulations and PS-MetaD simulations (38), as well as more expensive adaptive-sampling MD simulations (44). In contrast, only one intermediate state between fully inactive and active conformations was identified for the class B glucagon receptor (GCGR) (50), using parallel tempering WT-MetaD and CVs representing two linear combinations of the alpha-carbon root mean square deviation (RMSD) of TM6 to the inactive conformation of GCGR and to the active, closely related, glucagon-like peptide 1 receptor (69, 70). An interesting observation of this study based on a comparison between the computed conformational free-energy landscape associated with the activation of the receptor-agonist complex and that of the receptor-agonist-G protein complex was that the agonist stabilizes the receptor in a preactivated complex before its full activation is achieved by G protein binding (50).

MetaD was also recently used to propose the activation mechanisms of several prototypic class A (45) and class C (49) receptors. These mechanisms were class-specific and revealed an active role of the G protein in promoting conformational changes of the receptor.

Being able to predict the receptor conformations a drug can stabilize is considered the “holy grail” in the drug discovery field as it might inform the drug’s biological outcome. With this in mind, we developed a computational strategy to study the ligand-induced modulation of the free-energy landscape of GPCRs (46) (see Figure 1C for an example of output information). Specifically, we used WT-MetaD  (56) to identify ligand-specific metastable states of the β₂AR, along pre-determined activation pathways between its high-resolution inactive and active structures using adiabatic biased MD ( Table 1 ). The results confirmed a tendency of ligands to stabilize an inactive or active conformation of the receptor depending on their efficacy (46).

That ligands with different efficacies shift the equilibrium between inactive and active states was recently also reported for the human adenosine A2A receptor using MetaD (48). Specific conformational rearrangements of key structural elements were responsible for this shift, including rotameric changes of the conserved Trp^6.48 residue. This movement can be induced by agonists but not inverse agonists, and appears to correlate with the opening of the G protein binding site via disruption of hydrophobic packing which causes TM6 and TM5 to move away from TM3.

MetaD was also recently used to study the interplay between GPCR ligand binding and the coupling of arrestins or G proteins (47). In addition to estimate the binding free energies of ligands with different efficacies at the β₂AR in the presence or absence of Gs or β-arrestin, a combination of WT-MetaD and FM was used to study the free energy landscape of activation and transducer coupling as a function of the distance between the receptor alpha-carbons of Arg^3.50 and Leu^6.34 and the distance between Arg^3.50 and the alpha-carbons of Gαs Glu392 or β-arrestin Val71. This setup ( Table 1 ) allowed to quantify the ligand/transducer cooperative effects on the activation of the β₂AR, thus providing a simple way to predict the functional bias of a ligand (47).

Knowledge of GPCR ligand binding kinetics, and especially drug-target residence times and involved molecular determinants, is highly desirable because it can serve as an effective guiding principle to select candidate molecules for clinical development (10). However, this collective information is difficult to obtain because on and off rates depend on the height of the highest energy barrier of the transition state between bound and unbound conformations, and this height is difficult to estimate both experimentally and computationally. MD simulations can help study ligand binding at the active site, but binding is a rare event on microscopic timescales, and as such, it is difficult to sample. Although microsecond-scale MD simulations allow to obtain multiple binding events at a GPCR embedded in an explicit lipid-water environment, from which it is possible to derive k_on estimates (71), much longer timescales would be required for the dissociation of the ligand from a GPCR orthosteric site, making it very difficult to derive k_off estimates from unbiased, standard MD simulations. To attempt more efficient predictions of dissociation rates from a GPCR, we recently designed a computational strategy that combines machine learning and infrequent MetaD (37). Specifically, the strategy used the automatic mutual information noise omission (AMINO) algorithm (72) to enable the robust and automated selection of non-redundant molecular features from unbiased MD simulations that can then be used in the reweighted autoencoded variational Bayes for enhanced sampling (RAVE) method (73, 74) to learn optimal CVs. These CVs were used in infrequent MetaD simulations and allowed to efficiently study the unbinding kinetics of two opioid drugs with very different residence times (morphine and buprenorphine) and to predict dissociation rates for these drugs from the MOR that were within one order of magnitude from experimental values (75). The simulations also gave structural information about rate-limiting transition states and metastable poses that can in principle be used in the design of drugs with a desired kinetic profile. One important lesson from these simulations was that accurate rate estimates require consideration of ligand solvation, as also concluded by others (52, 65).

Recently, MetaD was also used to investigate the dissociation of a potent drug targeting the MOR, specifically fentanyl, from the receptor (54). Not only did these simulations provide information about fentanyl’s binding kinetics and mechanism, but they suggested a role for the protonation state of the conserved residue His^6.52 in modulating fentanyl’s affinity and binding pose. While a recent study (76) reporting cryo-EM structures of the MOR bound to fentanyl showed that a His(6.52)Ala mutation reduced fentanyl-induced MOR activation of G protein and β-arrestin signaling, it did not support the predicted binding poses of fentanyl using MetaD or other computational methods (53, 54, 77).

Ligand binding kinetics is however not the only temporal aspect that is important to understand GPCR signaling. To tackle the kinetics of GPCR activation, we designed a computational strategy that uses a combination of PS-MetaD and the maximum caliber principle to not only thoroughly characterize the conformational space sampled by the receptor along its activation path, but also to derive unbiased kinetic rates from simulations with transitions accelerated by a bias potential (38) ( Figure 1D ). We showed that the approach could efficiently sample conformational transitions between the crystal structures of inactive and active GPCRs, specifically the MOR, at a ∼2 orders of magnitude computational cost reduction with respect to a more expensive high-throughput MD adaptive sampling protocol run on distributed computational resources. We also demonstrated that the strategy yielded thermodynamic and kinetic properties of MOR activation, as well as information about ligand-induced allosteric communications across the receptor when combined with n-body information theory (34), that were in agreement with those obtained by more computationally expensive adaptive sampling protocols.

Identifying the molecular determinants that underlie the different dynamic and kinetic behaviors of GPCRs is likely to provide fundamental information for drug discovery, and is therefore the focus of much research. However, obtaining this information is not a trivial undertaking, either experimentally or computationally. While strategies combining different variants of MetaD with state-of-the-art algorithms for statistical analyses of simulation data and/or AI-based tools have made strides in the study of GPCR ligand binding, activation, and related kinetics, the majority of studies published to date are limited to the receptors alone, although recent studies are focusing more and more on receptor complexes and the combined effect of ligand and transducer on the dynamic behavior of the receptor. To better understand the underlying principles of GPCR biased signaling, attention will probably need to be shifted more towards the appropriate characterization of the dynamic behavior of interacting signaling proteins, thus providing important insight into receptor-transducer complex stability and/or receptor-transducer activation kinetics. These studies will most likely require developing more sophisticated MetaD strategies or hybrid approaches and a community effort to disseminate best practices for MetaD simulations, as well as system setups, parameters, protocols, etc. to more expeditiously advance knowledge of GPCR dynamics and their relation to functional selectivity.

MF wrote the first draft of the manuscript. LS-E and BF made the table and figure, respectively. All authors contributed to the article and approved the submitted version.