Abstract

# Communication over the network of binary switches regulates the activation of A$_{2A}$ adenosine receptor

Communication over the network of binary switches regulates the activation of A adenosine receptor

Yoonji Lee, Sun Choi, Changbong Hyeon

National Leading Research Laboratory (NLRL) of Molecular Modeling and Drug Design, College of Pharmacy, Graduate School of Pharmaceutical Sciences, and Global Top 5 Research Program, Ewha Womans University, Seoul 120-750, Korea

School of Computational Sciences, Korea Institute for Advanced Study, Seoul 130-722, Korea

E-mail: sunchoi@ewha.ac.kr (S.C.); hyeoncb@kias.re.kr (C.H.)

## Abstract

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 A adenosine receptor (AAR) 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 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 receptor in three distinct ligand-bound states using A adenosine receptor as a model system of GPCRs. Statistical analyses on the dynamics of and correlation among the 10 “binary switches” reveal that the three receptor states retain distinct dynamic properties. The antagonist- and agonist-bound forms of the receptors explore vastly different conformational space, and the apo form lies between them, yet located closer to the antagonist-bound form. In regard to the agonist-binding triggered activation mechanism, the correlation map among the 10 binary switches unequivocally shows that direct sensing of agonist ligand by the indole ring of W246 actuates the rest of intramolecular signaling.

## 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 residue-to-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 adenosine receptor; and 2.96 Å for adrenergic receptor whose active state structure is crystalized with G-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 long-range 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 adenosine receptor (AAR), 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 AAR, 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 sec all-atom molecular dynamics simulation of AAR in apo, antagonist-bound (ZM-241385),[40, 2] 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 microstates. By projecting the receptor dynamics on 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 cross-correlation 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 adenosine receptors could be useful for the discovery of effective drugs.

## Results

Among the 18 microswitches of AAR, 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, D52, D101, R102, W129, Y197, E228, W246, N284, P285, Y288 where the superscripts refer to the Ballesteros Weinstein numbering system [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 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) and forms a new H-bond between S281 and N280, 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 and HB as the switch 1 () and the switch 2 (), 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 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 ( Å) (Fig. 2a and the top panel of Fig. 2b); while it is rarely formed in the agonist-bound state ( Å). In the apo state, the ionic-lock maintains the identical distance with that of the antagonist-bound state, but occasionally breaks ( Å) and reforms ( Å) (Fig. 2 and Fig. S2).

The ionic-lock and the inter-residue distance () between L110 in ICL2 and A221 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 on the receptor state is clearly observed; Å (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 and (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) binds TM3 tightly with TM2. As HB and the ionic-lock in the DRY motif can be used to discern the ligand-dependent macrostate of GPCR, we define HB and the ionic-lock as and , 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 , , , , and , respectively: (i) In the agonist-bound state, W129 () shows asymmetric bimodal distribution of dihedral angle that has a major peak at and a minor peak at . In the apo and antagonist-bound forms, is unimodal with a peak at (Fig. 3a and Fig. S5a); (ii) The dihedral angles of Y197 () also show distinct distribution in the three receptor states (Fig. 3b and Fig. S5b); (iii) W246 (), known as a central switch located at the bottom of the orthosteric binding cleft [5], shows discrete transition between and with the receptor state-depedent angle value and frequency (Fig. 3c, Fig. S5c, and supporting movie M1); (iv) N284 () and Y288 () 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 and S281, 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 () and HB.

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 -th switch () 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, corresponding to 2, is disrupted in the agonist-bound form. In this case we set for disrupted H-bond and for intact H-bond. Also, for whose configuration is best described using rotameric angle, we consider to be in the ON state if the dihedral angle (C-C-C-C) of W246 (see Fig. S4) is equal or less than ; otherwise, it is in the OFF state (). 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].

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, possible microstates; each microstate is expressed using binary number from 0000000000 to 1111111111 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 microstates are shown in Fig. 5a and Fig. 5c, respectively. The average value of each switch, calculated in different receptor state as (Fig. 5b), where with 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 Fig. 6, (i) 80 % of antagonist-bound form are populated in the 0000000000 or 0000000001 state. (ii) 20.2 % of the switches in the agonist-bound form are in 1023th state (1111111111), and % are in 895th state (1101111111). (iii) Lastly, in the apo form of GPCR, on an average, , , are ON state, while , , , , are OFF state. Microstates that constitute the major population of the apo form are 1000000101 (31.32 %) and 1001000101 (10.56 %).

