Energy landscapes reveal agonist control of GPCR activation via microswitches

Agonist binding to the extracellular part of G protein-coupled receptors (GPCRs) leads to conformational changes in the transmembrane region that activate cytosolic signalling pathways. Although several high resolution structures of the inactive and active receptor states are available, the atomistic details of the allosteric coupling that transmits the signal across the membrane are not fully understood. We calculated free energy landscapes of β2 adrenergic receptor activation using atomistic molecular dynamics simulations in an optimized string of swarms framework, which sheds new light on the roles of microswitches in governing the equilibrium between conformational states. Contraction of the extracellular binding site in the presence of agonist is obli-gatorily coupled to conformational changes in a connector motif located in the core of the transmembrane region. In turn, the connector is probabilistically coupled to the conformation of the intracellular region: an active connector promotes desolvation of a buried solvent-filled cavity, a twist of the conserved NPxxY motif, and an interaction between two conserved tyrosines in transmembrane helices 5 and 7 (Y-Y motif), which leads to a larger population of active-like states at the G protein binding site. This coupling is further augmented by protonation of the strongly conserved Asp792.50, which locks the solvent cavity, NPxxY, and Y-Y motifs in active-like conformations. The agonist binding site hence communicates with the intracellular region via a cascade of locally connected switches and characterizing the free energy landscapes along the conformation of these microswitches contributes to understanding of how ligands can stabilize distinct receptor states. We demonstrate that the developed simulation protocol is transferable to other class A GPCRs and anticipate that it will become a useful tool in the design of drugs with specific signaling properties.


