Deep Learning Dynamic Allostery of G-Protein-Coupled Receptors

G-protein-coupled receptors (GPCRs) are the largest superfamily of human membrane proteins and represent primary targets of ~ 1/3 of currently marketed drugs. Allosteric modulators have emerged as more selective drug candidates compared with orthosteric agonists and antagonists. However, many X-ray and cryo-EM structures of GPCRs resolved so far exhibit negligible differences upon binding of positive and negative allosteric modulators (PAMs and NAMs). Mechanism of dynamic allosteric modulation in GPCRs remains unclear. In this work, we have systematically mapped dynamic changes in free energy landscapes of GPCRs upon binding of allosteric modulators using the Gaussian accelerated molecular dynamics (GaMD), Deep Learning (DL) and free energy prOfiling Workflow (GLOW). A total of 18 available high-resolution experimental structures of allosteric modulator-bound class A and B GPCRs were collected for simulations. A number of 8 computational models were generated to examine selectivity of the modulators by changing their target receptors to different subtypes. All-atom GaMD simulations were performed for a total of 66 μs on 44 GPCR systems in the presence/absence of the modulator. DL and free energy calculations revealed significantly reduced conformational space of GPCRs upon modulator binding. While the modulator-free GPCRs often sampled multiple low-energy conformational states, the NAMs and PAMs confined the inactive and active agonist-G protein-bound GPCRs, respectively, to mostly only one specific conformation for signaling. Such cooperative effects were significantly reduced for binding of the selective modulators to “non-cognate” receptor subtypes in the computational models. Therefore, comprehensive DL of extensive GaMD simulations has revealed a general dynamic mechanism of GPCR allostery, which will greatly facilitate rational design of selective allosteric drugs of GPCRs.


