Structural Basis for Binding of Allosteric Drug Leads in the Adenosine A1 Receptor

Despite intense interest in designing positive allosteric modulators (PAMs) as selective drugs of the adenosine A1 receptor (A1AR), structural binding modes of the receptor PAMs remain unknown. Using the first X-ray structure of the A1AR, we have performed all-atom simulations using a robust Gaussian accelerated molecular dynamics (GaMD) technique to determine binding modes of the A1AR allosteric drug leads. Two prototypical PAMs, PD81723 and VCP171, were selected. Each PAM was initially placed at least 20 Å away from the receptor. Extensive GaMD simulations using the AMBER and NAMD simulation packages at different acceleration levels captured spontaneous binding of PAMs to the A1AR. The simulations allowed us to identify low-energy binding modes of the PAMs at an allosteric site formed by the receptor extracellular loop 2 (ECL2), which are highly consistent with mutagenesis experimental data. Furthermore, the PAMs stabilized agonist binding in the receptor. In the absence of PAMs at the ECL2 allosteric site, the agonist sampled a significantly larger conformational space and even dissociated from the A1AR alone. In summary, the GaMD simulations elucidated structural binding modes of the PAMs and provided important insights into allostery in the A1AR, which will greatly facilitate the receptor structure-based drug design.

(for neuropathic pain), but failed due to lack of efficacy 23,24 . Overall, these compounds still suffer from major limitations for pharmaceutical use, such as low solubility, affinity and cooperativity. To date, structural information regarding PAM interactions with the A 1 AR has been largely unavailable to guide previous drug design efforts aimed towards the development of therapeutically effective A 1 AR PAMs 5 . Mutagenesis and molecular modeling studies have suggested that the A 1 AR allosteric site may reside within the second extracellular loop (ECL2) 26,27 , however the precise location of the allosteric site and the molecular mechanisms underlying the allosteric modulation of PD81723 and other PAMs remain unclear. Recently, the first X-ray crystal structure of the A 1 AR (PDB: 5UEN) was determined by Christopoulos and coworkers 28 . In the structure, the A 1 AR was bound to an irreversible antagonist DU172, which forms a covalent bond with Tyr271 7. 36 in the transmembrane (TM) helix 7; superscript denoting Ballesteros-Weinstein residue numbering 29 . Compared with previous X-ray structures of the A 2A AR 30,31 , the A 1 AR exhibits a significantly wider extracellular cavity with a distinct conformation of the ECL2. Another similar X-ray structure was determined for the A 1 AR bound by the antagonist in the inactive state 32 . The X-ray structures serve as an excellent starting point for computational modeling and structure-based drug discovery of the A 1 AR.
Molecular dynamics (MD) is a powerful computational technique for simulating biomolecular dynamics on an atomistic level 33 . For GPCRs, MD has been applied to simulate binding of both orthosteric and allosteric ligands 28,34,35 . Using the specialized supercomputer Anton, Dror et al. performed microsecond-timescale MD simulations on the β 1 and β 2 adrenergic receptors (β 1 AR and β 2 AR) 34 . These simulations showed that antagonist and agonist ligands entered the receptor orthosteric site through an opening between ECL2-ECL3, which was suggested to be a dominant binding pathway of GPCR drugs. Subsequent Anton MD simulations captured the same orthosteric ligand binding pathway for the M 2 and M 3 muscarinic acetylcholine GPCRs (mAChRs) 28 . Binding of several known negative allosteric modulators (NAMs) to the extracellular vestibule of the M 2 receptor was observed in further Anton MD simulations 35 . The modulators formed cation-π interactions with aromatic residues in the receptor extracellular vestibule. The extracellular allosteric binding mode was confirmed by mutation experiments and later by the X-ray structure of the active M 2 receptor that is recognized by a PAM 36 . Despite these successes, direct MD simulations are computationally expensive for studying protein-ligand binding. They often suffer from insufficient sampling of slow ligand binders 28 and cannot capture ligand dissociation due to limited simulation timescales.
During the last several decades, many enhanced sampling methods have been developed to improve MD simulations [37][38][39][40][41][42][43] . Among these methods, metadynamics 44,45 , random acceleration MD 46,47 , steered MD 48 , temperature accelerated MD (TAMD) 28,49 , accelerated MD (aMD) 50,51 and Gaussian aMD (GaMD) [52][53][54] have been applied to simulate ligand binding to GPCRs 55 . Metadynamics was applied to simulate binding of a PAM to the δ-opioid receptor in the presence of an agonist 56 and calculate ligand binding free energies in the β 2 AR 57 . We performed aMD simulations on binding of the tiotropium antagonist, acetylcholine agonist and arecoline partial agonist to the M 3 muscarinic receptor 58 . In comparison with the previous Anton MD simulations 28 , aMD captured a similar ligand binding pathway, but with significant speedup (~80 times faster for agonist binding to the receptor orthosteric ligand-binding site). Using GaMD that provides unconstrained enhanced sampling and improved free energy calculations 52-54 , we also captured spontaneous binding of the agonist acetylcholine and identified its low-energy binding sites in the M 3 receptor 53 . The energetically preferred pathway of agonist binding identified from the GaMD simulation was similar to that found in previous long-timescale cMD 34 and aMD 58 simulations. Furthermore, we successfully applied GaMD to capture both dissociation and binding of the arecoline partial agonist in the M 2 receptor 59 . Therefore, GaMD is well suited for investigating ligand binding of large biomolecules such as GPCRs.
In this study, we have applied GaMD to simulate binding of allosteric drug leads to the A 1 AR. Extensive GaMD simulations using the AMBER and NAMD simulation packages at different acceleration levels captured spontaneous binding of two prototypical PAMs to the A 1 AR. The GaMD simulations also allowed free energy calculations to identify low-energy binding modes of the PAMs at the putative allosteric site formed by ECL2 of the receptor, which is highly consistent with the mutation experimental data. Furthermore, PAM binding was found to stabilize agonist binding at the receptor orthosteric site. Therefore, GaMD simulations have provided a greater understanding of the structural binding modes and allosteric effects of PAMs at the A 1 AR.

