Introduction

Supervised molecular dynamics^{1, 2} (SuMD) is an efficient technique for studying ligand-receptor binding and unbinding pathways; here, we present the multiple walker enhancement (mwSuMD) to study a significantly wider range of structural transitions relevant to the drug design. We validated the method by applying it to G protein-coupled receptors (GPCRs), since their inherent flexibility is essential to their function and because these are the most abundant family of membrane receptors in eukaryotes³ and the target for more than one-third of drugs approved for human use⁴.

Vertebrate GPCRs are subdivided into five subfamilies (Rhodopsin or class A, Secretin or class B, Glutamate or class C, Adhesion, and Frizzled/Taste2) according to function and sequence^{5, 6}. X-ray and cryo-electron microscopy (cryo-EM) show that GPCRs possess seven transmembrane (TM) helices connected by three extracellular loops (ECLs) and three intracellular loops (ICLs), with an extended and structured N-terminal extracellular domain (ECD) in all subtypes, but class A. The primary function of GPCRs is transducing extracellular chemical signals into the cytosol by binding and activating four G protein families (G_s/olf, G_i/o, G_12/13, and G_q/11) responsible for decreasing (G_i/o) or increasing (G_s/olf) cyclic adenosine-3’,5’- monophosphate (cAMP), and generating inositol-1,4,5-triphosphate (IP₃) and diacylglycerol (DAG) to increase Ca²⁺ intracellular levels (G_q)⁷.

Many aspects of GPCR pharmacology remain elusive: for example, the structural determinants of the selectivity displayed towards specific G proteins or the ability of certain agonists to drive a preferred intracellular signaling pathway over the others (i.e. functional selectivity or bias)⁸. GPCRs are challenging proteins to characterize experimentally due to their inherent flexibility and the transitory nature of the complexes formed with extracellular and intracellular effectors. Importantly, agonists can allosterically modify the receptor selectivity profile by imprinting unique intracellular conformations from the orthosteric binding site. The mechanism behind these phenomena is one of the outstanding questions in the GPCR field⁹.

Molecular dynamics (MD) is a powerful computational methodology that predicts the movement and interactions of (bio)molecules in systems of variable complexity, at atomic detail. However, classic MD sampling is limited to the microsecond or, in the best conditions, the millisecond time scale^{10, 11}. For this reason, different algorithms have been designed to speed up the simulation of rare events such as ligand (un)binding and conformational transitions. Amongst the most popular and effective ones, are metadynamics¹², accelerated MD (aMD)¹³, and Gaussian-accelerated MD (GaMD)¹⁴, which introduce an energy potential to overcome the energy barriers preventing the complete exploration of the free energy surface, thus de facto bias the simulation. Such methods have been used to study GPCR activation^{15, 16}. Energetically unbiased MD protocols, on the other hand, comprise weighted ensemble MD (weMD)¹⁷ and SuMD^{1, 18}, which have largely been applied to (unbinding) small molecules, peptides, and small proteins^{1, 18–22}. Since SuMD is optimized for (un)bindings, we have designed mwSuMD to address more complex conformational transitions and protein-protein associations. GPCRs preferentially couple to very few G proteins out of 23 possible counterparts^{9, 23}. It is increasingly accepted that dynamic and transient interactions determine whether the encounter between a GPCR and a G protein results in productive or unproductive coupling²⁴. MD simulations are considered a useful orthogonal tool for providing working hypotheses and rationalizing existing data on G protein selectivity. However, so far, it has not delivered as expected. Attempts so far have employed energetically biased simulations or have been confined to the Gα subunit^{15, 16, 23}.