Introduction
G-protein-coupled receptors (GPCRs) are the largest superfamily of human membrane proteins with > 800 members.GPCRs play key roles in cellular signaling and mediate various physiological activities, including vision, olfaction, taste, neurotransmission, endocrine, and immune responses 1 .They represent primary targets of ~ 1/3 of currently marketed drugs 2 .GPCRs can be classi ed into six different classes, including class A (Rhodopsin-like), B (secretin receptors), C (metabotropic glutamate receptors (mGluRs)), D (fungal mating pheromone receptors), E (cyclic AMP receptors), and class F (frizzled/TAS2 receptors) 3,4 .GPCRs share a characteristic structural fold of seven transmembrane (TM) helices (TM1-TM7) connected by three extracellular loops (ECL1-ECL3) and three intracellular loops (ICL1-ICL3).For decades, the primary endogenous agonist-binding ("orthosteric") site has been targeted for drug design of GPCR agonists, antagonists and inverse agonists 5 .However, the orthosteric site is usually highly conserved in different subtypes of GPCRs.An orthosteric drug often binds and activates/deactivates multiple GPCRs simultaneously with poor selectivity, thereby causing toxic side effects 6 .Alternatively, allosteric modulators have been discovered to bind topographically distant ("allosteric") sites of GPCRs with advantages [7][8][9][10][11][12] .They are able to modulate the binding a nity and signaling of orthosteric ligands, including positive and negative allosteric modulators (PAMs and NAMs) 13 .The allosteric effect has been shown to depend on the orthosteric probe 14 , with a "ceiling level" determined by the magnitude and direction of cooperativity between the orthosteric and allosteric ligands.Because the allosteric site is usually more divergent in residue sequences and conformations, allosteric modulators offer higher receptor selectivity than the orthosteric ligands.They serve as important chemical probes and promising selective therapeutics of GPCRs.
Important insights have been obtained using X-ray crystallography and cryo-electron microscopy (cryo-EM) about structural changes induced by allosteric modulator binding in certain GPCRs 7,9,10,15 .For class A GPCRs, binding of the LY2119620 PAM in the M 2 muscarinic receptor (M 2 R) led to sidechain rotation of residue W7. 35 and slight contraction of the receptor extracellular pocket, which was pre-formed in the active agonist-bound structure 16,17 .GPCR residues are numbered according to the Ballesteros-Weinstein scheme 18 .Binding of the muscarinic toxin MT7 NAM to the antagonist bound M 1 muscarinic receptor (M 1 R) resulted in conformational changes in the ECL2, TM1, TM2, TM6 and TM7 extracellular domains, as well as the TM2 and TM6 intracellular domains 19 .In the free fatty acid receptor GPR40 (FFAR1), AgoPAM binding in a lipid-facing pocket formed by TM3-TM4-ICL2 induced conformational changes in the ICL2, TM4 and TM5 of the active receptor 20 .The ICL2 adopted a short helical conformation and the TM5 was shifted along its helical axis towards the extracellular side relative to the TM4 20 .A similar allosteric site was identi ed for binding of the NDT9513727 and Avacopan NAMs between TM3-TM4-TM5 on the lipid-exposed surface of the C5a 1 receptor (C5AR1) 21 .For class B GPCRs, the LSN3160440 PAM was found to bind between the extracellular domains of TM1 and TM2 of the GLP-1 receptor (GLP1R) 22 .In the glucagon receptor (GLR), NAM binding outside of the 7TM bundle between TM6-TM7 restricted the outward movement of the TM6 intracellular domain required for activation and G-protein coupling of the receptor 23 .The ECL2 stretched to the central axis of the TM helical bundle, allowing for interactions from TM3 to TM6 and TM7 in the inactive class B GPCRs 23 .Despite remarkable advances, the X-ray and cryo-EM structures represent static snapshots of GPCRs.Many GPCR structures exhibit small/negligible differences in the absence and presence of allosteric modulators, notably for the A 1 AR 24 , M 4 R 25,26 , β 2 -adrenoceptor (β 2 AR) [27][28][29] , C5AR1 30 , CB 1 cannabinoid receptor (CB 1 ) 31 , chemokine receptor CCR2 32 , dopamine receptor 1 (D 1 R) 33 , GPBA receptor (GPBAR) 34 , and GLP1R 22,35,36 .A dynamic review has been suggested for allosteric modulation of GPCRs 37 .However, the dynamic mechanism of GPCR allostery remains unclear.
Molecular dynamics (MD) is a powerful computational technique for simulating biomolecular dynamics on an atomistic level 38 .For GPCRs, MD has been applied to simulate binding of both orthosteric and allosteric ligands [39][40][41][42][43][44] .Binding of known NAMs to the M 2 R was observed in conventional MD (cMD) simulations using the specialized supercomputer Anton 45 .The modulators formed cation-π interactions with aromatic residues in the receptor extracellular vestibule, which was con rmed by mutation experiments and later by X-ray structure of M 2 R recognized by a PAM 16 .Microsecond-timescale cMD simulations revealed mechanistic insights into allosteric modulation by Na + in dopamine and opioid receptors 46,47 .Accelerated MD (aMD) simulations also captured Na + binding to the highly conserved D2.50 allosteric site, which stabilized a muscarinic GPCR in the inactive state 48 .Recently, spontaneous binding of prototypical PAMs to the putative ECL2 allosteric site of the A 1 AR was captured in Gaussian accelerated molecular dynamics (GaMD) simulations 49 .Moreover, metadynamics simulations captured binding of the BMS-986187 PAM to the δ-opioid receptor 50 .Metadynamics and GaMD enhanced sampling simulations revealed positive binding cooperativity between allosteric and orthosteric ligands of the CCR2 51 .Despite these exciting advances, MD simulations of allosteric modulation have been limited to mostly few selected class A GPCRs 44 .
Recently, we have developed the GaMD, Deep Learning (DL) and free energy pro ling work ow (GLOW) to predict molecular determinants and map free energy landscapes of biomolecules 52 .GaMD is an unconstrained enhanced sampling technique that works by applying a harmonic boost potential to smooth biomolecular potential energy surface 53 .Since this boost potential usually exhibits a near Gaussian distribution, cumulant expansion to the second order ("Gaussian approximation") can be applied to achieve proper energy reweighting 54 .GaMD allows for simultaneous unconstrained enhanced sampling and free energy calculations of large biomolecules 53 .GaMD has been successfully demonstrated on enhanced sampling of ligand binding, protein folding, protein conformational changes, protein-membrane/peptide/protein/nucleic acid/carbohydrate interactions 55 .In GLOW, DL of imagetransformed residue contact maps calculated from GaMD simulation frames allows us to identify important residue contacts by classic gradient-based pixel attribution in the saliency (attention) maps 52,56 .Finally, free energy pro les of these residue contacts are calculated through reweighting of GaMD simulations to characterize the biomolecular systems of interest 52 .
In this work, we have applied GLOW to systematically map dynamic changes in free energy landscapes of GPCRs upon binding of allosteric modulators (Fig. 1).A total of 18 different high-resolution experimental structures of class A and B GPCRs are collected for modeling and 8 computational models were generated by changing target receptors of the modulators to different subtypes.Our comprehensive DL analysis of extensive GaMD simulations has provided important mechanistic insights into the dynamic allostery of GPCRs.