Results
GaMD simulations captured spontaneous binding of PAMs. Using the first X-ray crystal structure of the A 1 AR (PDB: 5UEN, Fig. 1A) 60 , we have performed all-atom GaMD simulations to investigate binding of two prototypical PAMs, PD81723 11,12 and VCP171 25 (Fig. 1B). The antagonist was removed from the X-ray structure and the agonist 5′-N-ethylcarboxamidoadenosine (NECA) placed in the receptor with atomic coordinates extracted from the A 2A AR X-ray structure (PDB: 2YDV), after aligning the two receptor transmembrane domains. GaMD simulations of the NECA-bound A 1 AR were performed in the absence and presence of two PAMs, PD81723 and VCP171. Each PAM was initially placed at least 20 Å away from the receptor (Fig. 1C). Multiple independent GaMD simulations were performed using AMBER and NAMD at different acceleration levels to investigate the PAM binding processes (Table 1).
With AMBER, GaMD simulations boosted both the total and dihedral energetic terms ("dual-boost GaMD") 52 on the NECA-bound A 1 AR, NECA-bound A 1 AR in the presence of PD81723 and NECA-bound A 1 AR in the presence of VCP171 provided boost potentials of 17.89 ± 5.23 kcal/mol, 18.36 ± 5.29 kcal/mol and 17.66 ± 5.23 kcal/mol, respectively. In comparison, dual-boost GaMD simulations using NAMD showed boost potentials of 11.77 ± 3.07 kcal/mol and 11.14 ± 3.07 kcal/mol for the NECA-bound A 1 AR in the absence and presence of PD81723, respectively. Further GaMD simulations were performed by boosting the dihedrals only ("dihedral GaMD") 53 using NAMD. The boost potentials were 6.04 ± 2.23 kcal/mol and 5.65 ± 2.14 kcal/mol for the NECA-bound A 1 AR in the absence and presence of the PD81723, respectively ( average and standard deviation of the boost potential (∆V avg and σ ∆V ) lead to higher acceleration in the biomolecular structural fluctuations. Therefore, the AMBER version of GaMD appeared to provide higher acceleration than the NAMD version, due to slightly different algorithms implemented for computing the potential statistics in the two packages (SI Method). Despite the different acceleration levels, all the GaMD simulations successfully captured spontaneous binding of PD81723 and VCP171 to the A 1 AR. Traces of the diffusing PAMs obtained in GaMD simulations with PD81723 and VCP171 are shown in Figs S1 and S2, respectively. Overall, PD81723 and VCP171 bound to a pocket formed by the ECL2 with the highest probability, although the PAMs could also transiently visit other regions of the A 1 AR. This agrees with previous mutagenesis experiments that alanine substitutions of residues in the ECL2 affected binding of the PAMs 26,27 . Therefore, the GaMD simulations successfully captured spontaneous binding of the A 1 AR PAMs PD81723 and VCP171.
GaMD predicted binding poses of PD81723 were consistent with structure-function data. Structural clustering of the PAMs and calculated free energies of the resulting structural clusters were performed on the GaMD simulation trajectories (see details in Methods). The lowest-energy binding poses of PD81723 at the ECL2 allosteric site obtained from dual-boost GaMD simulations using AMBER, dual-boost GaMD simulations using NAMD and dihedral GaMD simulations using NAMD are shown in Fig. 2A-C, respectively. Overall, PD81723 exhibited similar binding poses at the ECL2 allosteric site in the GaMD simulations at different acceleration levels. The trifluoro-phenyl group all pointed in the same direction towards Trp156, although the 2-amino-thiophene group was able to rotate slightly in the bottom part of the ECL2 pocket. The carbonyl oxygen always pointed towards the solvent, favoring hydrophilic interactions. The phenyl ring aligned in parallel with the short helix of ECL2, forming favorable hydrophobic interactions with protein residues Phe77, Val152, Ala155, Pro165 and Ile167. Common residues that were identified within 5 Å of the bound PD81723 in After removal of antagonist from the A 1 AR X-ray structure, NECA was placed in the orthosteric pocket with atomic coordinates copied from the A 2A AR X-ray structure (PDB: 2YDV) after aligning the two receptor transmembrane domains. Four molecules of each PAM were placed >20 Å away from the receptor. the lowest-energy poses from the GaMD simulations included Ala155, Ala159, Trp156, Pro165, Ile167, Phe77 and Lys173. Moreover, residues Asn148, Ala151, Val152, Glu153, Ser161, Val166, Lys168, Glu172 and Val174 appeared within 5 Å of PD81723 in one or two of the binding poses ( Fig. 2A-C). The lowest-energy binding poses of PD81723 obtained from the GaMD simulations were highly consistent with site-directed mutagenesis experiments 26,27 . Particularly, alanine substitution of Asn147, Asn148, Ser161, Pro165, Val166, Ile167, Glu170, Lys173 or Ile175 significantly enhanced the binding affinity (pK B ) of PD81723, while mutation of Glu172 to alanine had the opposite effect (Fig. 2D). Alanine substitutions of residues Asn147, Asn148, Gly160, Ser161, Lys173 and Glu170 decreased the binding cooperativity (logα I ) between PD81723 and NECA at the A 1 AR (Fig. 2E). Alanine substitutions of Asn148, Glu153, Ser161 and Ile167 decreased PAM efficacy logτ B(C) (Fig. 2F). Finally, mutation of Glu172 to alanine increased the functional cooperativity (logαβ) between PD81723 and NECA for A 1 AR-mediated inhibition of cAMP accumulation (Fig. 2G).
GaMD predicted binding mode of VCP171 was consistent with structure-function data. The lowest-energy binding mode of VCP171 identified from dual-boost GaMD simulations using AMBER is shown in Fig. 3A. VCP171 also recognized the putative allosteric site formed by the ECL2 and adopted a similar orientation to PD81723. The trifluoro-phenyl group was parallel to the short helix in ECL2. The thiol group pointed towards Asn148 and the additional phenyl group formed hydrophobic interactions with the Phe77 and Ile167 side chains (Fig. 3A). Residues found within 5 Å of the bound VCP171 included Phe77, Asn148, Ala151, Val152, Ala155, Ile167, Lys168, Glu172 and Lys173 (Fig. 3A). The majority of interacting residues were shared between the PAMs PD81723 and VCP171 (Fig. 2).
Residues found within 5 Å of the bound VCP171 were highly consistent with site-directed mutagenesis experiments 26,27 . Notably, mutation of Glu172 to alanine significantly decreased binding affinity of VCP171 (Fig. 3B). Alanine substitution of Trp146 and Trp156 decreased binding cooperativity between VCP171 and NECA at the A 1 AR (Fig. 3C). Alanine substitution of Asn148, Glu153, Arg154, Ser161, Ile167 and Ile175 decreased the receptor efficacy (Fig. 3D). Finally, alanine substitution of Arg154, Ser161, Trp156, Val174 and Ile175 decreased the functional cooperativity between VCP171 and NECA for A 1 AR-mediated inhibition of cAMP accumulation (Fig. 3E). Therefore, the binding mode of VCP171 obtained from the GaMD simulations was supported by experimental structure-function analysis 26,27 . assessed the influence of PAMs on NECA binding within the A 1 AR orthosteric site. Structural clusters of the agonist NECA were obtained in the absence and presence of PAM binding (Fig. 4). In the NECA-bound A 1 AR simulation system in the absence of an allosteric ligand, NECA sampled a large conformational space in the orthosteric pocket. Numerous structural clusters of NECA were identified from simulations of the NECA-bound A 1 AR system using dual-boost GaMD with AMBER (51 clusters, Fig. 4A), dual-boost GaMD with NAMD (4 clusters, Fig. 4B) and dihedral-boost GaMD with NAMD (6 clusters, Fig. 4C). In dual-boost GaMD simulations with AMBER, which provided the highest acceleration (Table 1), the NECA agonist was able to explore the entire orthosteric pocket (Fig. 4A). Whereas in the dual-boost and dihedral GaMD simulations using NAMD with lower acceleration levels, NECA sampled a smaller conformational space (Fig. 4B,C). Accordingly, the RMSD of NECA relative to the crystal conformation obtained from the 2YDO X-ray structure of the A 2A AR with two receptor TM domains aligned exhibited large variations during the AMBER dual-boost GaMD simulations (Fig. S3A), while the NAMD dual-boost and dihedral GaMD simulations showed smaller variations in the agonist RMSDs ( Fig. S3B and S3C). During GaMD simulations of the NECA-bound A 1 AR in the presence of PD81723, RMSD of NECA typically remained <5 Å with small variations (Fig. S4), except that it reached ~9 Å in one of the five dual-boost GaMD simulations using NAMD ("Sim5" in Fig. S4B). This suggested that PD81723 was able to stabilize NECA binding in the orthosteric site. We tracked PD81723 diffusion and identified structural clusters of NECA without and with PD81723 bound at the ECL2 allosteric site. Different numbers of structural clusters for the orthosteric agonist NECA with no PD81723 bound at the ECL2 allosteric site were obtained during the dual-boost GaMD simulations with AMBER (2 clusters), dual-boost GaMD simulations with NAMD (12 clusters) and dihedral  Table 1). In comparison, PD81723 binding to the ECL2 allosteric site led to fewer structural clusters and smaller conformational space of NECA in the orthosteric pocket ( Fig. 4 and Table 1). During the AMBER dual-boost GaMD, NAMD dual-boost GaMD and NAMD dihedral GaMD simulations, the number of structural clusters identified for NECA was 2, 9 and 1, respectively (Table 1). Upon binding of PD81723 at the ECL2 allosteric site, NECA sampled a smaller conformational space. Structural clusters of NECA identified in simulations of the "A 1 AR + NECA + PD81723" system with PD81723 bound at the ECL2 allosteric site using dual-boost GaMD with AMBER, dual-boost GaMD with NAMD and dihedral-boost GaMD with NAMD are shown in Fig. 4G-I. Movement of the NECA agonist was greatly reduced in the presence of PD81723. Therefore, the PD81723 PAM stabilized agonist binding at the orthosteric site of the A 1 AR.
Agonist dissociation was observed in the absence of PAM binding. Dual-boost GaMD simulations using AMBER were also performed on the system of NECA-bound A 1 AR in the presence of the VCP171 PAM.
While NECA showed small movements with mostly <10 Å RMSD relative to the 2YDV crystal conformation (Fig. S5), it escaped out of the receptor with >20 Å RMSD in one of the five GaMD simulations ("Sim3"). The structural clusters of NECA were identified from the GaMD simulations with and without VCP171 bound at the ECL2 allosteric site (Fig. 5). VCP171 binding to ECL2 greatly limited the conformational space of NECA agonist with only four structural clusters identified in the receptor orthosteric pocket ( Fig. 5A and Table 1). In the absence of VCP171 binding to the ECL2 allosteric site, NECA sampled a significantly larger conformational space and even dissociated from the A 1 AR (Fig. 5B). A large number of structural clusters (73 clusters) were identified for the diffusing NECA ( Table 1). The lowest energy clusters depict an agonist dissociation pathway, which connect the receptor orthosteric site, extracellular opening between the ECL2/ECL3, the allosteric site formed by the ECL2 and finally the solvent (Fig. 5B). This is consistent with previous simulation findings that ligand binding through the ECL2/ECL3 opening is an energetically preferred pathway of class A GPCRs 34,49,58,59 .
Furthermore, the open pocket formed by only the ECL2 serves as an additional metastable binding site of the agonist, as well as the target site of PAMs in the A 1 AR.