We tested mwSuMD on a series of increasingly difficult structural transitions involving both class A and class B1 GPCRs. Firstly, we validated the method on the nonapeptide arginine vasopressin (AVP) by simulating binding (dynamic docking) and unbinding paths from the vasopressin 2 receptor (V₂R). Dynamic docking, although more computationally demanding than standard molecular docking, provides insights into the binding mode of ligands in a fully hydrated and flexible environment. Moreover, it informs about binding paths and the complete mechanism of formation leading to an intermolecular complex, delivering in the context of binding kinetics²⁵ and structure-kinetics relationship (SKR) studies²⁶. We then show the G_s and G_i G proteins binding to the adrenoreceptor β₂ (β₂ AR) and the adenosine 1 receptor (A₁R₎, respectively, starting from different conditions. GPCRs preferentially couple to very few G proteins out of 23 possible counterparts^{9, 23}.

A further case study is the activation of the class B1 GPCR glucagon-like peptide-1 receptor (GLP-1R) by the small molecule PF06882961. GLP-1R is a validated target in type 2 diabetes and probably the best-characterized class B1 GPCR from a structural perspective. GLP-1R is the only class B1 receptor with structurally characterized non-peptidic orthosteric agonists, which makes it a model system for studying the druggability of the entire B1 subfamily. After GLP-1R agonist binding and activation, the coupling of G_s and the release of GDP, which is the rate-limiting step of the G protein activation, was simulated for the first time using an energy-unbiased method.

The last GPCR key process addressed by mwSuMD is the heterodimerization in the membrane between the adenosine receptor A₂ (A_2AR) and the dopamine receptor D₂ (D₂R). The A_2AR:D₂R heterodimer²⁷ is a therapeutic target for neurodegenerative diseases, Parkinson’s disease, and schizophrenia^28–30 due to the reciprocal antagonistic allosteric effect between monomers³¹. A_2AR activation reduces the binding affinity of D₂R agonists, while A_2AR antagonists enhance the dopaminergic tone by decreasing the adenosine negative allosteric modulation on D₂R. Heterobivalent ligands able to inhibit A_2AR and activate D₂R represent a valuable pharmacological tool³² and, in principle, therapeutic options for conditions characterized by a reduction of dopaminergic signaling in the central nervous system. The successive dynamic docking of the heterobivalent ligand compound 26³³ to the heterodimer suggested by mwSuMD produced a ternary complex stabilized by lipids and is indicative of the usefulness of mwSuMD for illuminating the molecular events involved in GPCR function.

Results and Discussion

Short mwSuMD time windows improve the AVP dynamic docking prediction

Arginine vasopressin (AVP) is an endogenous hormone (Figure S1a) that mediates antidiuretic effects on the kidney by signaling through three class A GPCR subtypes: V_1a and V_1b receptors activate phospholipases via G_q/11, while the V₂ receptor (V₂R) activates adenylyl cyclase by interacting with G_s ³⁴ and is a therapeutic target for hyponatremia, hypertension, and incontinence³⁵. AVP is amphipathic and in the bound state interacts with both polar and hydrophobic V₂R residues located on TM helices and ECLs (Figure S1b). Although AVP presents an intramolecular C1-C6 disulfide bond that limits the overall conformational flexibility of the backbone, it has a high number of rotatable bonds, making dynamic docking complicated³⁶. We assessed the performance of mwSuMD and the original version of SuMD in reconstructing the experimental V₂R:AVP complex using different settings, simulating a total of 92 binding events (Table S1). As a reference, the AVP RMSD during a classic (unsupervised) equilibrium MD simulation of the X-ray AVP:V₂R complex was 3.80 ± 0.52 Å (Figure S2). SuMD^{1, 18} produced a minimum root mean square deviation (RMSD) to the cryo-EM complex of 4.28 Å, with most of the replicas (i.e., distribution’s mode) having an RMSD close to 10 Å (Figure S3a). mwSuMD, with the same settings (Figure S3b, Table S1) in terms of time window duration (600 ps), metric supervised (the distance between AVP and V₂R), and acceptance method (slope) produced slightly more precise results (i.e., distribution’s mode RMSD = 7.90 Å) but similar accuracy (minimum RMSD = 4.60). Supervising the AVP RMSD to the experimental complex rather than the distance (Figure S3c) and using the SMscore (Equation 1) as the acceptance method (Figure S3d) worsened the prediction. Supervising distance and RMSD at the same time (Figure S3e), employing the DMscore (Equation 2), recovered accuracy (minimum RMSD = 4.60 Å) but not precision (distribution mode RMSD = 12.40 Å). Interestingly, decreasing the time window duration from 600 ps to 100 ps impaired the SuMD ability to predict the experimental complex (Figure 1a), but enhanced mwSuMD accuracy and precision (Figure 1b-d). The combination of RMSD as the supervised metric and SMscore produced the best results in terms of minimum RMSD and distribution mode RMSD, 3.85 Å and 4.40 Å, respectively (Figure 1d, Video S1), in agreement with the AVP deviations in the equilibrium MD simulation of the X-ray AVP:V₂R complex (Figure S2). These results confirm the inherent complexity of reproducing the AVP:V₂R complex via dynamic docking and suggest that short time windows can improve mwSuMD performance on this system. However, it is necessary to know the final bound state to employ the RMSD as the supervised metric, while the distance is required to dynamically dock ligands with unknown bound conformation as previously reported^{1, 23}. Both distance and RMSD-based simulations delivered insights into the binding path and the residues involved along the recognition route. For example, mwSuMD pinpointed V₂R residues E184^ECL2, P298^ECL3, and E303^ECL3 (Figure S4a) as involved during AVP binding, although not in contact with the ligand in the orthosteric complex. None of them are yet characterized through mutagenesis studies according to the GPCRdb³⁷.

Figure 1.

Further to binding, a SuMD approach was previously employed to reconstruct the unbinding path of ligands from several GPCRs^{2, 38}. We assessed mwSuMD’s capability to simulate AVP unbinding from V₂R. Five mwSuMD and five SuMD replicas were collected using 100 ps time windows (Table S1). Overall, mwSuMD outperformed SuMD in terms of time required to complete a dissociation (Figure S5, Video S2), producing dissociation paths almost 10-fold faster than SuMD. Such rapidity in dissociating inherently produces a limited sampling of metastable states along the pathway, which can be compensated by seeding classic (unsupervised) MD simulations from configurations extracted from the unbinding pathway^{39, 40}. Here, the set of V₂R residues involved during the dissociation was comparable to the binding (Figure S4b), though ECL2 and ECL3 were slightly more involved during the association than the dissociation, in analogy with other class A and B GPCRs^{20, 39}.

Binding simulations between G proteins and a class A GPCR

We tested the ability of mwSuMD to simulate the binding between a prototypical class A receptor, the β₂ adrenoreceptor (β₂ AR), and the stimulatory G protein (G_s), without energy input. mwSuMD simulations started from the intermediate, agonist-bound conformation of the β₂-AR and the inactive G_s to resemble pre-coupling conditions. Three mwSuMD replicas were performed by supervising the distance between G_s helix 5 (α5) and β₂ AR as well as the RMSD of the intracellular end of TM6 to the fully active state of the receptor (Table S1). To monitor the progression of the simulations, we computed the RMSD of the Cα atoms of the Gα and Gβ subunits to the experimental complex⁴¹ (Video S3, Figure 2ab). During two out of three replicas, both Gα and Gβ achieved distance values close to 5 Å (minimum RMSD = 3.94 Å and 3.96 Å respectively), in good agreement with the reference (the β₂ AR:G_s complex, PDB 3SN6, Figure 2c). The flexibility of G_sβ is backed by both MD and cryo-EM data suggesting G protein rocking motions around G_sα:receptor interactions^{20, 42}.

Figure 2.

Usually, ICL3 of the GPCR and the G_s protein loop hgh4, which is not found in other Gα⁴³, are masked out from deposited cryo-EM structures due to their high flexibility and therefore low resolution. During our simulations, these two loops formed polar intermolecular interactions through R239^ICL3, R260^ICL3, K235^ICL3, and E322^hgh4.12, D323^hgh4.13. Further transient interactions, not visible in the experiential structures, involved a mix of conserved and unique residues forming hydrogen bonds (Table S2): R63^ICL1-E392 ^α5.24, K232^5.71-D378^Α5.10, K235^5.74-D378 ^α5.10, K235^ICL3-D343^H4.13, K267^6.29-L394α5, R239^ICL3-E314^hgh4.04, and S137^3.56-D381 ^α5.13. None of the interactions reported in Table S2 is evident from the experimental β₂ AR:G_s complex, implying that mwSuMD can deliver useful working hypotheses for mutagenesis and spectroscopic experiments from out-of-equilibrium simulations. Results also suggest that the G_s binding is driven by a combination of conserved and unique transitory interactions with β₂ AR, possibly contributing to G protein selectivity. We speculate that the conserved interactions would be necessary for the binding, regardless of the receptor:G protein couple involved, while the transitory interactions due to the dynamics of G protein binding should produce an effective engagement of the G protein.

A possible pitfall of the above-reported G_s:β₂ AR mwSuMD binding simulation is that G proteins bear potential palmitoylation and myristoylation sites that anchor the inactive trimer to the plasma membrane^{44, 45}, de facto restraining possible binding paths to the receptor. Recently, the G_i binding to A₁R was simulated by combining the biased methods GaMD with SuMD⁴⁶ but without taking into account the role played by membrane-anchoring post-translational modifications on the G_i binding pathway. To address this point and test the possible system dependency of mwSuMD, we prepared a different class A GPCR, the adenosine A1 receptor (A₁R), and its principal effector, the inhibitory G protein (G_i) considering the G_iα residue C3 and G_γ residue C65 as palmitoylated and geranylgeranylated respectively and hence inserted in the membrane. Both classic (unsupervised) and mwSuMD simulations were performed on this system for comparison (Video S4, Figure 2d). In about 50 ns of mwSuMD, the G_iα subunit engaged its intracellular binding site on A₁R and formed a complex in good agreement with the cryo-EM structure (PDB 6D9H, RMSD ≍ 5 Å). For comparison, 1 μs of cMD did not produce a productive engagement as the G_iα remained at RMSD values > 40 Å (Figure 2d), suggesting the effectiveness of mwSuMD in sampling G protein binding rare events without the input of energy. The membrane anchoring affected the overall G_i binding and the final complex, which was rotated compared to the experimental structure due to the lipidation of G_iα and G_γ (Figure 2e). This suggests that future, more comprehensive studies of G protein binding and activation should consider several G protein orientations around the receptor as the starting points for mwSuMD simulations, to evaluate as many binding paths as possible

PF06882961 binding, GLP-1R activation, and G_s protein activation

The GLP-1R has been captured by cryo-EM in both the inactive apo (ligand-free) and the active (G_s-bound) conformations and in complex with either peptide or non-peptide agonists^47–52. In the inactive apo GLP-1R, residues forming the binding site for the non-peptide agonist PF06882961 are dislocated and scattered due to the structural reorganization of the transmembrane domain (TMD) and extracellular domain (ECD) (Figure S6) that occurs on activation. Moreover, GLP-1R in complex with GLP-1 or different agonists present distinct structural features, even amongst structurally related ligands (Figure S7). This complicates the scenario and suggests divergent recognition mechanisms amongst different agonists. We simulated the binding of PF06882961 using multistep supervision on different metrics of the system (Figure 3) to model the structural hallmark of GLP-1R activation (Video S5, Video S6).

Figure 3.

Several metrics were supervised in a consecutive fashion. Firstly, the distance between PF06882961 and the TMD as well as the RMSD of the ECD to the active state (stage 1); secondly, the RMSD of ECD and ECL1 to the active state (stage 2); thirdly, the RMSD of PF06882961 and ECL3 to the active state (stage 3); lastly, only the RMSD of TM6 (residues I345-F367, Cα atoms) to the active state (stage 4). The combination of these supervisions produced a conformational transition of GLP-1R towards the active state (Figure 3, Video S6). Noteworthy, mwSuMD, like any other CV-based technique, requires some degree of knowledge of the system simulated. The sequence of these supervisions was arbitrary and does not necessarily reflect the right order of the steps involved in GLP-1R activation. This kind of planned multistep approach is feasible when the end-point receptor inactive and active structures are available, and the inherent flexibility of different domains is known. In class B GPCRs, the ECD is the most dynamic sub-structure, followed by the ECL1 and ECL3 which display high plasticity during ligand binding^{20, 53}. For this reason, we first supervised these elements of GLP-1R, leaving the bottleneck of activation, TM6 outward movement, as the last step. However, the protocol employed can be tweaked to study how each conformational transition takes place and influences the receptor domains. Structural elements not directly supervised, such as TM1 or TM7, displayed an RMSD reduction to the active state because they were influenced by the movement of supervised helixes or loops. For example, the supervision of ECL3 (stage 3) and TM6 (stage 4) facilitated the spontaneous rearrangement of the ECD to an active-like conformation after the ECD had previously experienced transient high flexibility during stages 2 and 3 (Figure 3).

It is worth noting that the only inactive, almost full-sequence, GLP-1R structure available, was solved with an antibody bound to the ECD. This crystallization strategy for stabilizing the receptor might have forced the ECD in a closed conformation that engages the EC vestibule of GLP-1R and possibly restrained the whole TMD in an altered conformation that deviates from the active conformation. This might explain why the RMSD of the TMD elements monitored during the simulation rarely reach values lower than 3 or 4 Å. During the supervision of ECL3 and PF06882961 (stage 3), we observed a loosening of the intracellular polar interactions that stabilize GLP-1R TM6 in the inactive state. As a result, the subsequent supervision of TM6 (residues I345-F367, Cα atoms) rapidly produced the outward movement of TM6 towards the active state, in the last step of the mwSuMD simulation (stage 4). Taken together, these results suggest a concerted conformational transition for ECD and ECL1 during the binding of PF06882961 and an allosteric effect between ECL3 and the bottom of TM6. While the intracellular polar interactions were destabilized by the ECL3 transition to an active-like conformation (stages 2 and 3), the outward movement of TM6 (stage 4) did not favor the closure of ECL3 towards PF06882961, which appear to be driven by direct interactions between the ligand and R310^5.40 or R380^7.35. Interestingly, the mwSuMD simulation during Stage 4 (TM6 supervision) sampled a counterclockwise rotation of the helix (Figure 4a), consistent with the GLP-1R cryo-EM structures in the active, G_s-coupled state^{50, 54}. To the best of our knowledge, this is the first time an MD simulation captures the TM6 rotation upon receptor activation as results reported so far are largely limited to the TM6 opening and kinking⁵⁵.

Figure 4.

Encouraged by results previously obtained on G_i binding to β₂ AR and A₁R (Figure 3), we extracted the GLP-1R active conformation described above and simulated the G_s binding to its intracellular side. Starting from the inactive, membrane-anchored G_s, we performed three independent mwSuMD replicas by supervising the distance between G_s helix α5 and GLP-1R residues located at the intracellular binding interface. All three mwSuMD replicas showed the G_s approaching GLP-1R, with two out of three reaching an RMSD of the G_sα subunit close to or less than 10 Å, compared to the experimental complex 7LCI (Figure 4b). Replica 2, in particular, well reproduced the cryo-EM GLP-1R:G_s complex. According to the model of G protein activation, the G protein binding has the effect of allosterically stabilizing the orthosteric agonist in complex with the receptor⁵⁶ and destabilizing the guanosine-diphosphate (GDP) bound to Gα, triggering its release and exchange with guanosine-triphosphate (GTP)⁵⁷, upon opening of the G protein alpha-helical domain (AHD). In accordance with this model, PF06882961 and GDP were respectively stabilized and destabilized during the simulated G_s association (Figure 4c,d). The analysis of atomic contacts along the binding path of G_s to GLP-1R highlights a few persistent interactions not observed in the equilibrium MD simulations of GLP-1R:G_s cryo-EM complexes⁵⁶; in particular, we propose the hidden interaction between D344^6.33 and R385 ^α5 to be pivotal for G_s coupling (Table S3).

We further investigated the G_s activation mechanism by extending Replica 2 (Figure 4b) and supervising the opening of the G_s alpha-helical domain (AHD), which is considered a necessary step to allow GDP release from the Ras-like domain⁵⁸. We first easily obtained the opening of AHD in the extended Replica 2 and successively supervised the GDP unbinding in further three replicas, seeded after the AHD opening. In one of these three mwSuMD simulations, the nucleotide dissociated from G_s (Figure 4e-h). Video S7 shows the full G_s binding, AHD opening, and GDP release. Interestingly, no conformational changes of the β6-α5 loop happened before or during the GDP dissociation, suggesting that its conformational change as captured in the nucleotide-free GLP-1R:G_s complex (Figure S8) occurs after the GDP release as a result of the loss of binding stabilization, rather than being an initiator of the GDP dissociation. Our simulation is consistent with a major role played by α5 tilting and allosteric communication through internal structural elements like the β2-β3 strands, which weaken both GDP phosphate and purine binding sites as previously suggested by other groups^59–61.

The heterodimerization between A_2A and D₂R, and binding simulations of the heterobivalent ligand compound 26

The current structural model of the A_2AR:D₂R heterodimer is that TM4 and TM5 from both receptors contribute to the primary dimer interface, although the involvement of TM7 is not ruled out⁶². Following this interaction model, we first dynamically docked A_2AR and D₂R in an explicit 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane model, then simulated the binding of the heterobivalent compound 26³³ (CP26) to the preformed A_2AR:D₂R heterodimer (Video S8). Since membrane proteins are characterized by slow lateral diffusion⁶³, we encouraged the encounter between A_2AR and D₂R by inputting energy using metadynamics and adiabatic MD during mwSuMD (hybrid mwSuMD/metadynamics/aMD), followed by 1.5 μs of classic MD (cMD) to relax the system and check the stability of the A_2AR:D₂R interactions.

During the first 200 ns of simulation with energy bias (Figure 5a,c and Figure S9a), A_2AR and D₂R rapidly moved close to each other and reached a distance of about 30 Å (computed between centroids), before stabilizing at around 40 Å (Figure 5a). The computed molecular mechanics combined with the Poisson–Boltzmann and surface area continuum solvation (MM-PBSA) binding energy suggested two energy minima (Figure 5c) at these distances. The successive cMD simulation did not produce remarkable changes in the distance between receptors (Figure 5b), although the energy fluctuated before reaching about −10 kcal/mol, at the end of the simulation (Figure 5d). The sharp energy minima after 25 and 150 ns were due to the high number of direct contacts between A_2AR and D₂R (Figure S9), favored by the energy added to the system. When the input of energy bias was stopped (Figure 5b,d) the POPC residues re-equilibrated at the interface between proteins and mediated intracellular polar interactions between R150^{4.40 D2R}, Y146^{4.36 D2} and R199^{5.60 A2A}, Y103^{3.51 A2A} as well as extracellular polar interactions between the top of TM4^D2, TM5^D2 and TM5^A2A, TM6^A2A (Figure 5f), suggesting that the A_2AR:D₂R heterodimerization relies on lipids to mediate short-range interactions between receptors.

Figure 5.

The dynamic docking of the heterobivalent ligand C26 further stabilized the A_2AR:D₂R dimer (Figure 5e), in line with experimental data³³. C26 reached the bound state rapidly, inserting the agonist pharmacophore within the D₂R orthosteric site (Figure S10, Video S8), while the pyrazole-triazole-pyrimidine scaffold remained in a metastable complex with A_2AR, before completely binding the orthosteric site at the end of the simulation (Figure S11, Video S8). In the final state, the long linker between pharmacophores extended over the top of the interface formed by A_2AR and D₂R at the level of the receptors’ ECL2 (Figure 5g). A network of polar interactions between POPC, Y179^A2A, and Y192^D2 contributed to stabilizing this ternary complex. Interestingly, the latter residues were pinpointed as important for A_2AR:D₂R interactions⁶². From a binding energy perspective, C26 reached the most stable configurations between 80 and 100 ns (Figure S12), before the pyrazole-triazole-pyrimidine component of the ligand completed the binding to A_2AR. This suggests some contribution of the linker to the overall stability of the ternary complex with A_2AR and D₂R. Two out of four mwSuMD replicas produced A_2AR:D₂R:C26 ternary complexes with C26 engaged both by the orthosteric site of A_2AR and D₂R, while in the remaining two replicas the A₂AR pharmacophore remained stacked on the extracellular vestibule of the receptor, although in the proximity of the binding site (Figure S10). Considering these results, we suggest that the medicinal chemistry optimization of bivalent ligands should involve the tuning of the on-rate kinetic constant k_on measured on each monomer of the heterodimer, to differentiate which pharmacophoric scaffold is kinetically more affected by the chemical linking and therefore aiding the surgical improvement of the bivalent ligand’s affinity.

Discussion

Classic MD simulations sample phase space with an efficiency that depends on the energy barrier between neighboring minima. Processes like (un)binding and protein activation require the system to overcome numerous energy barriers, some of which create bottlenecks that slows the transition down to the millisecond, or second, time scale. To overcome some of these limits, we have developed and tested on complex structural events characterizing GPCRs an energetically unbiased adaptive sampling algorithm, namely multiple walker mwSuMD, which is based on traditional SuMD, while drawing on parallel multiple replica methods^{64, 65},

mwSuMD performed similarly to SuMD for the dynamic docking of AVP to V2R when time windows of 600 ps were employed. Time windows of 100 ps remarkably improved mwSuMD. Usually, dynamic docking is performed to either predict the geometry of complexes or sample the binding path of an already-known intermolecular complex, or both. The RMSD of AVP to the experimental coordinates as the supervised metric produced the best results. Consequently, the RMSD should be the metric of choice to study the binding path of well-known intermolecular complexes. The distance, on the other hand, is necessary when limited structural information about the binding mode is available. In the absence of structural information regarding the final bound state, it is possible to sample numerous binding events employing mwSuMD and evaluate the final bound states rank by applying end-point free energy binding methods like the molecular mechanics energies combined with the Poisson– Boltzmann or generalized Born and surface area continuum solvation (MM/PBSA and MM/GBSA⁶⁶) models. Our simulations propose that remarkable predictivity can be obtained with distance-driven mwSuMD, as demonstrated by the lowest deviation from the experimental AVP:V2R complex. Remarkably, the dissociation of AVP from V2R was simulated much more rapidly by mwSuMD than by SuMD, suggesting it is an efficient tool for studying the dissociation of ligands from GPCRs.

mwSuMD was further tested on the G_s and G_i binding to β₂ AR and A₁R, respectively. mwSuMD produced G protein:GPCR complexes in remarkable agreement with experimental structural data without the input of energy in a few hundred nanoseconds when starting from inactive G_s and the intermediate active β₂ AR, or a few tens of nanoseconds when considering the active-state A₁R and G_i anchored to the plasma membrane through the palmitoylation and the geranylgeranylation of Gαγ^{44, 45, 67}.

We increased the complexity of binding simulations by considering GLP-1R and the non-peptide agonist PF06882961. Using mwSuMD, we obtained a binding of PF06882961in good agreement with the cryo-EM structure, followed by an active-like conformational transition of GLP-1R. The choice of the metrics supervised was driven by structural data available⁵⁰ and extensive preparatory MD simulations, however, binding routes are possible from either the bulk solvent or the membrane^{39, 68, 69}. Future studies on GLP-1R and other class B1 GPCRs should consider different starting points for the ligand and alternative apo receptor conformations to improve the sampling.

mwSuMD was able to simulate the G_s binding to the active GLP-1R and the subsequent GDP release. To our knowledge, this is the first time the whole sequence of events leading from an inactive GPCR to the GDP release is simulated. Remarkably, this was obtained without biasing the simulations with external energy. Our results suggest that the full rotation and elongation of α5 as in cryo-EM structures could occur after the GDP release, supporting the role of hidden, metastable interactions as the driving force of G protein coupling and selectivity, in accordance with recent work on GLP-1R⁷⁰.

The final case study was the dimerization process between A_2AR and D₂R in a membrane model. To speed up the encounter between receptors, we introduced an energy bias in the form of abMD and metadynamics. Although mwSuMD is an unbiased adaptive sampling method, it can be easily coupled to many forms of bias to favor the simulation of energy-requiring processes. Our results suggest a fundamental contribution of the phospholipids on the stabilization of the heterodimer, in agreement with experiments^{71, 72} but confuting protein-protein molecular docking results frequently predicting extended interfaces between monomers⁷³. mwSuMD was able to dynamically dock the heterobivalent ligand CP26, supporting a stabilizing effect on the A_2AR:D₂R heterodimer. A complete characterization of the possible interfaces between GPCR monomers, which falls beyond the goal of the present work, should be achieved by preparing different initial unbound states characterized by divergent relative orientations between monomers to dynamically dock.

mwSuMD is suited for a mechanistic description of structural events. For this reason, it well integrates with mutagenesis and kinetics experiments informed by the method, and required for the validation of its results in line with any other computational technique. We note that mwSuMD trajectories, since can describe the sequential states along a transition pathway, can represent a precious backbone for further MD sampling aimed at quantifying or predict the kinetics of the event. Approaches such as path collective variable metadynamics⁷⁴, Markov State Models⁷⁵, or machine learning models^{76, 77} can be informed by mwSuMD. Future work will be addressed in this direction.

In summary, we showcased the applicability domain of mwSuMD to key, but scarcely understood, aspects of GPCR structural biology, pharmacology, and drug design hitherto unaddressed by unbiased simulations. Given the generality and simplicity of its implementation, we anticipate that mwSuMD can be employed to study a wide range of structural phenomena characterizing both membrane and cytosolic proteins. mwSuMD is an ongoing project updated on https://github.com/pipitoludovico/mwSuMD_multiEngine and in the recent implementation can exploit ACEMD⁷⁸, NAMD⁷⁹, GROMACS⁸⁰ and OPENMM⁸¹ as MD engine.

Methods

Multiple walker SuMD (mwSuMD) protocol

The supervised MD (SuMD) is an adaptive sampling method⁸² for speeding up the simulation of binding events between small molecules (or peptides^{83, 84}) and proteins^{1, 18} without the introduction of any energetic bias. Briefly, during SuMD, a series of short unbiased MD simulations are performed, and after each simulation, the distances between the centers of mass (or the geometrical centers) of the ligand and the predicted binding site (collected at regular time intervals) are fitted to a linear function. If the resulting slope is negative (showing progress towards the target) the next simulation step starts from the last set of coordinates and velocities, otherwise, the simulation is restarted by randomly assigning the atomic velocities.

In the implementation for ACEMD used in this work, mwSuMD needs as input the initial coordinates of the system as a pdb file, the coordinates, and the atomic velocities of the system from the equilibration stage, the topology file of the system, and all the necessary force field parameters. The user can decide to supervise one (X) or two metrics (X’, X’’) of the simulated system over short simulations seeded in batches, called walkers. In the former case, either the slope of the linear function interpolating the metric values or a score can be adopted to decide whether to continue the mwSuMD simulation. When the user decides to supervise two metrics, then a specific score is used. In the present work, distances between centroids, RMSDs, or the number of atomic contacts between two selections were supervised (Table S1). The choice of the metrics is system and problem-dependent, as the RMSD is most useful when the final state is known, while the distance is required when the target state is unknown; details on the scores are given below. The decision to restart or continue mwSuMD after any short simulation is postponed until all the walkers of a batch are collected. The best short simulation is selected and extended by seeding the same number of walkers, with the same duration as the step before.

For each walker, the score for the supervision of a single metric (SMscore) is computed as the square root of the product between the metric value in the last frame (X_{last frame}) and the average metric value over the short simulation (X̅):

If the metric is set to decrease (e.g. binding or dimerization) the walker with the lowest SMscore is continued, otherwise (e.g. unbinding or outwards opening of domains), the walker with the highest score is continued. Using the SMscore rather than the slope should give more weight to the final state of each short simulation, as it is the starting point for the successive batch of simulations. Considering the average of the metric should favor short simulations consistently evolving in the desired direction along the metric.

If both X’ and X’’ are set to increase during the mwSuMD simulations, the score for the supervision of two metrics (DMscore) on each walker is computed as follows:

Where X’_{last frame} and X’’_{last frame} are the metrics values in the last frame, while X̅ ^’_{batch walkers} and X̅ ^’’_{batch walkers} represent the average value of the two metrics over all the walkers in the batch. Subtracting the value 1 to the metric ratio ensures that if one of the two metrics from the last frame (X’_{last frame} or X^’’_{last frame}) is equal to the average (X̅ ^’’_{batch walkers} or X̅ ^’’_{batch walkers}) then that metric addend is null and DMscore depends only on the remaining metric. If any of the two metrics is set to decrease, then the corresponding component in Equation 2 is multiplied by −1 to maintain a positive score. Considering the average value of the two metrics over all the walkers rather than only over the considered walker should be more representative of the system evolution along the defined metric. In other words, the information about the metric is taken from all the walkers to better describe the evolution of the system.

The DMScore is designed to preserve some degree of independence between the two metrics supervised. Indeed, if the variation of one of them slows down and gets close to zero, the other metric is still able to drive the system’s evolution. It should be noted that DMScore works at its best if the two metrics have similar variations over time, as in the case of distance and RMSD (both of which are distance-based). Notably, differently from SuMD, when a walker is extended by seeding a new batch of short simulations and the remaining walkers are stopped, the atomic velocities are not reassigned. This allows the simulations to be as short as a few picoseconds if desired, without introducing artifacts due to the thermostat latency to reach the target temperature (usually up to 10-20 ps when a simulation is restarted reassigning the velocities of the atoms).

The current implementation of mwSuMD is for python3 and exploits MDAnalysis⁸⁵ and MDTRaj⁸⁶ modules.

Force field, ligands parameters, and general systems preparation

The CHARMM36^{87, 88}/CGenFF 3.0.1^89–91 force field combination was employed in this work. Initial ligand force field, topology and parameter files were obtained from the ParamChem webserver⁸⁹. Restrained electrostatic potential (RESP)⁹² partial charges were assigned to all the non-peptidic small molecules but adrenaline and guanosine-5’-diphosphate (GDP) using Gaussian09 (HF/6-31G* level of theory) and AmberTools20.

Six systems were prepared for MD (Table S1). Hydrogen atoms were added using the pdb2pqr⁹³ and propka⁹⁴ software (considering a simulated pH of 7.0); the protonation of titratable side chains was checked by visual inspection. The resulting receptors were separately inserted in a 1-palmitoyl-2-oleyl-sn-glycerol-3-phosphocholine (POPC) bilayer (previously built by using the VMD Membrane Builder plugin 1.1, Membrane Plugin, Version 1.1. at: http://www.ks.uiuc.edu/Research/vmd/plugins/membrane/), through an insertion method⁹⁵. Receptor orientation was obtained by superposing the coordinates on the corresponding structure retrieved from the OPM database⁹⁶. Lipids overlapping the receptor transmembrane helical bundle were removed and TIP3P water molecules⁹⁷ were added to the simulation box by means of the VMD Solvate plugin 1.5 (Solvate Plugin, Version 1.5. at <http://www.ks.uiuc.edu/Research/vmd/plugins/solvate/). Finally, overall charge neutrality was reached by adding Na⁺/Cl^- counter ions up to the final concentration of 0.150 M), using the VMD Autoionize plugin 1.3 (Autoionize Plugin, Version 1.3. at <http://www.ks.uiuc.edu/Research/vmd/plugins/autoionize/).

System equilibration and general MD settings

The MD engine ACEMD3⁷⁸ was employed for both the equilibration and productive simulations. The equilibration was achieved in isothermal-isobaric conditions (NPT) using the Berendsen barostat⁹⁸ (target pressure 1 atm) and the Langevin thermostat⁹⁹ (target temperature 300 K) with low damping of 1 ps^-1. For the equilibration (integration time step of 2 fs): first, clashes between protein and lipid atoms were reduced through 1500 conjugate-gradient minimization steps, then a positional constraint of 1 kcal mol^-1 Å^-2 on all heavy atoms was gradually released over different time windows: 2 ns for lipid phosphorus atoms, 60 ns for protein atoms other than alpha carbon atoms, 80 ns for alpha carbon atoms; a further 20 ns of equilibration was performed without any positional constraints.

Productive trajectories (Table S1) were computed with an integration time step of 4 fs in the canonical ensemble (NVT). The target temperature was set at 300 K, using a thermostat damping of 0.1 ps^-1; the M-SHAKE algorithm^{100, 101} was employed to constrain the bond lengths involving hydrogen atoms. The cut-off distance for electrostatic interactions was set at 9 Å, with a switching function applied beyond 7.5 Å. Long-range Coulomb interactions were handled using the particle mesh Ewald summation method (PME)¹⁰² by setting the mesh spacing to 1.0 Å.

Vasopressin binding simulations

The vasopressin 2 receptor (V₂R) in complex with vasopressin (AVP) and the G_s protein¹⁰³ was retrieved from the Protein Data Bank¹⁰⁴ (PDB ID 7DW9). The G_s was removed from the system and the missing residues on ECL2 (G185-G189) were modeled from scratch using Modeller 9.19¹⁰⁵. AVP was placed away from V₂R in the extracellular bulk and the resulting system was prepared for MD simulations and equilibrated as reported above.

During SuMD simulations, the distance between the centroids of AVP residues C1-Q4 and V₂R residues Q96, Q174, Q291, and L312 (Cα atoms only) was supervised over time windows of 600 ps or 100 ps (Table S1). MwSuMD simulations considered the same distance, the RMSD of AVP residues C1-Q4 to the experimental bound complex or the combination of the two during time windows of 600 ps (3 walkers) or 100 ps (10 walkers) (Table S1). Slope, SMscore, or DMscore (see Methods section MwSuMD protocol) was used in the different mwSuMD replicas performed (Table S1). Simulations were stopped after 300 ns (time window duration = 600 ps) or 50 ns (time window duration = 100 ps) of total SuMD or mwSuMD simulation time.

Vasopressin unbinding simulations

The V₂R:AVP complex was prepared for MD simulations and equilibrated as reported above. During both SuMD and mwSuMD simulations (Table S1), the distance between the centroids of AVP residues C1-Q4 and V₂R residues Q96, Q174, Q291, and L312 (Cα atoms only) was supervised over time windows of 100 ps (10 walkers seeded for mwSuMD simulations). Replicas were stopped when the AVP-V₂R distance reached 40 Å.

G_s protein:β₂ AR binding simulations

The model of the adrenergic β₂ receptor (β₂ AR) in an intermediate active state was downloaded from GPCRdb (https://gpcrdb.org/). The full agonist adrenaline (ALE) was inserted in the orthosteric site by superposition with the PDB ID 4LDO (fully active β₂ AR)¹⁰⁶. The structure of the inactive, GDP-bound G_s protein¹⁰⁷ was retrieved from the Protein Data Bank¹⁰⁴ (PDB ID 6EG8) and placed in the intracellular bulk. The resulting system (G_s > 50 Å away from (β₂ AR) was prepared for MD simulations and equilibrated as reported above. The PDB ID 3SN6 (fully active β₂ AR in complex with G_s⁴¹) was used as the reference for RMSD computations. Three mwSuMD replicas (Table S1) were performed supervising at the same time the distance between the helix 5 (α5) G_αs residues R385-L395 and the β₂ AR residues V31-P330 as well as the RMSD of β₂ AR TM6 residues C265-I278 (Cα atoms only) to the fully active state, during 100 ps time windows (5 walkers).

Membrane-anchored Gi protein:A₁R simulations

Since the full-length structure of the inactive human Gi protein has not been yet resolved by X-ray or cryo-EM, it was modeled by superimposing the AlphaFold2¹⁰⁸ AI models of the Gαi (P63096-F1), Gβ (Q9HAV0-F1), and Gγ (P50151-F1) subunits to the PDB file 6EG8 (a G_s heterotrimer). The resulting homotrimer (without GDP) was processed through Charmm-GUI¹⁰⁹ to palmitoylate residue C3^Gαi and geranylgeranylate residue C65^{Gγ 44, 110}. The side chains of these two lipidated residues were manually inserted into a 120 x 120 Å POPC membrane and previously built by using the VMD Membrane Builder plugin 1.1, Membrane Plugin, Version 1.1. at: http://www.ks.uiuc.edu/Research/vmd/plugins/membrane/. Lipids overlapping the palmitoyl and geranylgeranyl groups were removed and TIP3P water molecules⁹⁷ were added to the simulation box by means of the VMD Solvate plugin 1.5 (Solvate Plugin, Version 1.5. at <http://www.ks.uiuc.edu/Research/vmd/plugins/solvate/). Finally, overall charge neutrality was reached by adding Na⁺/Cl^- counter ions up to the final concentration of 0.150 M), using the VMD Autoionize plugin 1.3 (Autoionize Plugin, Version 1.3. at <http://www.ks.uiuc.edu/Research/vmd/plugins/autoionize/). The first stage of equilibration was performed as reported above (Methods section System equilibration and general MD settings) for 120 ns, followed by a second stage in the NVT ensemble for a further 1 μs without any restraints to allow the membrane-anchored heterotrimeric Gi protein to stabilize within the intracellular side of the simulation box. After this two-stage, long equilibration, the active state A₁R in complex with adenosine (PDB 6D9H) was manually inserted into the equilibrated membrane above the Gi protein using the corresponding structure retrieved from the OPM database as a reference, and the system further equilibrated for 120 ns as reported above (Methods section System equilibration and general MD settings). The A₁R-Gi system was then subjected to both a 1 μs-long classic MD simulation and a mwSuMD simulation (Table S1). During the mwSuMD simulation, the RMSD of helix 5 (α5) G_αs residues 329-354 to the PDB 6D9H was supervised, seeding three walkers of 100 ps each until the productive simulation time reached 50 ns (total simulation time 150 ns).

GLP-1R:PF06882961 binding simulations

The inactive, ligand-free glucagon-like peptide receptor (GLP-1R) was retrieved from the Protein Data Bank¹⁰⁴ (PDB ID 6LN2)¹¹¹. Missing residues in the stalk and ICL2 were modeled with Modeller 9.19. The PF06882961 initial conformation was extracted from the complex with the fully active GLP-1R⁵⁴ (PDB ID 7LCJ) and placed away from GLP-1R in the extracellular bulk. The resulting system was prepared for MD simulations and equilibrated as reported above. CGenFF dihedral force field parameters of PF06882961 with the highest penalties (dihedrals NG2R51-CG321-CG3C41-CG3C41 (penalty=143.5) and NG2R51-CG321-CG3C41-OG3C51 (penalty=152.4)) were optimized (Figure S13) employing Gaussian09 (geometric optimization and dihedral scan at HF/6-31g(d) level of theory) and the VMD force field toolkit plugin¹¹².

Four classic MD replicas, for a total of 8 μs, were performed on the inactive, ligand-free receptor (prepared for MD simulations and equilibrated as reported above) to assess the possible binding path to the receptor TMD and therefore decide the initial position of PF06882961 in the extracellular bulk of the simulation box. A visual inspection of the trajectories suggested three major conformational changes that could allow ligand access to the TMD (Figure S13). Transitory openings of the ECD (distance Q47^ECD - S310^ECL2), TM6- TM7 (distance H363^6.52 - F390^7.45), and TM1-ECL1 (distance E138^1.33 and W214^ECL1) were observed. Since the opening of TM1-ECL1 was observed in two replicas out of four, we placed the ligand in a favorable position for crossing that region of GLP-1R.

MwSuMD simulations (Table S1) were performed stepwise to dock the ligand within GLP-1R first and then relax the receptor towards the active state. The PF06882961 binding was obtained by supervising at the same time the distance between the ligand and GLP-1R TM7 residues L379^7.34-F381^7.36, which are part of the orthosteric site (Cα atoms only), and the RMSD of the ECD (residues W33^ECD-W120^ECD, Cα atoms only) to the active state (PDB ID 7LCJ) until the former distance reached 4 Å. In the second phase of mwSuMD, the RMSD of the ECD (residues W33^ECD-W120^ECD, Cα atoms only) and the ECL1 to the active state (PDB ID 7LCJ) Cα atoms of residues M204^2.74-L224^3.27) were supervised until the latter reached less than 4 Å. During the third phase, the RMSD of PF06882961, as well as the RMSD of ECL3 (residues A368^6.57-T378^7.33, Cα atoms), were supervised until the former reached values lower than 3 Å. In the last mwSuMD step, only the RMSD of TM6 (residues I345^6.34-F367^6.56, Cα atoms) to the active state (PDB ID 7LCJ) was supervised until less than 5 Å.

Membrane-anchored G_s protein:GLP-1R simulations and GDP dissociation

The PDB 6EG8 was processed through Charmm-GUI¹⁰⁹ to palmitoylate residue C3^Gαi. The Gβ and Gγ subunits prepared as reported above were superimposed to the PDB file 6EG8 and the original Gβ and Gγ were removed. This step produced a heterotrimeric G_s protein palmitoylated at residue C3^Gαs and geranylgeranylated at residue C65^Gγ. The resulting system was inserted into a 120 x 120 Å POPC membrane and previously built by using the VMD Membrane Builder plugin 1.1, Membrane Plugin, Version 1.1. at: http://www.ks.uiuc.edu/Research/vmd/plugins/membrane. Lipids overlapping the palmitoyl and geranylgeranyl groups were removed and TIP3P water molecules⁹⁷ were added to the simulation box by means of the VMD Solvate plugin 1.5 (Solvate Plugin, Version 1.5. at <http://www.ks.uiuc.edu/Research/vmd/plugins/solvate/). Finally, overall charge neutrality was reached by adding Na⁺/Cl^- counter ions up to the final concentration of 0.150 M), using the VMD Autoionize plugin 1.3 (Autoionize Plugin, Version 1.3. at <http://www.ks.uiuc.edu/Research/vmd/plugins/autoionize/). The first stage of equilibration was performed as reported above (Methods section System equilibration and general MD settings) for 120 ns, followed by a second stage in the NVT ensemble for a further 1 μs without any restraints to allow the membrane-anchored heterotrimeric G_s protein to stabilize within the intracellular side of the simulation box. After this two-stage long equilibration, GLP-1R from the final frame of the activation simulation (in complex with PF06882961) was manually inserted into the equilibrated membrane above the G_s protein using the corresponding structure retrieved from the OPM database as a reference, and the system further equilibrated for 120 ns as reported above (Methods section System equilibration and general MD settings). The GLP-1R-G_s system was then subjected to three simulations (Table S1). Each mwSuMD replica was interrupted by 500 ns of classic MD twice, to relax the system during the transition. In the supervised stages, the distance between residues M386-L394 ^Gαs of helix 5 (α5) and the GLP-1R intracellular residues R176^2.46, R348^6.37, S352^6.41, and N405^7.60 was monitored, seeding three walkers of 200 ps each.

The AHD opening was simulated starting from the best GLP-1R:G_s binding mwSuMD replica (Replica 2 in Figure 4) by supervising the distance between AHD residues G70-R199 ^Gαs and K300-L394^Gαs during three walkers of 100 ps each. 300 ns of classic MD was performed to relax the system. Finally, the GDP unbinding was supervised as the distance between GDP and residues E50, K52, T55, K293, and V367 of G_αs; five walkers were used in a 50 ps long mwSuMD simulation.

A2A:D2R heterodimerization

The inactive state A_2AR and D₂R were retrieved from the Protein Data Bank¹⁰⁴ (PDB 5NM4 and 6LUQ, respectively)^{113, 114}. Antagonists bound to the orthosteric site were removed and no modeling of the missing IC loops was attempted. A_2AR and D₂R were manually placed roughly 40 Å away from each other, on the plane of the membrane, orienting the two receptors to favor the dimerization through the interface formed by TM4 and TM5, as suggested by Borroto-Esquela D. O. et al.⁶² The resulting system was prepared for MD simulations and equilibrated as reported above.

The heterodimerization between A_2AR and D₂R was simulated with mwSuMD, seeding batches of three walkers with a duration of 100 ps each (Table S1). During each walker, the distance between TM5 of A_2AR and D₂R was supervised. At the same time, the distance between the centroids of A_2AR and D₂R was used as a collective variable for adiabatic MD¹¹⁵ (abMD) and well-tempered metadynamics^{116, 117} (wtMetaD) performed with Plumed 2.6¹¹⁸. For abMD, a distance target of 30 Å and a force constant of 10000 kJ*mol^-1*Å^-1) was used, while wtMetaD was performed by seeding Gaussian functions every 1 ps (sigma=1 Å; height=0.837 kJ*mol^-1; T=310K) with a bias factor of 30. When the A_2AR - D₂R distance reached values lower than 40 Å and the first contacts between proteins were formed, the abMD was stopped and wtMetaD continued with a harmonic energy wall at 30 Å to avoid artificial crushing between the receptors due to the added energy bias. When the distance between A_2AR and D₂R was stable at about 30 Å, the collective variable biased by wtMetaD was set as the number of atomic contacts between A_2AR and D₂R, until reaching 200 ns of simulation. Finally, to relax the system and challenge the stability of the heterodimer formed during the biased mwSuMD simulation, a 1.5 μs classic MD simulation was performed.

A_2AR-D₂R heterobitopic ligand binding simulations

The A₂AR-D₂R heterobivalent ligand compound 26³³ was parameterized as reported above and placed in the bulk solvent of the A₂AR:D₂R complex from the classic MD. Four mwSuMD replicas were collected supervising at the same time the distance between the A_2A antagonist pyrazole-triazole-pyrimidine scaffold and the centroid of A_2AR residues F168, N253, and A277 (Cα atoms) as well as the distance between the D₂ antagonist 4-fluorobenzyl scaffold and the centroids of the Cα of D₂R residues C118, F198, and V115 (Cα atoms). Ten walkers of 100 ps were simulated for every mwSuMD batch of replicas.

MD Analysis

Interatomic distances were computed through MDAnalysis⁸⁵; root mean square deviations (RMSD) were computed using VMD¹¹⁹ and MDAnalysis⁸⁵.

Interatomic contacts and ligand-protein hydrogen bonds were detected using the GetContacts scripts tool (https://getcontacts.github.io), setting a hydrogen bond donor-acceptor distance of 3.3 Å and an angle value of 120° as geometrical cut-offs. Contacts and hydrogen bond persistency are quantified as the percentage of frames (over all the frames obtained by merging the different replicas) in which protein residues formed contacts or hydrogen bonds with the ligand.

The MMPBSA.py¹²⁰ script, from the AmberTools20 suite (The Amber Molecular Dynamics Package, at http://ambermd.org/), was used to compute molecular mechanics energies combined with the generalized Born and surface area continuum solvation (MM/GBSA) method or the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) approach, after transforming the CHARMM psf topology files to an Amber prmtop format using ParmEd (documentation at <http://parmed.github.io/ParmEd/html/index.html). Supplementary Videos were produced employing VMD and avconv (at https://libav.org/avconv.html). Molecular graphics images were produced using the UCSF Chimera¹²¹ (v1.14).

Numbering system

Throughout the manuscript, the Ballesteros-Weinstein residues numbering system for class A¹²² and the Wootten residues numbering system for class B GPCRs¹²³ are adopted.

Associated content

Supporting Information (PDF)

Supporting Videos S1-S8 (mpg).

All the MD trajectories (stripped of POPC, water molecules, and ions) and topology files (psf and pdb) are available here:

https://zenodo.org/record/7944479

The mwSuMD software version used in this study is available at:

https://github.com/pipitoludovico/mwSuMD

Author Contributions

GD and CAR conceived and supervised the project; GD designed and implemented the software and planned the simulations; GD, LP, RMR, TW, and PG carried out the simulations; GD analyzed the data; GD, AC, SM, and CAR interpreted the results, GD wrote the manuscript with input from CAR, AC, and SM; all the authors edited and reviewed the final version of the manuscript.

Acknowledgements

GD is a member of the GPCRs-focused European COST action ERNEST. CAR is grateful for a Royal Society Industry Fellowship. GD and CAR are grateful for support from the BBSRC (BB/W016974/1) and from Diabetes UK (BDA 20/0006307).

Competing Interest

All authors declare that they have no conflicts of interest.