Simulation Protocols
All-atom dual-boost GaMD simulations 53 were performed on the GPCR structures and models with and without allosteric modulators (Fig. 1C and Supplementary Table 1).The simulations of the A 1 AR with and without the MIPS521 PAM were obtained from a previous study 24 .GaMD simulations of the other systems followed a similar protocol.The systems were energetically minimized for 5000 steps using the steepest-descent algorithm and equilibrated with the constant number, volume, and temperature (NVT) ensemble with 310 K.They were further equilibrated for 375 ps at 310 K with the constant number, pressure, and temperature (NPT) ensemble.The cMD simulations were then performed for 10 ns using the NPT ensemble with constant surface tension at 1 atm pressure and 310 K temperature.GaMD implemented in GPU version of AMBER 20 53,70,71 was applied to simulate the GPCR systems.The simulations involved an initial short cMD of 4-10 ns to calculate GaMD acceleration parameters and GaMD equilibration of added boost potential for 16-40 ns.Three independent 500-ns GaMD production simulations with randomized initial atomic velocities were performed for each system at the "dual-boost" level, with one boost potential applied to the dihedral energetic term and the other to the total potential energetic term.The reference energy was set to the lower bound E = V max , and the upper limit of the boost potential standard deviation, σ 0 , was set to 6.0 kcal/mol for both the dihedral and total potential energetic terms.The GaMD simulations were summarized in Supplementary Table 1.

Deep Learning And Free Energy Pro ling Of GaMD Simulations With
The GLOW Work ow GLOW 52 was applied to systematically analyze GPCR allostery (Fig. 1D-1E).Residue contact maps of 6,600,000 GaMD simulation frames obtained from 66µs GaMD simulations were calculated and transformed into images for DL (Fig. 1D).A contact de nition of 4.5 Å between any heavy atoms in two protein residues was used.For DL, 80% of the residue contact map images were randomly assigned to the training set, while the remaining 20% were put in the validation set for each GPCR.The residue contact map images of each GPCR system were separated into four different classes for DL analysis based on the absence and presence of NAMs and PAMs, including NAM-free ("Antagonist"), NAM-bound ("AntagonistNAM"), PAM-free ("Agonist"), and PAM-bound ("AgonistPAM").DL models of two-dimensional (2D) convolutional neural networks were built to classify the frame images with and without the allosteric modulator bound for each GPCR subfamily.Saliency (attention) maps of residue contact gradients were calculated through backpropagation by vanilla gradient-based pixel attribution 56 using the residue contact map of the most populated structural cluster of each GPCR system (Fig. 1E).The hierarchical agglomerative clustering algorithm was used to cluster snapshots of receptor conformations with all GaMD production simulations combined for each system 52 .
Furthermore, root-mean-square uctuations (RMSFs) of the receptors and orthosteric ligands within the GPCR complexes were calculated by averaging the RMSFs calculated from individual GaMD simulations of each GPCR system.Changes in the RMSFs (∆RMSF) upon binding of allosteric modulators were calculated by subtracting the RMSFs of GPCRs without modulators from those with modulators bound (Fig. 1F).If the absolute average of ∆RMSF calculated from three simulations of a residue was smaller than the corresponding standard deviation, the exibility change for that residue was considered not signi cant and related residue pairs were neglected for further analysis.Important residue contacts were selected with the highest contact gradients in the attention maps from DL and signi cant changes in the GPCR residue exibility upon modulator binding (Fig. 1G-1H).They were nally used as reaction ≤ coordinates (RCs) to calculate free energy pro les by reweighting the GaMD simulations using the PyReweighting toolkit [52][53][54] , with bin sizes of 0.5-1.0Å and cutoff of 100-500 frames in one bin.