PAM binding promotes the formation of a salt bridge E172 ECL2 -K265 ECL3 in the A 1 AR extracellular vestibule.
We have examined protein residue interactions in the A 1 AR extracellular domains to understand the allosteric mechanism of stabilized agonist interactions within the orthosteric site in the presence of PAM binding to ECL2. Simulation analysis revealed a salt bridge between Glu172 ECL2 -Lys265 ECL3 in the extracellular mouth of the A 1 AR. In particular, the favorable hydrophobic interactions between the PAM and ECL2 within the allosteric pocket positioned Glu172 ECL2 to extend its side chain towards residue Lys265 ECL3 . The Glu172 ECL2 -Lys265 ECL3 salt bridge in the A 1 AR extracellular mouth was then closed upon PAM binding to the ECL2 allosteric site (Fig. 6). This predicted interaction is consistent with mutation experimental data that substitution of Glu172 to alanine significantly affected binding affinity of PD81723 and functional cooperativity between the PAM and agonist 26,27 (Fig. 2).
A 2D PMF profile of the Glu172 ECL2 -Lys265 ECL3 distance and RMSD of the NECA agonist relative to the starting bound conformation was calculated from AMBER dual-boost GaMD simulations of the NECA-bound A 1 AR system (Fig. 6A). Three low-energy conformational states, "open", "intermediate" and "closed", were identified from the free energy profile. Notably, this salt bridge adopted the open conformation in the X-ray structure of antagonist DU172-bound A 1 AR (PDB: 5UEN) 61 and intermediate conformation in the cryo-EM structure of adenosine-G i -bound A 1 AR (PDB: 6D9H) 62 , for which the Glu172 ECL2 -Lys265 ECL3 distance is 14.82 Å/15.14 Å (dimer in the 5UEN structure) and 7.12 Å, respectively. Upon PAM binding to the A 1 AR, the salt bridge changed to the closed conformation with 3.0 Å distance between the Glu172 ECL2 and Lys265 ECL3 (Fig. 6B). In addition, 2D PMF profiles of the Glu172 ECL2 -Lys265 ECL3 distance and the occupancy of PAMs at the ECL2 allosteric site were calculated from AMBER dual-boost GaMD simulations of the "A 1 AR + NECA + PD81723" and "A 1 AR + NECA + VCP171" systems ( Fig. 6C,D). Time courses of the PAM occupancy showed that PD81723 bound typically fast (within ~30 ns) to the ECL2 allosteric site in all the five GaMD simulations (Fig. S6A). Accordingly, the E172 ECL2 -K265 ECL3 distance decreased to ~3 Å (Fig. S7B) and the salt bridge stayed mostly closed in GaMD simulations of the "A 1 AR + NECA + PD81723" system (Fig. 6C). In comparison, VCP171 rarely bound to the ECL2 allosteric site during two of the five GaMD simulations (Sim3 and Sim4 as shown in Fig. S6B). The salt bridge sampled closed, intermediate and open conformations in the "A 1 AR + NECA + VCP171" system without VCP171 binding to the ECL2, but it was confined to the closed state upon binding of VCP171 to the ECL2 (Figs 6D and S7C). Therefore, PAM binding to the ECL2 allosteric site biased conformational ensemble of the Glu172 ECL2 -Lys265 ECL3 salt bridge towards the closed state, leading to stabilized agonist binding at the orthosteric site. This provided important insights into the mechanism of allosteric modulation in the A 1 AR.