Introduction
G protein-coupled receptors (GPCRs) are membrane proteins that activate cellular signaling pathways in response to extracellular stimuli. There are more than 800 GPCRs in the human genome 1 and these recognize a remarkably large repertoire of ligands such as neurotransmitters, peptides, proteins, and lipids. This large superfamily plays essential roles in numerous physiological processes and has become the most important class of drug targets 2 .
All GPCRs share a common architecture of seven transmembrane (TM) helices, which recognize the cognate ligand in the extracellular region and triggers intracellular signals via a more conserved cytosolic domain ( Fig. 1) 3,4 . GPCRs are inherently flexible proteins that exist in multiple conformational states and drug binding alters the relative populations of these. Agonists will shift the equilibrium towards active-like receptor conformations, which promote binding of G proteins and other cytosolic proteins (e.g. arrestin), leading to initiation of signaling via multiple pathways. In the apo state, GPCRs can still access active-like conformations and thereby exhibit a smaller degree of signaling, which is referred to as basal activity.
Breakthroughs in structure determination of GPCRs during the last decade have provided insights into the process of activation at atomic resolution (Fig. 1). In particular, crystal structures of the β 2 adrenergic receptor (β 2 AR) in both active and inactive conformational states [5][6][7][8][9][10][11][12][13][14] have revealed hallmarks of GPCR activation. The observations made for this prototypical receptor have recently been reinforced by crystal and cryogenic electron microscopy (cryo-EM) structures for other family members 15 . The most prominent features of GPCR activation are a large ∼ 1.0 − 1.4 nm outward movement of TM6 and a slight inward shift of TM7 on the intracellular side ( Fig. 1), which create a large cavity for binding of cytosolic proteins.
Conserved changes in the extracellular part are more difficult to discern due to the lower sequence conservation in this region. In general, the structural changes are relatively subtle and only involve a small contraction of the orthosteric site [5][6][7] . In the case of the β 2 AR, the catechol group of adrenaline forms hydrogen bonds with Ser207 5.46 (superscripts denote Ballesteros-Weinstein numbering 16 ), which leads to a ∼ 0.2 nm inward bulge of TM5. These structural changes then propagate through the receptor via several conserved motifs. The rearrangement of TM5 influences a connector region (PI 3.40 F 6.44 motif), which is in contact with the highly conserved Asp79 2.50 and NP 7.50 xxY 7.53 motif via a water-filled cavity. Upon activation, the cavity dehydrates and the NPxxY motif twists, which reorients Tyr326 7.53 to a position where it can form a water-mediated interaction with Tyr219 5.58 ("Y-Y motif") 17 and enables formation of the G protein binding site. Characterization of the role these individual microswitches play in determining the conformational ensemble of structures could guide the design of drugs with tailored signaling profiles.
Indeed, the allosteric control of GPCR activation by extracellular ligands cannot be fully understood from the static structures captured by crystallography or cryo-EM. Mutagenesis and spectroscopy studies [18][19][20][21] have suggested that the efficacy of a ligand is determined by a complex interplay between different microswitches and population of distinct states lead to specific functional outcomes. Molecular dynamics (MD) simulations are well suited to study the conformational landscape of GPCRs as this method can provide an atomistic view of the flexible receptor in the presence of membrane, aqueous solvent, and ligands [22][23][24][25][26][27][28][29][30] . The Figure 1: The activation mechanism of GPCRs involves a series of microswitches. Agonist binding leads to an inward bulge of TM5 (quantified as the distance between Ser207 5.46 and Gly315 7.41 ), which leads to a conformational change in the connector region (Ile121 3.40 and Phe282 6.44 ). The transmembrane cavity surrounding Asp79 2.50 dehydrates and the NPxxY motif twists upon activation, leading the Y-Y interaction to form. TM6 moves outwards to create the G protein binding site. The inactive (PDB ID 2RH1 5 ) and active (PDB ID 3P0G 7 ) structures are represented in white and orange, respectively. seminal paper by Dror et al 27 provided insights into the activation mechanism of the β 2 AR by monitoring how a crystallized active state conformation relaxed to an inactive state upon removal of the intracellular binding partner, a G protein mimicking nanobody, using MD simulations. Although key conformational changes involved in the transition from active to inactive conformations were identified in these simulations, this approach has inherent limitations. Indeed, understanding the roles of different microswitches in activation and the strength of the coupling between these requires mapping of the relevant free energy landscapes, which are still too costly to calculate using brute-force MD simulations. Enhanced sampling methods, on the other hand, provide a means to explore the conformational landscapes of proteins to a relatively low computational cost 31 and have been exploited to study the activation mechanism of GPCRs [22][23][24]26,[28][29][30] .
In this work, we aimed to identify the most probable path describing the transition between inactive and active states of β 2 AR and to characterize the full conformational ensemble along this pathway at an atomistic level of resolution. A new version of the string method with swarms of trajectories was first developed for this purpose. The free energy landscapes associated with β 2 AR activation revealed that whereas agonist binding is only loosely connected to outward movement of TM6, the coupling between microswitch pairs in the transmembrane region ranges from weak to strong and is influenced by the protonation of Asp79 2.50 . Finally, we demonstrated that our approach can be transferred to study free energy landscapes of any class A GPCR at a modest computational cost.