GaMD simulations on effects of allosteric modulator binding to GPCRs
All-atom GaMD simulations were obtained for systems of the A , GPBAR, GLP1R, and GLR.GaMD simulations performed in this study recorded overall similar averages (~ 13-16 kcal/mol) and standard deviations (~ 4-5 kcal/mol) of boost potentials across the different GPCR systems, except for the A 1 AR simulations that were obtained from a previous study 24 using a different force eld parameter set 69 (Supplementary Table 1).We rst examined structural dynamics of the GPCR orthosteric and allosteric ligands.Time courses of the orthosteric and allosteric ligand RMSDs relative to the simulation starting structures were plotted in Supplementary Figs.1-3.In most of the GPCR systems, the orthosteric ligands underwent similar structural deviations in the absence and presence of allosteric modulators during the GaMD simulations (Supplementary Figs.1-3).This was consistent with previous ndings that modulator binding mostly does not cause large changes in the X-ray and cryo-EM structures of the GPCRs 22,24−36 .However, a number of orthosteric ligands, including adenosine in A 1 AR and PE5 in GLR exhibited signi cantly smaller structural deviations in the presence of the MIPS521 PAM and MK-0893 NAM, respectively (Supplementary Figs.2A and 3D).
In general, binding of allosteric modulators reduced uctuations of the orthosteric ligands and GPCRs as shown in Figs.2-3 and Supplementary Figs.4-5.NAM binding primarily stabilized the allosteric binding pockets, with additional reduced exibility observed in the extracellular and/or intracellular domains of the receptors (Fig. 2 and Supplementary Fig. 4).In particular, binding of MT7 signi cantly reduced uctuations of the extracellular mouth between ECL2 and ECL3 in M 1 R (Fig. 2A).Binding of Cmpd-15 in the β 2 AR and GTPL9431 in the CCR2 reduced uctuations of the allosteric pocket formed by TM6, TM7, ICL4, and H8 (Fig. 2B, 2G).Notably, AS408, NDT9513727, and Avacopan bound to a similar TM3-TM4-TM5 region on the lipid-facing surface of the β 2 AR and C5AR1.They reduced uctuations of the intracellular domains of TM3, TM4, and TM5 (Fig. 2C-2E), as well as ICL2 in β 2 AR (Fig. 2C).This was consistent with recent nding that the NDT9513727 NAM stabilized TM5 through the hydrophobic stacking between TM4 and TM5 72 .Binding of NNC0640 to GLP1R and MK-0893 to GLR stabilized the lipid-facing pocket on the intracellular domains of TM7 and TM6, respectively (Fig. 2H, 2J and Supplementary Fig. 4H, 4J).While PF-06372222 also bound to a similar region in GLP1R, no signi cant exibility change was observed in the receptor likely because the modulator-free receptor was already stable (Fig. 2I and Supplementary Fig. 4I).Lastly, binding of ORG27569 to CB 1 reduced uctuations of the TM4, TM6, TM7 intracellular domains and ICL2 (Fig. 2F).
Binding of PAMs to GPCRs generally reduced uctuations of the extracellular domains, orthosteric agonist-binding pocket, and/or intracellular G-protein coupling domains (Fig. 3 and Supplementary Fig. 5).Speci cally, binding of MIPS521 to the A 1 AR signi cantly reduced uctuations of the adenosine agonist, the orthosteric pocket constituted by extracellular domains of TM2, ECL1, ECL2, TM3, TM5, TM6, ECL3, and TM7, and the intracellular domains of TM5, TM7, and ICL2 (Fig. 3A and Supplementary Fig. 5A), which was consistent with the previous observation that the PAM stabilized the adenosine-A 1 AR-G-protein complex 24 .LY2119620 binding to the M 2 R reduced exibility of TM2, TM7, ECL1-ECL3 and H8 (Fig. 3B and Supplementary Fig. 5B).In the M 4 R, LY2119620 binding signi cantly reduced exibility of the G-protein coupling domains of ICL1, TM3, TM5, TM6, and H8, while stabilized ECL3 to a lesser extent (Fig. 3C and Supplementary Fig. 5C).In the β 2 AR, Cmpd-6FA binding reduced uctuations of orthosteric residues in the ECL1 and TM7 extracellular end as well as the intracellular domains of ICL1, TM3, and ICL2 (Fig. 3D and Supplementary Fig. 5D).Binding of LY3154207 to the D 1 R reduced uctuations of the TM1 extracellular domain and ECL3, as well as the TM6 intracellular domain and ICL2 (Fig. 3E and Supplementary Fig. 5E).Similar to the 6N48 PDB structure of β 2 AR, the 5TZY PDB structure of AgoPAMbound FFAR1 did not have a G-protein bound.Binding of AgoPAM reduced exibility of the ECL1, ECL2, ECL3 and TM1, TM2, TM3, and TM5 extracellular ends as well as the ICL2, and TM3 and TM4 intracellular ends, which constituted the G-protein binding pocket (Fig. 3F and Supplementary Fig. 5F).The stabilization of ICL2 observed in the AgoPAM-bound FFAR1 was in good agreement with the structural data of the PAM-bound 5TZY and PAM-free 5TZR PDB structures 20 .In particular, the ICL2 adopted a short helical conformation in the PAM-bound 5TZY PDB and was completely missing in the PAM-free 5TZR PDB structure of FFAR1 20 .Two molecules of INT777 with − 1 charge served as the orthosteric and allosteric ligands of GPBAR, which potentially caused electrostatic repulsion and destabilized the orthosteric INT777 and most of the orthosteric residues (Fig. 3G and Supplementary Fig. 5G).In fact, Dror et al. uncovered that electrostatic repulsion between orthosteric and allosteric ligands could weaken the binding of one in the presence of the other 45 .Even so, binding of INT777 to the allosteric site reduced uctuations of ECL2 and the TM4 extracellular end, as well as the TM5 and TM6 intracellular ends (Fig. 3G and Supplementary Fig. 5G).Finally, the LSN3160440 PAM binding signi cantly reduced uctuations of the large orthosteric pocket of GLP1R formed by TM2, ECL1, ECL2, TM4, and ECL3, and the G-protein coupling domains in TM5 and TM6 (Fig. 3H and Supplementary Fig. 5H).