Discussion
In this study, we have determined structural binding modes of prototypical PAMs in the A 1 AR through extensive GaMD enhanced simulations. The GaMD simulations have been performed using the AMBER and NAMD simulation packages at different acceleration levels. In the GaMD simulations, the A 1 AR PAMs bound to an allosteric site formed by ECL2, a finding that was highly consistent with experimental data 26,27 . Many of the A 1 AR residues identified within 5 Å of bound PAMs in the GaMD simulations, including ECL2 residues Asn148, Glu153, Ser161, Ile167 and Glu172, have previously been suggested to be important for PAM binding in structure-function studies. These studies demonstrated that alanine substitution of these residues significantly influenced PAM affinity, cooperativity and/or efficacy 26,27 . Therefore, these findings suggest the A 1 AR PAM allosteric site resides within the extracellular vestibule, predominantly involving interactions with ECL2. Furthermore, the GaMD simulations showed that the PAMs stabilized agonist binding at the A 1 AR orthosteric site. Specifically, compared to simulations performed in the absence of an allosteric ligand, agonist movement within the orthosteric site decreased in the presence of a PAM bound to the ECL2 allosteric site. This was correlated with conformational change in the Glu172 ECL2 -Lys265 ECL3 distance in the A 1 AR extracellular mouth. In GaMD simulations of the NECA-bound A 1 AR system, the Glu172 ECL2 -Lys265 ECL3 salt bridge sampled three low-energy conformational states ("Open", "Intermediate" and "Closed"). Importantly, the open and intermediate conformations have been determined in the X-ray structure of antagonist DU172-bound A 1 AR (PDB: 5UEN) 61 and the cryo-EM structure of adenosine-G i -bound A 1 AR (PDB: 6D9H) 62 , respectively. New high resolution A 1 AR structures co-bound with a PAM and agonist are required to confirm the predicted PAM-mediated stabilization of the Glu172 ECL2 -Lys265 ECL3 salt bridge in the closed conformation. Moreover, further GaMD simulations on PAM binding to the active A 1 AR using the recent cryo-EM structure of the adenosine-bound A 1 AR-Gi complex are subject to future study. Nevertheless, the present GaMD simulations have provided important insights into the mechanism of allosteric modulation at the A 1 AR.
In the absence of PAM binding to ECL2, the orthosteric agonist explores a significantly larger conformational space and could even dissociate from the A 1 AR through an opening between ECL2 and ECL3. This pathway connecting the orthosteric site and ECL2/ECL3 opening has also been suggested as an energetically preferred ligand binding pathway of other class A GPCRs, including the β 2 AR 34 , M 2 59 and M 3 mAChRs 53 . An orthosteric antagonist ZM241385 has also been observed to dissociate from the A 2A AR through a similar pathway in previous Anton simulations using TAMD 49 .
The extracellular allosteric site formed by only the ECL2 appears to be unique in the A 1 AR. Such an allosteric target site has not been identified in other GPCRs so far 63 . In the M 2 receptor, a PAM LY2119620 binds to the receptor extracellular vestibule formed by the TM2, TM6 and TM7 in addition to ECL2 as identified in X-ray crystallography 36 and MD simulations 35 . Sequence alignment of the four subtypes of ARs showed that while the seven TM helix bundle of the A 1 AR shares high similarity with the A 2A AR (71%), A 2B AR (70%) and A 3 AR (77%), the similarity is significantly reduced in the ECLs, being 43% for A 2A AR, 45% for A 2B AR and 35% for A 3 AR when compared with the A 1 AR. The entire ECL2 has low sequence conservation among the ARs (Fig. S8A).
Comparison of X-ray structures of the A 1 AR (PDB: 5UEN) 60 with the adenosine-bound A 2A AR (PDB: 2YDO) 64 also showed significant differences in the ECL2 conformations. The ECL2 forms a longer helix in the A 1 AR than in the A 2A AR and the helix adopts distinct orientations in the two receptors 60 . Two disulfide bonds, SSB1 and SSB2 that anchor ECL2 to ECL1 and Cys 3.22 in A 2A AR, respectively, are not conserved in the A 1 AR (Fig. S8B). This likely leads to higher flexibility of ECL2 in the A 1 AR and could play an important role in the binding of selective PAMs.
In summary, we have successfully identified a structural binding mode of A 1 AR PAMs through extensive all-atom GaMD simulations that is consistent with previous experimental structure-function analysis. The GaMD simulations provide important insights into the allosteric modulation mechanism of the A 1 AR, predicting that a salt bridge between Glu172 ECL2 -Lys265 ECL3 in the A 1 AR extracellular mouth is closed upon PAM binding to the ECL2, leading to stabilized agonist binding within the orthosteric site. The ECL2 appears to serve as the target site for A 1 AR PAMs, as well as an additional metastable binding site for orthosteric agonists. With remarkable divergence of residue sequences and conformations, the ECL2 presents an exciting target site for designing selective allosteric drugs of the A 1 AR. The GaMD simulations, together with mutagenesis data, will greatly facilitate future structure-based computer-aided drug design of novel A 1 AR PAMs.