To quantify the statistical similarity explored by two different receptor states, say and () in the 10-switch representation, we employ Hamming distance:

 dαβ=10∑i=1|⟨sα,i⟩−⟨sβ,i⟩| (1)

where is the absolute value (or modulus) of , and is the average value of -th switch () in the receptor form , which is calculated in Fig. 5b. Since , it is expected that . The more similar, the smaller the value of should be. We obtain , , .

Next, the “complexity” of each macrostate, defined with the ensemble of microstates where , is quantified using the Shannon entropy:

 I=−Ns∑i=1pilog2pi (2)

where is the probability of occupying the -th microstate as in Fig. 6. When the receptor explores only a single state ( and ), the value of should be ; and if all the states are uniformly explored ( for all ) then . Thus, the larger the value of , the more diverse microstates are explored, which is also gleaned from Fig. 6. We obtain , , , 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 () by using the conformational ensemble from the simulations.

 Cij=⟨δsiδsj⟩√⟨(δsi)2⟩√⟨(δsj)2⟩ (3)

where is the variation of the switch value from its mean. assesses the extent of coherence in the “change” in switch dynamics between the -th and -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 , and many other negative correlations ; (ii) By contrast, in the antagonist-bound state (the left panel in Fig. 7), a positive correlation is detected only between the and , 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, (W246), a central rotameric switch located at the deep bottom of ligand binding cleft, displays direct couplings with 6 other switches , , , , , and , and additional positive correlations are observed in , , , . This suggests that intramolecular signaling over the entire structure can be initiated by stimulating the .

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 Å and Å (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 - 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 , , , , , , and 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 () and were simulated for ns. The second two traces (cases 3 & 4) were generated at 595 ns and 625 ns when the ionic-lock was disrupted (), and were simulated for 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 () is compared (Fig. 5b middle panel versus Fig. 8b). This is not so surprising given that time scale of our simulation ( sec) 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 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 (); in particular, the value of 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 (, 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 and (Fig. 8d, e) than the cases 1 and 2 (, ). 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 AAR.

## Discussion

Given that typical time scale associated with the transition from inactive to active GPCR conformers is 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 sec [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 sec), 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 sec 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 AAR 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 AAR 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, 41, 4, 54].

Until recently, there are some other studies that have compared the conformational differences of AAR depending on the receptor states [34, 35, 36]. Pang et al. studied dynamic behaviors of AAR induced by the binding of two distinct antagonists [35]. Lee et al. simulated the ligand-dependent cholesterol interactions in AAR [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 AAR, and made explicit that there are marked differences among the ligand-dependent 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 G 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-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 AAR 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 () in the GPCR activation. Although the importance of W246 residue were largely documented based on experimental or theoretical studies[5, 53, 34], 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 (1, 2, 3, 5, 8, 9) 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 AAR, 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 AAR 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 AAR 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 AAR 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 AAR 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 AAR 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 AAR 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 (agonist-bound) water molecules. We added 38 K and 49 Cl ions to make mM salt condition. As a result, the simulation system contains 68,799 (apo), 70,052 (antagonist-bound), and 69,449 (agonist-bound) atoms in 85 Å88 Å99 Å periodic box. Nonbonded interactions were smoothly switched off between 10 and 12 Å. To handle electrostatic interactions, the particle mesh Ewald algorithm was employed with a grid spacing smaller than 1 Å.

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 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.

## Acknowledgments

This work was partly supported by the grant from the Basic Science Research Program (Grant number 2013R1A6A3A01066055) (to Y.L.) and the National Leading Research Lab (NLRL) program (2011-0028885) (to S.C.) funded by the Ministry of Science, ICT & Future Planning (MSIP) and the National Research Foundation of Korea (NRF). C.H. thanks the KITP at the University of California, Santa Barbara, for support during the preparation of the manuscript (NSF PHY11-25915). We thank KIAS and KISTI Supercomputing Center for providing computing resources.

## Supporting Figures

### References

1. Rosenbaum D, Rasmussen S, Kobilka B (2009) The structure and function of G-protein-coupled receptors. Nature 459: 356–363.
2. Jaakola V, Griffith M, Hanson M, Cherezov V, Chien E, et al. (2008) The 2.6 angstrom crystal structure of a human A2A adenosine receptor bound to an antagonist. Science 322: 1211-1217.
3. Rasmussen S, DeVree B, Zou Y, Kruse A, Chung K, et al. (2011) Crystal structure of the adrenergic receptor-gs protein complex. Nature 477: 549–555.
4. Lebon G, Warne T, Edwards P, Bennett K, Langmead C, et al. (2011) Agonist-bound adenosine A2A receptor structures reveal common features of GPCR activation. Nature 474: 521–525.
5. Nygaard R, Frimurer T, Holst B, Rosenkilde M, Schwartz T (2009) Ligand binding and micro-switches in 7TM receptor structures. Trends Pharmacol Sci 30: 249–259.
6. Monod J, Changeux JP, Jacob F (1963) Allosteric proteins and cellular control systems. J Mol Biol 6: 306–329.
7. Edward WY, Koshland DE (2001) Propagating conformational changes over long (and short) distances in proteins. Proc Natl Acad Sci USA 98: 9517–9520.
8. Gandhi PS, Chen Z, Mathews FS, Di Cera E (2008) Structural identification of the pathway of long-range communication in an allosteric enzyme. Proc Natl Acad Sci USA 105: 1832–1837.
9. Weinkam P, Pons J, Sali A (2012) Structure-based model of allostery predicts coupling between distant sites. Proc Natl Acad Sci USA 109: 4875–4880.
10. Hyeon C, Lorimer GH, Thirumalai D (2006) Dynamics of allosteric transition in GroEL. Proc Natl Acad Sci USA 103: 18939-18944.
11. Horovitz A, Willison KR (2005) Allosteric regulation of chaperonins. Curr Opin Strct Biol 15: 646-651.
12. Hyeon C, Onuchic JN (2011) A Structural Perspective on the Dynamics of Kinesin Motors. Biophys J 101: 2749-2759.
13. Vale RD, Milligan RA (2000) The way things move: Looking under the hood of molecular motor proteins. Science 288: 88-95.
14. Whitford PC, Miyashita O, Levy Y, Onuchic JN (2007) Conformational transitions of adenylate kinase: Switching by cracking. J Mol Biol 366: 1661-1671.
15. Lee Y, Jeong LS, Choi S, Hyeon C (2011) Link between Allosteric Signal Transduction and Functional Dynamics in a Multisubunit Enzyme: S-Adenosylhomocysteine Hydrolase. J Am Chem Soc 133: 19807-19815.
16. Bahar I, Lezon TR, Yang LW, Eyal E (2010) Global Dynamics of Proteins: Bridging Between Structure and Function. Annu Rev Biophys 39: 23-42.
17. Zheng W, Brooks BR, Thirumalai D (2006) Low-frequency normal modes that describe allosteric transitions in biological nanomachines are robust to sequence variations. Proc Natl Acad Sci USA 103: 7664-7669.
18. Cui Q, Karplus M (2008) Allostery and cooperativity revisited. Protein Science 17: 1295–1307.
19. Popovych N, Sun S, Ebright RH, Kalodimos CG (2006) Dynamically driven protein allostery. Nature Struct Mol Biol 13: 831–838.
20. Balabin I, Yang W, Beratan D (2009) Coarse-grained modeling of allosteric regulation in protein receptors. Proc Natl Acad Sci USA 106: 14253–14258.
21. Lee Y, Choi S, Hyeon C (2014) Mapping the intramolecular signal transduction of G-protein coupled receptors. Proteins: Struct Func Bioinfo 82: 727-743.
22. Cooper A, Dryden DTF (1984) Allostery without conformational change. Eur Biophys J 11: 103–109.
23. Buchli B, Waldauer SA, Walser R, Donten ML, Pfister R, et al. (2013) Kinetic response of a photoperturbed allosteric protein. Proc Natl Acad Sci USA 110: 11725–11730.
24. Motlagh HN, Wrabl JO, Li J, Hilser VJ (2014) The ensemble nature of allostery. Nature 508: 331–339.
25. Itoh K, Sasai M (2010) Entropic mechanism of large fluctuation in allosteric transition. Proc Natl Acad Sci 107: 7775-7780.
26. Yao X, Parnot C, Deupi X, Ratnala VRP, Swaminath G, et al. (2006) Coupling ligand structure to specific conformational switches in the beta2-adrenoceptor. Nat Chem Biol 2: 417–422.
27. Hoffmann C, Zürn A, Bünemann M, Lohse M (2008) Conformational changes in G-protein-coupled receptorsâthe quest for functionally selective conformations is open. British J Pharmacol 153: S358–S366.
28. Ambrosio M, Zürn A, Lohse MJ (2011) Sensing G protein-coupled receptor activation. Neuropharmacology 60: 45–51.
29. Altenbach C, Yang K, Farrens DL, Farahbakhsh ZT, Khorana HG, et al. (1996) Structural features and light-dependent changes in the cytoplasmic interhelical E-F loop region of rhodopsin: a site-directed spin-labeling study. Biochemistry 35: 12470–12478.
30. Altenbach C, Klein-Seetharaman J, Hwa J, Khorana HG, Hubbell WL (1999) Structural features and light-dependent changes in the sequence 59-75 connecting helices I and II in rhodopsin: a site-directed spin-labeling study. Biochemistry 38: 7945–7949.
31. Niesen MJM, Bhattacharya S, Vaidehi N (2011) The Role of Conformational Ensembles in Ligand Recognition in G-Protein Coupled Receptors. J Am Chem Soc 133: 13197-13204.
32. Kohlhoff KJ, Shukla D, Lawrenz M, Bowman GR, Konerding DE, et al. (2014) Cloud-based simulations on Google Exacycle reveal ligand modulation of GPCR activation pathways. Nat Chem 6: 15–21.
33. Dror RO, Arlow DH, Maragakis P, Mildorf TJ, Pan AC, et al. (2011) Activation mechanism of the 2-adrenergic receptor. Proc Natl Acad Sci USA 108: 18684–18689.
34. Li J, Jonsson AL, Beuming T, Shelley JC, Voth GA (2013) Ligand-dependent activation and deactivation of the human adenosine A(2A) receptor. J Am Chem Soc 135: 8749–8759.
35. Pang X, Yang M, Han K (2013) Antagonist binding and induced conformational dynamics of GPCR A2A adenosine receptor. Proteins 81: 1399–1410.
36. Lee JY, Patel R, Lyman E (2013) Ligand-dependent cholesterol interactions with the human A(2A) adenosine receptor. Chem Phys Lipids 169: 39–45.
37. Freeman L (1979) Centrality in social networks conceptual clarification. Social networks 1: 215–239.
38. Goh KI, Kahng B, Kim D (2001) Universal behavior of load distribution in scale-free networks. Phys Rev Lett 87: 278701.
39. Del Sol A, Fujihashi H, Amoros D, Nussinov R (2006) Residues crucial for maintaining short paths in network communication mediate signaling in proteins. Mol Sys Biol 2: 2006.0019.
40. Dall’lgna OP, Porciúncula LO, Souza DO, Cunha RA, Lara DR (2003) Neuroprotection by caffeine and adenosine A2A receptor blockade of -amyloid neurotoxicity. British J Pharmcol 138: 1207–1209.
41. Xu F, Wu H, Katritch V, Han G, Jacobson K, et al. (2011) Structure of an agonist-bound human A2A adenosine receptor. Science 332: 322.
42. Ballesteros JA, Weinstein H (1995) Integrated methods for modeling G-protein coupled receptors. Methods Neurosci 25: 366–428.
43. Mono J, Wyman J, Changeux JP (1965) On the nature of allosteric transitions: a plausible model. J Mol Biol 12: 88-118.
44. Weis W, Kobilka B (2008) Structural insights into G-protein-coupled receptor activation. Curr Opin Struct Biol 18: 734–740.
45. Dahl ACE, Chavent M, Sansom MSP (2012) Bendix: intuitive helix geometry analysis and abstraction. Bioinformatics 28: 2193–2194.
46. Zwanzig R (1995) Simple model of protein folding kinetics. Proc Natl Acad Sci USA 92: 9801-9804.
47. Bissantz C, Kuhn B, Stahl M (2010) A Medicinal Chemistâs Guide to Molecular Interactions. J Med Chem 53: 5061–5084.
48. Nygaard R, Zou Y, Dror RO, Mildorf TJ, Arlow DH, et al. (2013) The Dynamic Process of -Adrenergic Receptor Activation. Cell 152: 532–542.
49. Vilardaga JP, Bünemann M, Krasel C, Castro M, Lohse MJ (2003) Measurement of the millisecond activation switch of G protein coupled receptors in living cells. Nat Biotech 21: 807–812.
50. Grossfield A, Feller SE, Pitman MC (2007) Convergence of Molecular Dynamics Simulations of Membrane Proteins. Proteins: Struct Funct Bioinfo 67: 31-40.
51. Hyeon C, Thirumalai D (2011) Capturing the essence of folding and functions of biomolecules using coarse-grained models. Nat Commun 2: 487.
52. Dror RO, Dirks RM, Grossman J, Xu H, Shaw DE (2012) Biomolecular simulation: a computational microscope for molecular biology. Annu Rev Biophys 41: 429–452.
53. Crocker E, Eilers M, Ahuja S, Hornak V, Hirshfeld A, et al. (2006) Location of Trp265 in metarhodopsin II: implications for the activation mechanism of the visual receptor rhodopsin. J Mol Biol 357: 163–172.
54. Doré AS, Robertson N, Errey JC, Ng I, Hollenstein K, et al. (2011) Structure of the Adenosine A Receptor in Complex with ZM241385 and the Xanthines XAC and Caffeine. Structure 19: 1283–1293.
55. Smock RG, Gierasch LM (2009) Sending signals dynamically. Science 324: 198–203.
56. Yuan S, Filipek S, Palczewski K, Vogel H (2014) Activation of G-protein-coupled receptors correlates with the formation of a continuous internal water pathway. Nat Commun 5: 4733.
57. Bhattacharya S, Vaidehi N (2014) Differences in Allosteric Communication Pipelines in the Inactive and Active States of a GPCR. Biophys J 107: 422–434.
58. Yao XJ, Ruiz GV, Whorton MR, Rasmussen SGF, DeVree BT, et al. (2009) The effect of ligand efficacy on the formation and stability of a GPCR-G protein complex. Proc Natl Acad Sci USA 106: 9501–9506.
59. Galés C, Van Durm JJ, Schaak S, Pontier S, Percherancier Y, et al. (2006) Probing the activation-promoted structural rearrangements in preassembled receptor–G protein complexes. Nat Struct Mol Biol 13: 778–786.
60. Philip F, Sengupta P, Scarlata S (2007) Signaling through a G Protein-coupled receptor and its corresponding G protein follows a stoichiometrically limited model. J Biol Chem 282: 19203–19216.
61. Qin K, Dong C, Wu G, Lambert NA (2011) Inactive-state preassembly of G(q)-coupled receptors and G(q) heterotrimers. Nat Chem Biol 7: 740–747.
62. Jakubík J, Janíčková H, Randáková A, El-Fakahany EE, Doležal V (2011) Subtype differences in pre-coupling of muscarinic acetylcholine receptors. PLoS one 6: e27732.
63. Maurice P, Kamal M, Jockers R (2011) Asymmetry of GPCR oligomers supports their functional relevance. Trends Pharmacol Sci 32: 514–520.
64. Ernst OP, Gramse V, Kolbe M, Hofmann KP, Heck M (2007) Monomeric G protein-coupled receptor rhodopsin in solution activates its G protein transducin at the diffusion limit. Proc Natl Acad Sci U S A 104: 10859–10864.
65. Filizola M (2010) Increasingly accurate dynamic molecular models of G-protein coupled receptor oligomers: Panacea or Pandora’s box for novel drug discovery? Life Sci 86: 590–597.
66. Whorton MR, Bokoch MP, Rasmussen SGF, Huang B, Zare RN, et al. (2007) A monomeric g protein-coupled receptor isolated in a high-density lipoprotein particle efficiently activates its g protein. Proc Natl Acad Sci U S A 104: 7682–7687.
67. Terrillon S, Bouvier M (2004) Roles of G-protein-coupled receptor dimerization. EMBO Rep 5: 30–34.
68. Vidi PA, Chen J, Irudayaraj JMK, Watts VJ (2008) Adenosine A(2A) receptors assemble into higher-order oligomers at the plasma membrane. FEBS Lett 582: 3985–3990.
69. Fanelli F, Felline A (2011) Dimerization and ligand binding affect the structure network of A(2A) adenosine receptor. Biochim Biophys Acta 1808: 1256–1266.
70. Kasai RS, Kusumi A (2014) Single-molecule imaging revealed dynamic GPCR dimerization. Curr Opin Cell Biol 27: 78–86.
71. Shen My, Sali A (2006) Statistical potential for assessment and prediction of protein structures. Protein Sci 15: 2507–2524.
72. Im W, Lee MS, Brooks CL (2003) Generalized born model with a simple smoothing function. J Comp Chem 24: 1691–1702.
73. Zoete V, Cuendet MA, Grosdidier A, Michielin O (2011) Swissparam: a fast force field generation tool for small organic molecules. J Comput Chem 32: 2359–2368.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters