Communication over the Network of Binary Switches Regulates the Activation of A2A Adenosine Receptor

Dynamics and functions of G-protein coupled receptors (GPCRs) are accurately regulated by the type of ligands that bind to the orthosteric or allosteric binding sites. To glean the structural and dynamical origin of ligand-dependent modulation of GPCR activity, we performed total ~ 5 μsec molecular dynamics simulations of A2A adenosine receptor (A2AAR) in its apo, antagonist-bound, and agonist-bound forms in an explicit water and membrane environment, and examined the corresponding dynamics and correlation between the 10 key structural motifs that serve as the allosteric hotspots in intramolecular signaling network. We dubbed these 10 structural motifs “binary switches” as they display molecular interactions that switch between two distinct states. By projecting the receptor dynamics on these binary switches that yield 210 microstates, we show that (i) the receptors in apo, antagonist-bound, and agonist-bound states explore vastly different conformational space; (ii) among the three receptor states the apo state explores the broadest range of microstates; (iii) in the presence of the agonist, the active conformation is maintained through coherent couplings among the binary switches; and (iv) to be most specific, our analysis shows that W246, located deep inside the binding cleft, can serve as both an agonist sensor and actuator of ensuing intramolecular signaling for the receptor activation. Finally, our analysis of multiple trajectories generated by inserting an agonist to the apo state underscores that the transition of the receptor from inactive to active form requires the disruption of ionic-lock in the DRY motif.