Results
Optimizing the enhanced sampling protocol for increased efficiency We aimed to compute the most probable transition pathway linking the inactive and active states of β 2 AR and the relative free energy of the states lining this pathway. For this purpose we chose the string of swarms method 32 . In this framework, the minimum free energy path in an N-dimensional collective variable (CV) space is estimated iteratively from the drift along the local gradient of the free energy landscape (Fig. S1). From each point on the string, a number of short trajectories are launched (a swarm), which enables us to calculate the drift. The string is then updated considering this drift and reparametrized to ensure full sampling of the configurational space along the pathway. Convergence is reached when the string diffuses around an equilibrium position. The method allows to sample a highdimensional space at a relatively inexpensive computational cost since it only samples along the one-dimensional path of interest.
We characterized the pathway linking equilibrated conformations originating from active and inactive structures of the the β 2 AR (PDB IDs 3P0G and 2RH1, respectively), adding two short strings spanning the active and in the inactive regions to increase sampling of the end state environments (see Material and Methods and SI Methods). First, we characterized a CV set that embeds receptor activation by automatically inferring the inter-residue distances that are related to activation. Based on one of Dror et al.'s simulation trajectories of spontaneous deactivation of the β 2 AR 27 , we identified metastable states by clustering simulation configurations, followed by classification of these by training a fully connected neural network to identify states 33 . The most important input features for classification were identified via deep Taylor decomposition 34 and taken as CVs ( Fig. 2 and S2). The set of CVs we inferred was a network of interatomic distances between all seven TM helices ( Fig. 2). Encouragingly, this five-dimensional CV set captured the main degrees of freedom implicated in β 2 AR activation, including the large outward movement of TM6 and a smaller inward shift of TM7 at the intracellular face of the receptor.
To speed up convergence of the string optimization, we initiated our string simulation from a rough guess of the minimum free-energy landscape. The latter was obtained by estimating the density of points from Dror et al.'s trajectory in this CV space (Fig. S3). We also introduced algorithmic improvements to the string of swarms method: we adaptively chose the number of trajectories in a swarm, gradually increased the number of points on the string and introduced a reparametrization algorithm that improves performance as well as promotes exchanges of configurations between adjacent string points (Supplementary Methods, Fig. S1, S4). We carried out 300 iterations of string optimization for each system, considering a number of points on the string ranging from 20 to 43 and swarms consisting of 16 to 32 10 ps trajectories. Because we sought to only sample the vicinity of the most probably activation path, derivation of a converged free energy landscape required a mere total of ∼2-3 µs simulation time (Fig. S5, S6, S7, S8 and S9).