Methods
Gaussian Accelerated Molecular Dynamics. Gaussian accelerated molecular dynamics (GaMD) is an enhanced sampling technique that works by adding a harmonic boost potential to reduce the system energy barriers 52 . GaMD accelerates biomolecular simulations by orders of magnitude. GaMD does not require predefined collective variables. Compared with the enhanced sampling methods that rely on careful selection of the collective variables, GaMD is of particular advantage for studying "free" protein-ligand binding processes 37,52 . Moreover, because the boost potential follows a Gaussian distribution, biomolecular free energy profiles can be properly recovered through cumulant expansion to the second order 52 . GaMD thus solves the energetic reweighting problem as encountered in the previous aMD method 50,65 for free energy calculations. GaMD has been implemented in the widely used AMBER 52,66 and NAMD 53 packages. It has allowed us to characterize protein folding, protein-ligand binding, protein-protein binding and protein-nucleic acid interactions 52,53,59,67,68 . Details of the method have been described in previous studies 52,53 . A brief summary is provided here.
Consider a system with N atoms at positions r  = {  r 1 …,  r N }. When the system potential V( r  ) is lower than a reference energy E, the modified potential V * (  r ) of the system is calculated as: 2 where k is the harmonic force constant. The two adjustable parameters E and k are automatically determined based on three enhanced sampling principles 52 . The reference energy needs to be set in the following range: where V max and V min are the system minimum and maximum potential energies. To ensure that Eqn. (2) is valid, k must satisfy: , then 0 < k 0 ≤ 1. The standard deviation of ΔV needs to be small enough (i.e., narrow distribution) to ensure proper energetic reweighting 69 where V avg and σ V are the average and standard deviation of the system potential energies, σ ΔV is the standard deviation of ΔV with σ 0 as a user-specified upper limit (e.g., 10 k B T) for proper reweighting. When E is set to the lower bound E = V max , k 0 can be calculated as: Alternatively, when the threshold energy E is set to its upper bound = + E V min k 1 , k 0 is set to: if ″ k 0 is found to be between 0 and 1. Otherwise, k 0 is calculated using Eqn. (3). Similar to aMD, GaMD provides options to add only the total potential boost ΔV P , only dihedral potential boost ΔV D , or the dual potential boost (both ΔV P and ΔV D ). The dual-boost simulation generally provides higher acceleration than the other two types of simulations for enhanced sampling 51 . The simulation parameters comprise of the threshold energy values and the effective harmonic force constants, k 0P and k 0D for the total and dihedral potential boost, respectively.
For energetic reweighting of GaMD simulations to calculate potential of mean force (PMF), the probability distribution along a reaction coordinate is written as p * (A). Given the boost potential  ΔV r ( ) of each frame, p * (A) can be reweighted to recover the canonical ensemble distribution, p(A), as: where M is the number of bins, β = k B T and  βΔ ⟨ ⟩ e V r j ( ) is the ensemble-averaged Boltzmann factor of ΔV r ( )  for simulation frames found in the j th bin. The ensemble-averaged reweighting factor can be approximated using cumulant expansion: where the first two cumulants are given by: The boost potential obtained from GaMD simulations usually follows near-Gaussian distribution 68 . Cumulant expansion to the second order thus provides a good approximation for computing the reweighting factor 52,69 . The reweighted free energy F(A) = −k B T ln p(A) is calculated as: System Setup. The first X-ray crystal structure of the A 1 AR (PDB: 5UEN) 60 was used to set up the simulation system. After removal of antagonist, the NECA agonist was placed in the orthosteric pocket with atomic coordinates extracted from the A 2A AR X-ray structure (PDB: 2YDV) after aligning the two receptor transmembrane domains. Four molecules of each PAM (PD81723 and VCP171) were initially placed >20 Å away from the receptor.
Two systems "A 1 AR + NECA + PD81723" and "A 1 AR + NECA + VCP171" were prepared for the simulations. In addition, the system NECA-bound A 1 AR in the absence of PAMs was also added for comparison (Table 1) The protein that was fused into the receptor to replace intracellular loop 3 (ICL3) for crystallizing the receptor structure was omitted. All chain termini were capped with neutral groups (acetyl and methylamide). The disulphide bonds that were resolved in the crystal structure were maintained in the simulations. Using the psfgen plugin in VMD 70 , protein residues were set to the standard CHARMM protonation states at neutral pH. Then the receptor was inserted into a palmitoyl-oleoyl-phosphatidyl-choline (POPC) bilayer with all overlapping lipid molecules removed using the Membrane plugin in VMD 70 . The system charges were then neutralized at 0.15 M NaCl using the Solvate plugin in VMD 70  152 lipid molecules, ~15,600 water molecules and a total of ~72,300 atoms. Periodic boundary conditions were applied on the simulation systems.
Simulation Protocols. The CHARMM36 parameter set 71 was used for the protein and POPC lipids. For agonist NECA and PAMs PD81723 and VCP171, the force field parameters were obtained from the CHARMM ParamChem web server 27,72 . Initial energy minimization and thermalization of the A 1 AR system follow the same protocol as used in the previous GPCR simulations 58 . In the present GaMD simulation, the threshold energy E for adding boost potential is set to the lower bound 52,53 . The simulations included 2 ns cMD, 50 ns equilibration after adding the boost potential and then multiple independent production runs lasting 200-500 ns with randomized atomic velocities. The GaMD simulations are summarized in Table 1.
For the "A 1 AR + NECA" and "A 1 AR + NECA + PD81723" systems, GaMD simulations were performed using AMBER at the dual-boost level 52 , and NAMD at the dual-boost and dihedral acceleration levels 53 . For the "A 1 AR + NECA + VCP171" system, dual-boost GaMD simulations using AMBER were performed. The GaMD simulations were carried out using AMBER 16 52,73 and/or NAMD2.13 53,74 . GaMD production frames were saved every 0.1 ps for analysis.
Simulation Analysis. The VMD 70 and CPPTRAJ 75 tools were used for trajectory analysis. The Density Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm 76 was applied to cluster the diffusing ligand molecules for identifying their highly populated binding conformations. The frames were sieved at a stride of 200 for clustering. The remaining frames were assigned to the closest cluster afterwards. The distance cutoff for DBSCAN clustering was set to 4 Å for the PAMs and 1.5 Å for the NECA agonist. Finally, the PyReweighting toolkit 69 was applied to compute free energy values of the ligand structural clusters.