Author Summary
As the key signal transmitters of a number of physiological processes, G-protein coupled receptors (GPCRs) are arguably one of the most important therapeutic targets. Orchestration of the intra-molecular signaling across transmembrane domain is key for the function of GPCRs. To investigate the microscopic underpinnings of intramolecular signaling that regulates the activation of GPCRs, we performed molecular dynamics simulations of the Introduction G-protein coupled receptors (GPCRs) are one of the most versatile membrane proteins that mediate cellular responses to a myriad of extracellular signals associated with our perception, cardiovasicular, and immune functions [1]. GPCRs relay extracellular signals to the cytoplasmic domain. For a given extracellular signal, there is a corresponding GPCR subtype that processes the incoming signal to the cellular downstream [2]. Consisting of seven transmembrane (TM) helices, each of which is connected to the next TM helix by either an intracellular loop (ICL) or an extracellular loop (ECL), the interior of GPCR forms an interhelical residueto-residue interaction network that can transmit the signal specific to the ligand and/or receptor subtype.
Binding of an agonist to the orthosteric site leads to conformational rearrangement of TM helices, transforming the inactive conformation to the active one, which in turn enables accommodation of G-proteins and intracellular signal transductions [1,3,4]. It is widely appreciated that together with highly conserved residues in the TM region, the activation of the receptors belonging to the class A GPCR family is regulated by a set of fingerprint residues called "microswitches" [1,5], which commonly include three structural motifs: DRY, CWxP, and NPxxY motifs ("x" denotes any amino acid residue). Upon binding of an agonist, the microswitch residues change the orientation of their side chain, which transforms the global configuration of TM helices into the active form and helps the intracellular domain accommodate G-proteins [5]. In contast, binding of an inverse-agonist suppresses the GPCR function below the basal level, to which the apo or antagonist-bound forms of GPCRs is likely to be tuned [1].
The allostery, a long-range communication between two remote sites [6][7][8][9], is one of the key determinants for functions in many biomolecules. And it is of particlular interest for GPCRs because most of dynamic processes associated with GPCR activation and suppression involve signaling between extracellular and intracellular domains, established across the TM domain. Although the term "allosteric modulation" is often strictly distinguished from the orthosteric signaling in GPCR community, both can be considered allosteric signaling in that their signalings are both physically long-ranged. While high-resolution crystal structures in the agonist-bound or antagonist-bound states provide unprecedented view of ligand-dependent modulation of GPCR conformation, minor conformational difference between distinct receptor states make it difficult to glean the microscopic origin of intramolecular signaling and modulation of GPCR activity. Unlike molecular machines [10][11][12][13] and enzymes that exhibit large conformational changes [14,15], to which one can apply a method such as normal mode analysis and its variation [16,17], it is not straightforward to study the allosteric dynamics of biomolecules like GPCRs [18][19][20][21] whose structural changes between the inactive and active forms are relatively small (RMSD between inactive and active forms are only 1.78 Å for A 2A adenosine receptor; and 2.96 Å for b 2 adrenergic receptor whose active state structure is crystalized with G s -protein). It is noteworthy that there is a growing realization that a certain class of proteins can display allosteric responses by modulating conformational fluctuations, but with little conformational change [22][23][24][25]. Global conformational changes themselves are not the sole physical origin of protein allostery. Thus, to gain a better understanding of allostery or longrange intramolecular signaling of GPCR, it is imperative to probe the dynamical features of key structural components and their correlations. Although experimental methods, such as fluorescence resonance energy transfer (FRET), and site-directed spin labeling (SDSL), nuclear magnetic resonance (NMR) [26][27][28][29][30], are useful to monitor the conformational changes of specific key residues, it is generally difficult to simultaneously probe the multiple sites of molecule and elucidate dynamic origin leading to allosteric molecular responses. Much effort has recently been made to study the mechanism of GPCR activation via molecular dynamics (MD) simulations [31][32][33][34][35][36] by exploring the dynamic characteristics of GPCR conformations; however, systematic studies of cross-correlations among key structural motifs including microswtiches have not been carried out.
From sequence analysis, structural and biochemical studies using mutagenesis, there is a general consensus for the class A GPCRs that 18 hotspot residues, called microswitches, are responsible for the activation mechanism [1,3]. Recently we have shown that when GPCR structures are represented by a network of inter-residue contacts, many of those microswitches [1,5] retain high betweenness centrality values [21]. According to the network theory, vertices with high betweenness centralities in a given network are the sites that mediate flow of signal over the network [37,38]. A removal or alteration of such vertices from the network, which can be realized in the form of deletion [39] or more practically mutation into glycine [21], could impair the network communication and be deleterious to intramolecular signaling. Among the 18 microswitches, our network analysis in Ref. [21] could identify 11 of them in A 2A adenosine receptor (A 2A AR), and also showed that the majority of intra-molecular signaling pathways connecting the highly correlated residues between extra and intra-cellular domains pass through the 11 microswitches. However, our previous study provided only a static view of intramolecular allosteric wiring. Thus, to better understand the intramolecular signaling of A 2A AR, it would be of great interest to directly probe the dynamics associated with these microswitches and their cross-correlations.
In this study, we performed each of 1 msec all-atom molecular dynamics simulation of A 2A AR in apo, antagonist-bound (ZM-241385), [2,40] and agonist-bound (UK-432097) [41] forms (see Methods); and monitored the local dynamics of the microswitches. Based on the dynamics around the 11 microswitches, each of which displays molecular interactions that switches between two distinct states, we defined 10 ON/OFF binary switches as the microscopic components to describe the interhelical dynamics (note that the two microswitch residues of ionic-lock compose one on/off switch). Describing the dynamics of the receptor by employing 10 microscopic components amounts to projecting the entire dynamics into 10 local variables and considering 2 10 microstates. By projecting the receptor dynamics on 2 10 microstates, we show how the dynamics of each ligand-dependent macrostate differs from each other. Our simulations and analysis show that compared to that of the apo and antagonist-bound forms, the 10 binary switches of agonist-bound form retain greater crosscorrelation and coherence, which is consistent with the notion that GPCR in agonist-bound form has a functional fidelity to transmitting its intra-molecular signals across the TM domain [21]. Detailed knowledge on the local dynamics and their dynamic correlation elucidated in this study for A 2A adenosine receptors could be useful for the discovery of effective drugs.

Results
Among the 18 microswitches of A 2A AR, our previous study calculating the betweenness centrality of residue interaction network suggested that 11 of them are the putatively important spots for the intramolecular signal transmission [21], which include N24 1 [42]. Examining the local dynamics around the eleven allosteric hotspot residues [21] at different receptor states, we have defined 10 binary (on/off) switches. Below we first explain how we have defined each switch to be binary, which maximally separates the dynamic features of antagonist-bound state from those of agonist-bound state.
Ten binary switches defined from eleven hotspot residues for intra-molecular signaling.
Microswitches at the interfaces between TM1, TM2, and TM7 (Switches 1 & 2): In the class A GPCR family N24 and D52 are the most conserved microswitch residues in TM1 and TM2, respectively, with high betweenness centrality values [21]. Comparison of the H-bond network between three different receptor states reveals notable differences in the receptor configuration around them ( Fig. 1 and Fig. S1). First, in the apo form, H-bonds of N24 and D52 with S281 7.46 contribute to the stable conformation of TM helices. Binding of antagonist further stabilizes this conformation by incorporating an additional H-bond between N24 and D52. By contrast, binding of an agonist leads to the disruption of the H-bond between D52 and S281 (HB D52−S281 ) and forms a new H-bond between S281 and N280 7.45 , which is responsible for helix bending in TM7 (see below). Ligand-dependent switching dynamics of the H-bond of S281 from one side to the other is reminiscent of the salt-bridge switching, which is often observed at the interface between subunits in multisubunit proteins [10,43]. Hence, we designate HB N24-D52 and HB D52-S281 as the switch 1 (S1 ) and the switch 2 (S2 ), respectively.

DRY motif and ionic-lock (Switches 3 & 4):
The orientation of the TM5-ICL3-TM6 relative to the ICL2 is the hallmark of GPCR activation as the 10 o tilt of TM5-ICL3-TM6 enables the G-protein accommodation. Among others, the "ionic-lock", the salt-bridge between R102 and E228, which keeps the TM3 and TM6 in close proximity, is considered the key structural motif that directly regulates the orientation of TM6 helix [44]. Our simulations find that, in the antagonist-bound state, the ionic-lock is intact (d R102 Cz -E228C d ≈ 4 Å) ( Fig. 2a and the top panel of Fig. 2b); while it is rarely formed in the agonist-bound state (d R102-E228 ≈ 10 Å). In the apo state, the ionic-lock maintains the identical distance with that of the antagonist-bound state, but occasionally breaks (d R102 Cz -E228C d ≈ 7 Å) and reforms (d R102 Cz -E228C d ≈ 4 Å) ( Fig. 2 and Fig. S2).
The ionic-lock and the inter-residue distance (d L110-A221 ) between L110 C a in ICL2 and A221 C a in ICL3 were simultaneously probed in Fig. 2b. It highlights the correlation between the ionic-lock and the orientation of TM5-ICL3-TM6 relative to TM3. Dependence of d L110-A221 on the receptor state is clearly observed; d L110-A221 ≈ 22 Å (agonist-bound), 7 Å (antagonist-bound), and 14 Å (apo), which quantifies the position of TM6 helix relative to TM3 helix. For the apo form, the distance between ICL2 and ICL3 varies concomitantly with the state of ionic-lock (see the time traces after 500 nsec of Fig. 2b). The scatter plot of d R102-E228 and d L110-A221 (right panel of Fig. 2b) also shows a clear correlation between the dynamics of ionic-lock and the orientation of TM5-ICL3-TM6 domain. The status of ionic-lock affects the H-bond network around the DRY motif ( Fig. 2c and Fig. S3). In the agonist-bound form, E228 forms the H-bond with R107 instead of R102. Disruption of ionic-lock leads to release of TM6 from TM3, but as a counterbalance a new H-bond (HB T41-D101 ) binds TM3 tightly with TM2. As HB T41-D101 and the ionic-lock in the DRY motif can be used to discern the ligand-dependent macrostate of GPCR, we define HB T41-D101 and the ionic-lock as S3 and S4 , respectively.
Rotameric microswitches (Switches 5, 6, 7, 8, & 10): Ligand-dependent microstates of the five residues W129, Y197, W246, N284, and Y288 are best characterized by the rotameric angles of their side chains (Fig. S4), which allow us to define these residues as S5 , S6 , S7 , S8 , and S1 , respectively: (i) In the agonist-bound state, W129 (S5 ) shows asymmetric bimodal distribution of dihedral angle P(w 2 ) that has a major peak at ≈ −20 o and a minor peak at ≈ 90 o . In the apo and antagonist-bound forms, P(w 2 ) is unimodal with a peak at 90 o ( Fig [5], exhibits receptor state-dependent rotamer angle dynamics between ≈ +90 o and ≈ −90 o (Fig. 3c, Fig. S5c, and Supporting Movie M1); (iv) N284 (S8 ) and Y288 (S10 ) in the NPxxY motif show receptor-state dependent dihedral angle distribution ( Fig. 3d and Fig. S5d). Meanwhile, the dihedral angle distributions of P285, which also belongs to NPxxY motif, are little affected by the receptor state. Thus, we use the kink angle of P285 instead of its dihedral angle for defining the switch (see below).
Helix bending induced by proline kink (Switch 9): The geometry of proline residue constrains the backbone conformation and cause a sharp kink in alpha helix structures. Calculation of the helix bending angle [45] shows that in TM7 helix, the greatest kink is formed around P285, in particular, in the agonist-bound form (Fig. 3e). Comparison between structures for different receptor states suggests that the H-bonding between N280 7.45 and S281 7.46 , in the agonist-bound state, contributes to the helix bending in TM7. The scatter plot in Fig. 3e indicates that there is a direct correlation between TM7-helix bending (y TM7 ) and HB N280-S281 .
Projection of GPCR dynamics onto 10-binary switches. The dynamic features of the 11 microswitch residues that display two-state switch-like molecular interactions allow us to define 10 binary switches (Fig. 4). We choose the value of the i-th switch (s i ) as 1 or 0 (ON or OFF) in reference to the switch state in the agonist-bound form, so that the "similarity" of each switch can be assessed in reference to the agonist-bound form. For example, the HB D52-S281 , corresponding to S2 , is disrupted in the agonist-bound form. In this case we set s 2 = 1 for disrupted H-bond and s 2 = 0 for intact H-bond. Also, for S7 whose configuration is best described using are colored in black, blue, and red for the apo, antagonist-bound, and agonist-bound states, respectively, and their histograms are shown on the right side of the plots. On the right panel, scattered plot using the distances of R102-E228 and L110-A221 is shown with the yellow arrow depicting the conformational transition between three ligand forms. (c) Diagram of the H-bond network in the three receptor states. rotameric angle, we consider S7 to be in the ON state if the dihedral angle (C a -C b -C g -C d1 ) of W246 (see Fig. S4) is equal or less than −50 o ; otherwise, it is in the OFF state (s 7 = 0). Our 10 switch representation resembles the strategy in studying protein folding problem using "correctness" of the configuration of each residue with respect to the native state [46]. . S5 to S10 defined from the rotameric switches in TM4, TM5, TM6, and TM7. Rotameric states of (a) W129, (b) Y197, (c) CWxP motif, and (d) NPxxY motif are compared for the apo (white), antagonist-bound (cyan), and agonist-bound (pink) forms. (e) Helix bending in TM7. The helix bending angle (bendix) of TM7 was calculated using bendix program [45]. The helix is displayed as a cylinder marked with the heatmap ranging from 0 o to 20 o . The scatter plot on right side depicts the relationship between H-bond of N280-S281 and the bending angle of the TM7 helix (apo: black, antagonist-bound: blue, agonistbound states: red). The average bending angles are annotated with the symbols, X. Representing GPCR conformation in terms of the 10 binary switches amounts to "choosing" multiple progress coordinates (or multi-dimensional order parameters) to probe the allosteric dynamics of GPCR from the inactive to active state. The assumption that 10 binary switch can faithfully represent the dynamics of GPCRs leads to in total, 2 10 possible microstates; each microstate is expressed using binary number from 0000000000 (2) to 1111111111 (2) with each digit denoting the switch number from 1 to 10. These binary numbers can also be expressed with a decimal number from 0 and 1023 (Fig. S6). The time traces projected on the 10 binary switches and 2 10 microstates are shown in Fig. 5a and Fig. 5c, respectively. The average value of each switch, calculated in different receptor state x as 0 hs x,i i 1 (Fig. 5b), where hs x,i i with i denoting the switch index is the value of switch averaged over the simulation time, indicates that on average switches are ON in the agonist-bound form, OFF in the antagonist-bound form, and they lie in between in the apo form. The difference among the three receptor forms becomes more evident in terms of the population of microstates (Fig. 5d). The statistics of microstates shows that the receptor occupies different population of microstates depending on the type of ligand (Fig. 5). Among the entire microstates as summarized in  the switches in the agonist-bound form are in 1023th state (1111111111 (2) ), and ≈ 31.28% are in 895th state (1101111111 (2) ). (iii) Lastly, in the apo form of GPCR, on an average, S1 , S3 , S10 are ON state, while S2 , S4 , S6 , S8 , S9 are OFF state. Microstates that constitute the major population of the apo form are 1000000101 (2) (31.32%) and 1001000101 (2) (10.56%).
To quantify the statistical similarity explored by two different receptor states, say a and b (a ≠ b) in the 10-switch representation, we employ Hamming distance: where jxj is the absolute value (or modulus) of x, and hs x,i i is the average value of i-th switch  (i = 1,2. . .,10) in the receptor form x, which is calculated in Fig. 5b. Since 0 hs x,i i 1, it is expected that 0 d ab 10. The more similar, the smaller the value of d ab should be. We obtain d apo-ago = 4.74, d apo-antago = 2.68, d ago-antago = 7.31. Next, the "complexity" of each macrostate, defined with the ensemble of microstates {i = 1. . .N s } where N s = 2 10 , is quantified using the Shannon entropy: where p i is the probability of occupying the i-th microstate as in Fig. 6. When the receptor explores only a single state (p k = 1 and p i ≠ k = 0), the value of I should be I = I min = 0; and if all the 2 10 states are uniformly explored (p k = 2 −10 for all k) then I = I max = 10. Thus, the larger the value of I, the more diverse microstates are explored, which is also gleaned from Fig. 6. We obtain I antago ≈ 2.52, I apo ≈ 4.22, I ago ≈ 3.70, for antagonist-bound, apo, and agonist-bound state, resepctively; and hence the apo state explores the most diverse configurational space. Fig. 5e shows a schematic of relationship between the three receptor states combining the analyses using Eq. 1 and Eq. 2. The apo form is more similar to the antagonist-bound form than to the agonist-bound form, which is consistent with the general notion that agonist contributes to the active GPCRs while both apo and antagonist contribute to the inactive state.
Cross-correlations of the dynamics between binary switches. To identify the correlation between the ON/OFF dynamics of binary switches, we calculated their cross-correlation (C ij ) by using the conformational ensemble from the simulations.
where ds i = s i −hs i i is the variation of the switch value from its mean. C ij assesses the extent of coherence in the "change" in switch dynamics between the i-th and j-th switches. Marked differences of the correlation pattern are observed in the three distinct receptor states (Fig. 7a): (i) The apo state (the middle panel in Fig. 7) has only one positive correlation (C S4S8 > 0: 25 ), and many other negative correlations (C S3S4 , C S3S7 , C S3S8 , C S1S9 < À0:25 ); (ii) By contrast, in the antagonist-bound state (the left panel in Fig. 7), a positive correlation ( > 0.25) is detected only between the S3 and S5 , and the negative correlations present in the apo state are suppressed; (iii) The agonist-bound state has a greater number of the positive correlations between the switches. The diagrams in Fig. 7b illustrate how the allosteric couplings are established among the switches, especially highlighting many positive couplings among switches in the agonist-bound state. Most notably, S7 (W246), a central rotameric switch located at the deep bottom of ligand binding cleft, displays direct couplings with 6 other switches S1 , S2 , S3 , S5 , S8 , and S9 , and additional positive correlations are observed in C S1S2 , C S2S5 , C S3S8 , C S5S8 . This suggests that intramolecular signaling over the entire structure can be initiated by stimulating the S7 . As can be confirmed from simulations, the indole 6-ring of W246 is within the range of hydrophobic interaction (4-5 Å) [47] with the ethyl group of the agonist (UK-432097), while such direct interaction with W246 is missing in the antagonist (ZM-241385) (see Fig. 7c). Furthermore, there is a marked difference between the antagonist and agonist configurations in the orthosteric binding site; the antagonist is not stably poised as the agonist in the binding cleft (Fig. 7c, middle panel and Fig. S7). The dispersion of the center of position for each ligand is s antago = 2.204 Å and s ago = 0.854 Å (Fig. S7). Although the residues at the binding site adopt diverse configurations, some key residues retain persistent interaction with the receptor (Fig. S8a). Both the antagonist and agonist generally maintain the H-bondings with N253 and E169 residues (Fig. S8b). Also, the adenine rings of the ligands interact with the phenyl ring of F168 via p-p interaction. Fig. S8c shows that the adenyl ring of the agonist interact directly with F168 while that of antagonist does not. Occasionally, the antagonist spins or flips at the binding site, indicating the lack of interaction with F168 residue. The polar residues located at the bottom of the binding pocket, i.e., T88, S277, and H278, only interact with the agonist. Especially, T88, one of the hot spot residues in our earlier work [21], maintains the stable H-bonding with the agonist. Along with stable configuration of the agonist ligand (UK-432097) in the binding cleft (See Fig. S7), the widespread cross-correlation among the switch dynamics of S1 , S2 , S3 , S5 , S7 , S8 , and S9 in the agonist form suggests that the interaction between agonist and W246 actuates allosteric signaling and contributes to a stable activation of the receptor function.
Effect of inserting agonist to the apo state. Anticipating a detectable conformational change from inactive to active state, we generated 4 additional time traces by inserting agonist to the orthosteric binding site in the simulation trajectories of the apo state (Fig. 8a, Fig. S9. See Methods for the details of the simulation). The first two traces (cases 1 & 2) were generated by inserting agonist at 125 ns and 150 ns when the ionic-lock was still intact (s 4 = 0) and were simulated for ≈ 750 ns. The second two traces (cases 3 & 4) were generated at 595 ns and 625 ns when the ionic-lock was disrupted (s 4 = 1), and were simulated for ≈ 250 ns. The consequences of the insertion of agonist, summarized in Fig. 8, is still minor, which is evident when the average value of each switch hs i i (i = 1,. . . ,10) is compared (Fig. 5b middle panel versus Fig. 8b). This is not so surprising given that time scale of our simulation ( < 1 msec) is still too short to see a complete transition from inactive to active state in GPCRs, which is typically longer than milisecond time scale [48]. Although the overall trend of hs i i looks similar, each trace from the insertion of agonist explores distinct microstate population (Fig. 8c). Compared with the cases 1 & 2, the 10 switch states were closer to those of agonist-state in the cases 3 & 4 when the agonist was inserted to the receptor with disrupted ionic-lock (s 4 = 1); in particular, the value of s 4 for the cases 3 & 4 is greater and more variable (Fig. 8a,b) although the case 4 shows the stable rebinding of ionic-lock after 120 ns (s 4 = 1 ! 0, the rightmost panel in Fig. 8a). The complexity (Eq.2) of microstate ensemble explored in each simulation, are given in the table of Fig. 8d and the Hamming distances of the cases 1-4 from the three macrostates are calculated in Fig. 8e. The cases 3 and 4 are closer to the agoinist-bound state with d ago-case3 = 4.11 and d ago-case4 = 4.21 (Fig. 8d, e) than the cases 1 and 2 (d ago-case1 = 5.92, d ago-case1 = 6.60). It is noteworthy that the binding of agonist to a receptor with stable ionic-lock (cases 1 & 2) has driven the ensemble of microstates away from agonist-bound form and bring the ensemble close to the antagonist-bound state. Disruption of the ionic-lock in DRY motif is required for the activation of A 2A AR.

Discussion
Given that typical time scale associated with the transition from inactive to active GPCR conformers is ≳Oð1Þ msec [49], it is practically impossible to capture the complete evidence of transition from inactive to active state using a single time trace lasting only 1 msec [50]. Although the recent improvement in computational power and various computational strategies have ameliorated this time scale gap and have played significant roles in elucidating new facets of GPCR dynamics [31][32][33], the gap of time scale between all-atom molecular dynamics simulation and experiments is still a serious problem in linking the computation of biomolecular dynamics with experimental observation [51]. Rather short in time scale ( * 1 msec), the trajectories from our simulations, which essentially probe the dynamics of receptor in terms of the noisy local variables, enable us to map the intramolecular signaling of three different ligand states. When it comes to the sampling efficiency of each macrostate, 1 msec time scale is not too short to observe the individual dynamics of side chain flips that generally occur in picoseconds to nanoseconds timescale [52], and we tried to decipher collective signals out of those noisy individual signals. The dynamic characteristics of the three macrostates are well discerned in consistent with the general notion of GPCRs, and our analysis based on the dynamics of 10 binary switches quantifies the differences among the thre statese.
Our simulation results provide a comprehensive view of ligand-dependent conformational dynamics and show how the remote allosteric hotspots of A 2A AR are coupled each other, which lend support on the pre-existing experimental data. For example, W246 has long been proposed to function as a central rotamer switch in the activation of GPCRs [5]. While the rotameric change of W246 was confirmed in spectroscopy [53] and some MD simulation studies [5,34], X-ray crystal structures of A 2A AR in both antagonist-bound and agonist-bound states show little difference of the position and dihedral angles of this residue. Our simulation captures the transition of the indole ring of W246 in the presence of agonist (Supporting Movie M1). In addition, our simulations clearly demonstrate characteristics of local dynamics of microswitches such as ligand-dependent state of ionic-lock, and the ligand-dependent variation of the rotameric angles in Y197 and Y288 [2,4,41,54].
Until recently, there are some other studies that have compared the conformational differences of A 2A AR depending on the receptor states [34][35][36]. Pang et al. studied dynamic behaviors of A 2A AR induced by the binding of two distinct antagonists [35]. Lee et al. simulated the ligand-dependent cholesterol interactions in A 2A AR [36]. Exploring the distinct structural states that resemble the active and inactive states, Li et al. observed that the key structural elements change in a highly concerted fashion during the conformational transition [34]. Similar to our study, Li et al. highlighted the importance of the rotameric transition of W246 during the activation process. In contrast to these studies, which focused more on the ligand-receptor interactions, we tried to identify the communication among the microswitches. To best of our knowledge, our simulation is the first systematic study probing the dynamics of all microswtich residues of A 2A AR, and made explicit that there are marked differences among the liganddependent conformational ensembles; and thus each ligand-dependent macrostate is made of mutually exclusive population of microstates, which supports the view that different type of ligands effectively remodel the protein energy landscapes [24,55].
In describing the conformational transition from inactive to active states in GPCRs, the recent studies by Yuan et al. [56] and Bhattacharya et al. [57] showed consistent and complementary results to our work, although the used methodologies were different. Yuan et al. discovered that a hydrophobic layer located inside of the helix structure forms a gate that opens to form a continuous water channel only upon the agonist binding [56]. The configuration of the nearby NPxxY motif affects this water channel, and the Ga peptide stabilizes the active state of Y7.53 which is one of the 10 switches in our study. Also, they calculated the average volume of the intracellular G-protein binding site and found that the volume is significantly increased upon the activation process, so that the intracellular site can accommodate G-protein binding. The volume change might well be a consequence of the swinging motion of TM5-ICL3-TM6, and consistent with our result showing that the distance between the intracellular loop 2 (ICL2) and ICL3 varies depending on the ligand binding states. Our study revealed the correlation between the ICL3 motion and the ionic-lock formation (Fig. 2b). Bhattacharya et al. calculated the mutual information (MI) in internal coordinates from MD simulated trajectories, and they identified allosteric communication pipelines which are cenceptually similar to the long-range cross-correlation pathways discussed in our previous work [21]. They found that, in the inactive state, the allosteric pipelines mainly cross the TM6 helix, and as the state moves from intermediate to agonist-bound active state, these pipelines pass though the TM7 helix, suggesting that TM7 is increasingly correlated in the active state. As suggested in our study, diverse conformational space of GPCRs is dependent on the ligand binding states and is regulated by the allosteric pathways comprising of some hub residues. The works by these two groups and ours both surmise TM7 as the key helix in GPCR activation. While the most dynamic part is TM5-ICL3-TM6 region which accompanies the swinging motion, the hub residues, scattered throughout the helices, regulate the activation process in a concerted way.
The antagonist-bound form explores the rotamer angle space entirely different from those in the agonist-bound and apo forms. The microstates visited by the apo form show an overlap with those by agonist-bound form although the extent of such overlap in the entire microstate population is vanishingly small and the duration of such overlap is only transient [58]. Occasional 10 o -tilt of TM5-ICL3-TM6 relative to TM3, shown in the simulation of the apo state (Fig. 2b), may be related to the basal activities of the GPCRs. According to the pre-coupled theory, several class A GPCRs are suspected to be coupled to their corresponding G-proteins even in the inactive or basal states [59][60][61][62]. Inactive-state preassembly can facilitate the rapid and specific G protein activation [61]. Our simulation results suggest that the apo state of A 2A AR can also form an inactive-state preassembly by visiting the microstates that overlap with those of the antagonist-bound state. We expect that a simulation of the apo form conducted in the presence of a G-protein will amplify this overlap and change the dynamics and correlation between the hotspots as well.
Among several findings and predictions made in the present study, of particular note is the role of W246 (S7 ) in the GPCR activation. Although the importance of W246 residue were largely documented based on experimental or theoretical studies [5,34,53], our cross-correlation analysis unequivocally indicates that W246 can work as a hub in the communications among the important microswitch residues (Fig. 7). Positioned deep inside the binding cleft, a signal of rotameric change of W246, triggered by a direct hydrophobic interaction with an agonist, can be transmitted to the 6 other switches (S1 , S2 , S3 , S5 , S8 , S9 ) that W246 is in direct correlation with. Our study confirms that W246 plays pivotal roles in GPCR activation as both an agonist sensor and actuator of allosteric signaling.
The class A GPCRs are conventionally reported to function as monomers [63][64][65][66]; however, growing body of experimental evidence indicates that some GPCRs, including A 2A AR, can form homodimers, heterodimers, or even higher level of oligomers [67][68][69]. In recent years, single-molecule imaging techniques revealed that GPCRs undergo dynamic equilibrium between monomers and dimers [70], and studies of GPCR monomers and dimers are both meaningful and necessary. The dynamic features and correlations revealed in this study for monomer of A 2A AR could change when the receptor forms dimers or higher level of oligomers. Thus, it would be interesting to investigate how dimerization of GPCR alters the correlations between hotspots.
Faithful description of biomolecular dynamics, as a complex system, is a major challenge in both computational and experimental molecular biology. In this study, we try to simplify the description of each ligand-bound macrostate by selecting the 10 binary switches as the reaction coordinates. The conformational features of A 2A AR captured in our 10 binary switch representation confirmed existing knowledge on the receptor and made specific predictions amenable to a further experimental study.

Methods
Preparation of the apo, antagonist-, and agonist-bound structures. The X-ray crystal structures of A 2A AR bound with an agonist or antagonist were retrieved from the Protein Data Bank (PDB). Since some loop regions are not resolved in these crystal structures, homology modeling was performed using the MODELLER program implemented in Discovery Studio v.3.1 to prepare the full-length A 2A AR models including all the loop regions. We used the mutation-free X-ray crystal structures of PDB IDs 3QAK [41] and 3EML [2], which were available at the year 2011 we began this study, as the main templates for the agonist-bound and antagonist-bound models, respectively. To model the loop regions that were not determined in 3QAK and 3EML, the X-ray crystal structures with PDB IDs of 2YDV [4] and 3PWH [54] were used by retaining the conserved disulfide bridges connecting the loops, i.e., C71-C159, C74-C146, C77-C166, and C259-C262. These models were optimized by a simulated annealing step and selected based on the Discrete Optimized Protein Energy (DOPE) score [71]. The final structures were energy-minimized with the Conjugate Gradient method and the Generalized Born with simple SWitching (GBSW) implicit solvent model [72]. We obtained the apo form by minimizing the receptor strucuture after removing the ligand from the antagonist-bound form.
Simulations. By employing X-ray crystal structures of A 2A AR and homology modeling to resolve the missing residues in the loop as described above, we performed MD simulations using NAMD package v.2.8 with CHARMM22/CMAP force field. The topology and parameter files for the ligands were generated by SwissParam web server [73].
To construct the explicit membrane system, we first predicted the TM region of the A 2A AR based on the Orientations of Proteins in Membranes (OPM) database, and subsequently surround the TM region with 173 (apo), 182 (antagonist-bound), and 177 (agonist-bound) 1-Palmitoyl-2-oleoylphosphatidylcholine (POPC) lipid molecules, which form a bilayer embedding the receptor (length: 88 Å in x-axis, 91 Å in y-axis). The receptor in the membrane system was then solvated with 13,549 (apo), 13,551 (antagonist-bound), and 13,552 (agonistbound) water molecules. We added 38 K + and 49 Cl − ions to make * 150 mM salt condition.
The simulations were conducted via (i) energy-minimization of the initial system using the cojugate gradient method in the order of membrane, water molecules, and the entire molecules; (ii) gradual heating from 0 K to 300 K using a 0.01 K interval at each step; (iii) 50 nsec equilibration with NVT ensemble; and (iv) * 1 msec production run with NPT ensemble for each system with different ligands. No specific constraints, such as distance or angle, were applied during the simulations, except SETTLE algorithm for constraints in water molecules. We used the integration time step of 1 fsec, and for the analysis saved the simulated trajectories every 2 psec and sampled every 400 psec.
To simulate the effect of inserting the agonist into the simulated apo-state, we first aligned the agonist ligand structure (co-crystallized conformation in 3QAK.pdb) to the apo-state structure, and then deleted the water molecules around 5 Å from the region where the ligand is to be inserted. After inserting the agonist, we performed the energy minimization and equilibrated the system for 20 nsec, so that water can solvate the ligand binding site as well as the agonist. The effect of agonist binding on the receptor structure was monitored subsequently.
Supporting Information S1 Fig. Hydrogen bonding network of the micro-switch residues in TM1, TM2, and TM7. The time traces are colored in black, blue, and red for the apo, antagonist-bound, and agonistbound states, respectively, and their histograms are shown in the right side of the graphs.