Deep Learning Important Residue Contacts Underlying Allosteric Modulation Of GPCRs
DL models were built from image-transformed residue contact maps calculated from GaMD simulations of the modulator-free and modulator-bound systems for each GPCR subfamily.Classi cation of GPCRs bound by "Antagonist", "AntagonistNAM", "Agonist", "AgonistPAM" was carried out with high accuracies on both the training and validation sets after 15 epochs for all GPCRs (Supplementary Figs. 6 and 7).The saliency (attention) residue contact maps of gradients were calculated and plotted in Supplementary Fig. 8 for NAM-bound GPCRs and Supplementary Fig. 9 for PAM-bound GPCRs.Next, the residue contacts which exhibited the highest contact gradients ( 0.7) in the attention maps from DL and signi cant changes in uctuations of the involved residues upon allosteric modulator binding were considered important for allosteric modulation of GPCRs (Figs. 2 and 3).

Free Energy Pro ling Of Important Residue Contacts In GPCR Allostery
2D free energy pro les were calculated from GaMD simulations for the important residue contacts identi ed from DL and structural exibility analyses of GPCRs (Figs. 4 and 5).Overall, binding of NAMs and PAMs reduced conformational space of the inactive antagonist-bound and active agonist-bound GPCRs, respectively.Moreover, in case the modulator-free GPCRs were able to sample multiple low-energy conformational states, binding of the NAMs and PAMs con ned the GPCR residue contacts to fewer lowenergy states (only 1 for most GPCRs) (Figs. 4 and 5).The residue distances/RMSDs identi ed at the energy minima of these conformational states are listed in detail in Supplementary Tables 3 and 4.
Binding of the MT7 NAM to M 1 R con ned the TM4 extracellular domain and ECL2 from two conformational states ("S1" and "S2") to only the "S1" state, in which the extracellular mouth adopted the more closed conformation (Fig. 4A).Cmpd-15 binding to the β 2 AR reduced the conformational space from three states ("S1"-"S3") to only the "S3" state, in which the intracellular pocket formed by TM6, TM7, and H8 adopted the more open conformation to accommodate the NAM (Figs. 4B and 2B).Similarly, AS408 binding to the β 2 AR con ned the TM4 intracellular end and ICL2 from two conformational states ("S1" and "S2") to only state "S1", in which the intracellular pocket formed by the TM4 intracellular end and ICL2 adopted a more open conformation for stable NAM binding (Figs. 4C and 2C).Binding of NDT9513727 to the C5AR1 reduced the number of low-energy conformational states from two to one (the "S1" state), where the allosteric pocket located between TM3 and TM4 intracellular ends as well as the middle of TM5 and TM6 adopted a more open conformation to accommodate the NAM (Figs. 4D and  2D).Avacopan-binding to the C5AR1 con ned the TM4 and TM5 extracellular domains from two conformational states to only state "S1", where the TM4 and TM5 extracellular ends adopted the more closed conformation to stabilize NAM binding (Figs. 4E and 2E).The hydrophobic stacking found between TM4 and TM5 extracellular ends was again in good agreement with previous nding by Xiaoli et al. 72 .Binding of MK-0893 to the GLR con ned the conformational space from three states ("S1"-S3") to only state "S1", in which the middle of TM5 and TM6 as well as TM6 and TM7 extracellular ends adopted the more closed conformations (Fig. 4J).In the cases of ORG27569-bound CB 1 , GTPL9431-bound CCR2, NNC0640-bound GLP1R, and PF-06372222-bound GLP1R, the modulator-free GPCRs already sampled only one low-energy conformational state (Fig. 4F-4I).
In PAM-bound GPCRs, binding of MIPS521 to the A 1 AR con ned the TM2 extracellular domain and ECL2 as well as the middle of TM6 and TM7 from three conformational states ("S1"-"S3") to only the "S1" state, in which the extracellular mouth and orthosteric agonist-binding pocket adopted the more closed conformation to stabilize agonist binding (Figs. 5A and 3A).This nding was highly consistent with previous experimental and computational studies of the A 1 AR allosteric modulation 24,52,73,74 , where the PAM was found to signi cantly stabilize the adenosine agonist.The modulator-free GPCR in case of LY2119620-bound M 2 R sampled only one low-energy conformational state (Fig. 5B).LY2119620-binding to the M 4 R reduced the conformational space of the TM1 intracellular end and H8 from three states ("S1"-"S3") to only state "S1", in which the G-protein-binding region adopted the more closed conformation to stabilize G-protein binding (Fig. 5C).Binding of Cmpd-6FA to the β 2 AR con ned the TM2 and TM7 extracellular ends as well as TM2 and TM4 intracellular ends from three states ("S1"-S3") to the "S1" state, in which G protein-binding domain adopted a closed conformation (Fig. 5D).Binding of the LY3154207 PAM to the D 1 R con ned the TM1 extracellular end and ICL2 from two conformational states ("S1" and "S2") to only state "S1", reducing the exibility of both domains (Figs. 5E and 3E).AgoPAMbinding to the FFAR1 reduced the conformational space of the TM5 extracellular end, ECL3, and ICL2 from two states ("S1" and "S2") to state "S1", where the extracellular mouth between TM5 and ECL3 as well as the ICL2 adopted a closed and short helix conformation, respectively (Figs. 5F and 3F).Here, the nding that the ICL2 preferred a short helix conformation in the AgoPAM-bound FFAR1 resembled the structural data of 5TZY and 5TZR PDB structures well 20 .Binding of INT777 to the allosteric site of GPBAR reduced the number of low-energy conformational states from three ("S1"-"S3") to two ("S1" and "S2") (Fig. 5G).Lastly, LSN3160440-binding to the GLP1R con ned ECL1 and ECL2 from two conformational states ("S1" and "S2") to only the "S1" state, in which the extracellular mouth between ECL1 and ECL2 adopted the closed conformation to stabilize the peptide agonist (Figs.5H and 3H).

Selectivity Of GPCR Allosteric Modulators
Additional GaMD simulations were performed on arti cially generated computational models to examine binding selectivity of the MT7 and Cmpd-15 NAMs to the muscarinic and adrenergic receptors, respectively, as well as the LY2119620 and LY3154207 PAMs to the muscarinic and dopamine receptors, respectively.Flexibility changes were calculated by subtracting RMSFs of the cognate from the "noncognate" GPCRs of the modulators (Fig. 6).Furthermore, 2D free energy pro les of the heavy-atom RMSDs of orthosteric and allosteric ligands relative to their respective starting structures were calculated and shown in Supplementary Fig. 10.Overall, modulator binding in the "non-cognate" GPCRs resulted in higher complex uctuations and mostly larger conformational space compared to the cognate GPCRs, demonstrating the binding preference and selectivity of allosteric modulators towards their cognate subtypes.
Signi cantly higher uctuations were observed for binding of allosteric modulators to "non-cognate" GPCRs, especially in the allosteric pockets and various receptor domains, compared to their binding to the cognate GPCRs (Fig. 6).Compared to the MT7-bound M 1 R, the NAM showed moderately increased to much higher uctuations in the model M 2 R and M 4 R, respectively.Furthermore, NAM binding increased uctuations in the ECL2 of the M 2 R (Fig. 6A) and the TM4 extracellular end, ECL1, and ECL2 in the M 4 R (Fig. 6B).Compared to the Cmpd-15-bound β 2 AR, the NAM showed much higher uctuations in the "noncognate" subtypes of the α 1B AR, α 2A AR, and α 2C AR, and signi cantly increased uctuations of these three GPCR-antagonist complexes (Fig. 6C-6E).The exibility increase was smaller in the Cmpd-15-bound β 1 AR, likely due to the receptor similarity in its sequence and structure to the β 2 AR.Even so, binding of Cmpd-15 to the β 1 AR signi cantly increased uctuations in the TM2 extracellular end, ECL1, ECL2, TM4, ICL2, and H8 (Fig. 6F).Schober et al. uncovered that the binding preference of the LY2119620 PAM reduced from the M 2 R to the M 4 R and then M 1 R 75 .Here, LY2119620 binding to the M 1 R signi cantly increased uctuations in the TM2, TM3 and TM4 extracellular ends, ECL1, ECL2, and ECL3 compared to the LY2119620-bound M 2 R (Fig. 6G).In the LY2119620-bound M 4 R, PAM binding only slightly increased uctuations in the TM2 and TM3 extracellular ends, ECL1, and ECL2 compared to the M 2 R (Fig. 6H).Our simulation results were thus consistent with the experimental nding by Schober et al. 75 .Lastly, binding of the LY3154207 PAM to the D 2 R increased uctuations mostly in the TM1 and TM2 extracellular ends, TM4 intracellular end, ICL2, TM5, TM6, TM7, and the SKF-81297 agonist compared to the cognate D 1 R (Fig. 6I).

Conclusions
Allosteric modulators have emerged as more selective drug candidates than orthosteric agonist and antagonist ligands.However, many X-ray and cryo-EM structures of GPCRs resolved so far exhibit negligible differences upon binding of allosteric modulators.Consequently, mechanism of dynamic allosteric modulation in GPCRs remains unclear, despite their critical importance.In this work, we have integrated GaMD and DL in GLOW to map dynamic changes in free energy landscapes of GPCRs upon binding of allosteric modulators.By intersecting DL-predicted residue contacts with the highest gradient and residues with the largest exibility changes, characteristic residue contacts were selected for free energy pro ling to decipher the effects of allosteric modulator binding on GPCRs.The PAM and NAM binding primarily reduced uctuations of the GPCR complexes.NAMs stabilized the allosteric and antagonist-binding sites.PAMs stabilized the receptor extracellular domains, orthosteric agonist-binding pocket, and G protein coupling regions.Furthermore, the conformational space of GPCRs was signi cantly reduced upon modulator binding.The NAMs and PAMs con ned the GPCRs to mostly one speci c conformation for signaling.These effects transcended across class A and B GPCRs.
Furthermore, NAM and PAM binding were found selective towards their cognate receptor subtypes.
Signi cantly higher uctuations were observed for modulator binding to "non-cognate" GPCR subtypes, for which the orthosteric and allosteric ligands exhibited larger RMSDs.
While the effects of NAM and PAM binding on receptor dynamics were consistent across different GPCRs, exceptions were observed in the GPBAR bound by the INT777 PAM and the GLP1R bound by the PF-06372222 NAM.In the GPBAR, the fact that two molecules of the same charged ligand bound to the receptor could potentially create electrostatic repulsion, leading to increased uctuations in the orthosteric ligand and other parts of the receptor 45 (Fig. 3G).In addition, potential inaccuracies in especially the ligand force eld parameters could contribute to the inconsistencies observed in these two cases.Further ligand parameter optimization could be helpful to achieve more consistent results in these GPCR systems.In conclusion, we have deciphered the mechanism of dynamic allostery in class A and B GPCRs through DL of extensive GaMD simulations.Our ndings are expected to facilitate rational design of selective GPCR allosteric drugs.

Declarations
Author Contributions performed research, analyzed data, and wrote the manuscript.J.W. supervised GaMD simulations.Y.M. supervised the project, interpreted data, and wrote the manuscript.All authors contributed towards the nal version of the manuscript.GPCRs in the presence/absence of PAM were built (B).Three independent 500ns GaMD simulations were performed on each system (C).Residue contact maps were calculated for 150,000 x 44 GaMD simulation frames (D)and analyzed by Deep Learning, yielding saliency (attention) maps of residue contact gradients (E).Changes in root-mean-square uctuations (∆RMSFs) upon NAM/PAM binding in GPCRs were calculated from the GaMD simulations (F).If the absolute average ∆RMSF calculated from three simulations of a residue was smaller than the standard deviation of ∆RMSF, the exibility change for that residue was considered not signi cant and related residue pairs were neglected for further analysis.The characteristic residue contacts selected were those with ≥ 0.7 gradients and signi cant exibility changes upon modulator binding (G).They served as reaction coordinates for free energy pro ling of dynamic allostery of GPCRs (H).

Figure 1 Work
Figure 1