Minimum free energy pathway of β 2 AR activation
We derived the most probable transition path (   The CVs corresponding to distances TM2-TM7,  TM6-TM4, TM7-TM4, TM3-TM6 and TM6-TM5 defined in Table S1, are shown as purple, blue, green, yellow and red dashed lines, respectively. The change of these distance CVs from the inactive to the active state structures is reported in nm. a b Figure 3: The free energy landscapes are projected onto the first two CVs used to optimize the minimum free energy path. The minimum free energy pathways were optimized in the absence (a) and in the presence of a bound agonist (b). Characteristic events occurring during activation (see Table S2 for their definition) are located on the string. The active and inactive labels mark the regions close to the active and inactive structures 3P0G 7 and 2RH1 5 .
Binding of the agonist ligand modifies the order in which the helices rearrange, as can be seen when projecting the minimum free energy path along various CVs ( Fig. 3 and S5). The most probable activation pathway begins with an outward movement of TM6, followed by the inward shift of TM7. The connector and the bulge in TM5 are locked in their active-like states in the presence of agonist, which enables large outward movements of TM6. Similarly to the apo case, to reach the fully active conformation, the Asp79 2.50 cavity dehydrates and the Y-Y interaction forms shortly before the NPxxY twists.

Coupling between orthosteric and G-protein binding sites
As the final configurations of the string are at equilibrium in all dimensions, the trajectories from the last iterations can be used to compute the free energy landscape as a function of any variables (SI methods). This allowed a detailed analysis of how conformational changes induced by agonist binding propagate through the receptor to the G protein binding site.
The roles of conserved microswitches were assessed by comparing free energy landscapes in the absence and presence of a bound agonist. We first evaluated how the TM5 bulge, which reflects how the binding site contracts upon activation, influences the intracellular distance between TM6 and TM3 (Fig. 4a,b). In the absence of bound agonist, the receptor accessed both active and inactive conformations of the binding site, with the minimum of the free energy located close to the inactive state distance between TM3 and TM6. The TM3-TM6 distance could be as large as 1.5 nm when the ligand binding site was in the active conformation, an observation compatible with basal activity. Agonist binding led to the stabilization of the TM5 bulge in the orthosteric site (Fig. 4b). Although both inactive and active conformations remained accessible in the presence of the ligand, the minimum of the TM3-TM6 distance was shifted towards a more active-like state. A remarkable longrange allosteric coupling (>2 nm) between the orthosteric and G protein binding sites was hence captured by our simulations. The 0.5 nm outward movement of TM6 at the minimum of the landscape (Fig. 4a,b) is smaller than that observed in active crystal structures, in agreement with experiments demonstrating that the fully active conformation can only be stabilized in the presence of an intracellular partner 19,35 .

Propagation of activation through microswitches
Structural changes in the orthosteric site of the β 2 AR (Fig. 4s) have been proposed to propagate towards the intracellular part via a "connector" centered around Ile121 3.40 and Phe282 6.44 7 (Fig. 4t). Whereas we found the TM5 bulge and outward movement of TM6 to be loosely coupled, the free energy landscapes demonstrated that changes in the orthosteric site has a direct influence on the connector region (Fig. 4d,e). In the absence of agonist, both inactive and active conformations of the connector were populated and the connector region could assume an active-like state even if the TM5 bulge was inactive (Fig. 4d). In contrast, agonist binding resulted in a single free energy minimum where both the TM5 bulge and the connector were constrained to their active conformations (Fig. 4e). From the connector region, activation is propagated via several conserved motifs (Fig.   4g,h,j,k,m,n and u) 3 . To investigate the communication between microswitches in the core of the TM region, we analyzed if the connector was coupled to solvation of the cavity surrounding Asp79 2.50 , to the conformation of the NPxxY motif, and to the conformation of the Y-Y motif. In the apo state, an inactive connector region was tightly coupled to a hydrated cavity with ∼ 10 − 17 waters (Fig. 4g) and an inactive NPxxY motif (Fig. 4j). This result is consistent with a high-resolution crystal structure of an inactive β 1 adrenergic receptor (PDB ID 4BVN), in which a solvent network in this region was identified 36 . In contrast, the active conformation of the connector was compatible with both the fully hydrated cavity and a desolvated state with only 4 to 5 water molecules as well as an inactive and active NPxxY motif. Interestingly, the connector configuration was more tightly coupled to the Y-Y motif. An inactive connector was coupled to a broken Y-Y interaction whereas it was in intermediate state or fully formed with an active connector (Fig. 4m).
The free energy landscapes suggested that binding of agonist reduce the number of inactive-like states of the Asp79 2.50 cavity, NPxxY, Y-Y motifs (Fig. 4h,k,n) as well as the orientation of the Met82 2.53 (Fig. S9), a residue located one helix turn above Asp79 2.50 that has been the focus of NMR studies 37 . For the Asp79 2.50 cavity and NPxxY motif, both active and inactive conformations were accessible in the presence of agonist, but the active conformations were more favored energetically (Fig. S11). The presence of agonist resulted in an energy landscape with one active-like and two intermediate distances of the Y-Y motif. The two latter minima were shifted by approximately 0.1 nm towards the active conformation compared to the apo state (Fig. 4n).
The final combination of microswitches connected the Y-Y motif to the motion of TM6 (Fig. 4p,q). In the apo state, there was a relatively tight coupling between the Y-Y interaction and the TM3-TM6 distance, with several metastable states lining the minimum free energy pathway between inactive and active conformations (Fig. 4p). Ligand binding did not alter the coupling between these two microswitches, but generally tilted the free energy landscape towards more active states along both of these dimensions (Fig. 4q).
Having computed the free energy landscapes in absence and presence of ligand thus enables us to characterize the involvement of each microswitch in the transmission of allosteric coupling between the orthosteric ligand and the intracellular partner binding sites.

Impact of Asp79 2.50 protonation
Protonation of two ionizable residues, Asp79 2.50 and Asp130 3.49 , has been proposed to be involved in GPCR activation 27,38,39 . In particular, MD simulations have indicated that Asp79 2.50 , the most conserved residue among class A GPCRs, has a pK a value close to physiological pH and that the ionization state of this residue changes upon activation 38,40 .
As a previous simulation study of the β 2 AR showed that Asp79 2.50 (but not Asp130 3.49 ) may alter the activation pathway 27 , we repeated the calculations of the minimum free energy pathway of activation with this residue in its protonated (neutral) form (Fig. S7).
The free energy landscapes describing changes in the orthosteric site were similar to those obtained in simulations with Asp79 2.50 ionized (Fig. 4c,f). There was a weak coupling between the TM5 bulge and the intracellular region, with two major energy wells describing the conformation of TM6. Compared to the agonist-bound receptor with Asp79 2.50 ionized, the minima were shifted further towards active-like conformations for the protonated state ( Fig. 4c). The TM5 bulge remained strongly coupled to conformational changes observed in the connector region irrespective of the ionization state of Asp79 2.50 (Fig. 4e,f). The largest effects of Asp79 2.50 protonation were observed for the hydrated cavity surrounding this residue (Fig. 4h,i), NPxxY (Fig. 4k,l), and Y-Y motif (Fig. 4n,o): whereas the free energy landscapes showed that both active-and inactive-like conformations of these switches were populated in simulations with ionized Asp79 2.50 , the protonated state only resulted in energy wells close to the active conformation ( Fig. 4i,l,o). It was also evident that TM6 was stabilized in more active-like conformations by the protonated Asp79 2.50 (Fig. 4r).
Interestingly, intermediate states in which the Y-Y motif, as well as the NPxxY, adopted an active-like conformation and a low TM6 displacement were also populated ( Fig. 4r and   S11). Such an intermediate state has been observed, albeit rarely, in simulations considering a protonated Asp79 2.50 27 .
Comparison of representative structures from the simulations of ionized and protonated

Transferability of the methodology to other GPCRs
Efficient characterization of free energy landscapes with the string method relies on selection of appropriate CVs, which is a non-trivial task. Here, CVs were derived from a conventional MD trajectory of β 2 AR deactivation in a data-driven fashion. Considering that the conformational changes involved in class A GPCR activation are largely conserved 3 , we explored the possibility of transferring the CVs to the conformational sampling of other GPCRs. We mapped the CVs identified for β 2 AR to ten other class A GPCRs with active and inactive structures available (Fig. 6). Strikingly, the active and inactive structures clearly separate in two distinct clusters. This indicates that these CVs can describe the activation of many class A GPCRs and the protocol presented herein can be applied to other family members.   Table S3. Their clustering into two groups highlights that the activation mechanism of all class A GPCRs can likely by described by the CVs determined herein. and on methodologies suggested in previous molecular simulation studies 23,25,26,29 to further allow to assess the impact of agonist binding on microswitches involved in activation.

Discussion
Our work recapitulates a number of key findings of previous enhanced sampling 23,25,26,28,29 and long-timescale scale simulations of GPCRs 27 . In agreement with previous work, our simulations revealed a complex conformational landscape of the β 2 AR with only weak coupling between the othosteric site and G protein binding site. Through the analysis of the free energy landscapes we could quantify the coupling between spatially connected microswitches in the transmembrane region, which ranged from strong to relatively weak and was influenced by protonation of Asp79 2.50 . Our study also further illustrates that the energy landscapes depend on the variables chosen to project the conformational states. This is not an artifact of the protocol but rather an inherent limitation of dimensionality reduction. In particular, it is clear that considering only three states along the activation path (an active, intermediate and inactive state) does not allow to capture the complexity of the conformational changes induced by ligand binding 22,23,28,29 . The present work, through the resolution of the complete conformational ensemble lining the most probable activation pathway, can aid in determining the placement of spectroscopic probes to monitor the distribution of specific microswitches under different conditions.
Our work shows agreement with spectroscopy studies although a quantitative compar-ison remains challenging: whereas the β 2 AR crystallized in complex with agonists in the absence of intracellular partner (e.g. G protein or G protein-mimicking nanobody) are similar to identical to those determined with antagonists 45   On the other hand, a well-known limitation with the string with swarms method is that it only guarantees to converge to the most probable path closest to the initial path guess, and not necessarily the globally most probable path. The naive assumption of a straight initial path is not guaranteed to converge to the latter. Here we have proposed to alleviate this shortcoming by exploiting previous knowledge of the activation pathway and deriving an initial guess of the pathway likely to be close to the globally most probable path. An initial pathway can also be transferred from a similar system (as revealed in Fig. 6 Despite the major progress in structural biology for GPCRs, many aspects of receptor function remain misunderstood. Insights from atomistic MD simulations will continue to be valuable tools for interpreting experimental data. We expect our methodology to allow further insights into how binding of agonists influences the conformational landscape, poten-tially making it possible to design ligands with biased signalling properties 21 . The method is equally well suited to study the effect of allosteric modulators and the influence of different lipidic environments 54 . As the same approach can be transferred to other class A GPCRs, future applications will shed light on the common principles of activation as well as on the details that give each receptor a unique signalling profile, paving the way to the design of more effective drugs.

Methods
All swarm of trajectories simulations were instantiated with the coordinates from the 3P0G structure (the first two simulation systems in Table S4). The Asn187Glu mutation in 3P0G was reverted and Glu122 3.41 was protonated due to its localization in a particularly hydrophobic pocket, as has been common practice in other simulations of this receptor 45 . Residues To identify CVs, we performed clustering 61 on the frames from unrestrained MD simulation trajectory of β 2 AR (condition A in Dror et al.). The CVs were selected by training a multilayer perceptron classifier 62 using as input all the inter-residue distances and as output the cluster ID, followed by using Deep Taylor decomposition 34 to find key distances that could discriminate between clusters.
The endpoints of the main strings describing the transition between inactive and active states (subscript t ) were fixed to the output coordinates of equilibrated structures 2RH1 5 and 3P0G 7 . The initial path for simulations Holo t was guessed using data from Dror et al.: a rough estimate of the free energy landscape was calculated from the probability density landscape estimated using the Scikit-learn 62 kernel density estimator with automatic bandwidth detection (Fig. S3). Two additional short strings were set up to increase sampling in the active (subscript a ) and inactive (subscript i ) regions. The average, partially converged path between iterations 20-30 from Holo t was used as input path for simulations Apo t and Holo Ash79,t . All active (Apo a , Holo a and Holo Ash79,a ) and inactive (Apo i , Holo i and Holo Ash79,i ) substrings were initiated as straight paths between the endpoints. Even though the system was initiated from the nanobody bound structure 3P0G, the active state substring sampled conformations closer to the G protein bound structure 3SN6 (Fig. S12) 6 .
The swarms of trajectories simulations with optimizations (see SI and Fig. S4) were run for 300 iterations, at which point the strings had not changed on average for many iterations ( Fig. S5 and S7) and posterior distribution of free energy profiles given the data was narrow (Fig. S6).
By discretizing the system into microstates, or bins, it is possible to use the short trajectories from the swarms to create a transition matrix and derive the free energy distribution of the system 63 along a given variable (Fig. S10). In practice, the transition probabilities T ij of the transition matrix T can be estimated from the normalized number of transitions, N ij , from bin i to bin j: T ij = N ij / k N ik . The transition matrix of a physical system at equilibrium is constrained by detailed balance, such that for the stationary probability distribution, ρ: ρ i T ij = ρ j T ji . A Metropolis Markov chain Monte Carlo (MCMC) method was used to sample over the posterior distribution of transition matrices, given the unregularized elements of T ij 64 , and thereby obtain a distribution of free energy profiles for ρ (Fig. S6, S8 and S9). All code to run the simulations and reproduce the results in this paper is available for download 65 .