Dynamic allosteric networks drive adenosine A1 receptor activation and G-protein coupling

  1. Miguel A Maria-Solano
  2. Sun Choi

  • Curated by eLife

    eLife logo

    eLife assessment

    The authors describe the dynamics underlying allostery of the adenosine A1 receptor, providing valuable insights into the receptor's activation pathway. The enhanced sampling molecular dynamics simulations of available structural data, followed by network analysis, reveal transient conformational states and communication between functional regions. The authors carefully state the limitations of their work, including the restricted convergence of the free energy landscape and missing water-mediated hydrogen bond coordination. Collectively, they provide a convincing framework for advancing rational design strategies of specific modulators with desired modes of action.

    [Editors' note: this was originally reviewed and assessed by Biophysics Colab]

  • Curated by Biophysics Colab

    Biophysics Colab logo

    Evaluation statement (16 June 2023)

    Maria-Solano and Choi present the dynamics underlying allostery of the adenosine A1 receptor, providing valuable insights into the receptor's activation pathway. The enhanced sampling molecular dynamics simulations of available structural data, followed by network analysis, reveal transient conformational states and communication between functional regions. The authors carefully state the limitations of their work, including the restricted convergence of the free energy landscape and missing water-mediated hydrogen bond coordination. Collectively, the findings provide a convincing framework to advance rational design strategies of specific modulators with desired modes of action.

    Biophysics Colab considers this to be a convincing study and recommends it to scientists interested in the structural dynamics, allosteric pathway activations, and free energy landscapes of GPCRs.

    (This evaluation by Biophysics Colab refers to version 5 of this preprint, which has been revised in response to peer review of versions 3 and 4.)

Reviewed by

Read the full article See related articles

Listed in

Log in to save this article

Abstract

G-protein coupled receptors (GPCRs) present specific activation pathways and signaling among receptor subtypes. Hence, an extensive knowledge of the structural dynamics of the receptor is critical for the development of therapeutics. Here, we target the adenosine A 1 receptor (A 1 R), for which a negligible number of drugs have been approved. We combine molecular dynamics simulations, enhanced sampling techniques, network theory, and pocket detection to decipher the activation pathway of A 1 R, decode the allosteric networks, and identify transient pockets. The A 1 R activation pathway reveals hidden intermediate and pre-active states together with the inactive and fully-active states observed experimentally. The protein energy networks computed throughout these conformational states successfully unravel the extra and intracellular allosteric centers and the communication pathways that couple them. We observe that the allosteric networks are dynamic, being increased along activation and fine-tuned in the presence of the trimeric G-proteins. Overlap of transient pockets and energy networks uncovers how the allosteric coupling between pockets and distinct functional regions of the receptor is altered along activation. Through an in-depth analysis of the bridge between the activation pathway, energy networks, and transient pockets, we provide a further understanding of A 1 R. This information can be useful to ease the design of allosteric modulators for A 1 R.

Article activity feed

  1. Version published to 10.7554/elife.90773.2 on eLife
  2. Version published to 10.7554/elife.90773 on eLife
  3. Version published to 10.7554/elife.90773.1 on eLife
  4. eLife

    eLife assessment

    The authors describe the dynamics underlying allostery of the adenosine A1 receptor, providing valuable insights into the receptor's activation pathway. The enhanced sampling molecular dynamics simulations of available structural data, followed by network analysis, reveal transient conformational states and communication between functional regions. The authors carefully state the limitations of their work, including the restricted convergence of the free energy landscape and missing water-mediated hydrogen bond coordination. Collectively, they provide a convincing framework for advancing rational design strategies of specific modulators with desired modes of action.

    [Editors' note: this was originally reviewed and assessed by Biophysics Colab]

  5. eLife

    Joint Public Review:

    The objectives of the study:

    This paper aims to characterize the dynamics that drive allostery of the adenosine A1 receptor (A1R) via computational analysis of its activation free energy landscape and measurements of the appropriate geometrical parameters. This is done by focusing on the allosteric signaling pathways in different activation states, from inactive to active states via intermediate and pre-active ones, as well as the characterization of putative drug-binding pockets. The long-term objectives are to eventually be able to aid drug discovery efforts for this therapeutically important
    GPCR.

    Key findings and major conclusions:

    Conventional MD does not enable the sampling of the complete conformational landscape of receptor activation. Instead, enhanced sampling MD simulations are required to achieve this. Using metadynamics, the authors decipher the activation pathway of A1R, decode the allosteric networks and identify transient pockets. The protein energy networks computed throughout the inactive, intermediate active, pre-active and active conformational states unravel the extra and intracellular allosteric centers and the communication pathways that couple them, whereby the pathways are reinforced in the activated state. These conformations primarily differ in the dynamics of the ionic lock motif that couples TM3 to TM6 in the inactive conformation and reveal that G-proteins are required to fully stabilize the active conformation. Support for these findings comes from prior mutagenesis work on the A1R that identified key allosteric residues that in many cases map to identified communication nodes. Finally, the authors identified allosteric pockets throughout the A1R in four different conformational states that support prior experimental and MD studies on the mechanism of the positive allosteric modulator MIPS521 and which could be targeted for the design of new modulators. This indicates how energy networks are enhanced and redistributed by allosteric modulators and how this might explain their effect on receptor activity. Overall, these findings provide complementary support to a structure-based mechanism of activation and allosteric modulation of A1R, and extend the findings to incorporate dynamics across the full activation pathway.

    The perceived strengths and weaknesses:

    This preprint employs a combination of computational techniques to successfully reconstruct and analyze the conformational ensemble of the A1R activation. The metadynamics simulations supported the aim of the study, the results are clearly presented, and the work is very well written. The authors provide a valuable discussion of how the protein energy network analysis can contribute to the rational design of specific A1R modulators with desired mode of action. The employed computational approach does not capture communication pathways that involve water-mediated connections or interactions between ligands and residues. Moreover, full convergence of the free energy landscapes is not guaranteed. Overall, A1R is a good choice as the target for this study as there is existing structural and pharmacological data to support preliminary findings. Moreover, the framework presented herein could be adapted and scaled to other GPCRs with structural templates, which might enable comparison of allosteric pathways across families and classes.

  6. Version published to 10.1101/2022.08.08.502889v6 on bioRxiv
  7. Biophysics Colab

    Evaluation statement (16 June 2023)

    Maria-Solano and Choi present the dynamics underlying allostery of the adenosine A1 receptor, providing valuable insights into the receptor's activation pathway. The enhanced sampling molecular dynamics simulations of available structural data, followed by network analysis, reveal transient conformational states and communication between functional regions. The authors carefully state the limitations of their work, including the restricted convergence of the free energy landscape and missing water-mediated hydrogen bond coordination. Collectively, the findings provide a convincing framework to advance rational design strategies of specific modulators with desired modes of action.

    Biophysics Colab considers this to be a convincing study and recommends it to scientists interested in the structural dynamics, allosteric pathway activations, and free energy landscapes of GPCRs.

    (This evaluation by Biophysics Colab refers to version 5 of this preprint, which has been revised in response to peer review of versions 3 and 4.)

  8. Biophysics Colab

    Authors' Response (2 June 2023)

    GENERAL ASSESSMENT

    The objectives of the study: This paper aims to characterize the dynamics that drive allostery of the adenosine A1 receptor (A1R) via computational analysis of its activation free energy landscape and measurements of the appropriate geometrical parameters. This is done by focusing on the allosteric signaling pathways in different activation states, from inactive to active states via intermediate and pre-active ones, as well as the characterization of putative drug-binding pockets. The long-term objectives are to eventually be able to aid drug discovery efforts for this therapeutically important GPCR.

    Key findings and major conclusions: Conventional MD does not enable the sampling of the complete conformational landscape of receptor activation. Instead, enhanced sampling MD simulations are required to achieve this. Using metadynamics, the authors decipher the activation pathway of A1R, decode the allosteric networks and identify transient pockets. The protein energy networks computed throughout the inactive, intermediate active, pre-active and active conformational states unravel the extra and intracellular allosteric centers and the communication pathways that couple them, whereby the pathways are reinforced in the activated state. These conformations primarily differ in the dynamics of the ionic lock motif that couples TM3 to TM6 in the inactive conformation and reveal that G-proteins are required to fully stabilize the active conformation. Support for these findings comes from prior mutagenesis work on the A1R that identified key allosteric residues that in many cases map to identified communication nodes. Finally, the authors identified allosteric pockets throughout the A1R in four different conformational states that support prior experimental and MD studies on the mechanism of the positive allosteric modulator MIPS521 and which could be targeted for the design of modulators. Overall, these findings provide complementary support to a structure-based mechanism of activation and allosteric modulation of A1R, and extend the findings to incorporate dynamics across the full activation pathway.

    The perceived strengths and weaknesses: This preprint employs a combination of computational techniques to successfully reconstruct and analyze the conformational ensemble of the A1R activation. The metadynamics simulations supported the aim of the study, the results are clearly presented and the work is very well written. The authors could improve the discussion of how the protein energy network analysis could further advance rational design of specific modulators with a desired mode of action. The computational approach needs to be refined to be robust, with a focus on characterizing the convergence of the free energy landscapes. Overall, A1R is a good choice as the target for this study as there is existing structural and pharmacological data to support preliminary findings. Moreover, the framework presented herein could be adapted and scaled to other GPCRs with structural templates, which might enable comparison of allosteric pathways across families and classes.

    We thank the reviewers for contextualizing the findings of this study and for highlighting that our work provides complementary support to the mechanism of activation and allosteric modulation of A1R. We also thank the reviewers for their comments and suggestions, which had a great contribution to improve the quality and significance of the revised version of the preprint.

    RECOMMENDATIONS

    Revisions essential for endorsement:

    1. The paper could better demonstrate how the insights gained herein will or could lead to progress in the rational design of specific modulators with a desired effect. The authors should outline and discuss how they envision the modeling pipeline they have designed will be used towards this goal and tone-down or explain why "this information is essential to ease the design of allosteric modulators for A1R.". A recent study on FFAR1, where the authors targeted a specific dynamic pocket could be helpful in this respect (https://www.pnas.org/doi/full/10.1073/pnas.1811066116). Specifically, this might entail: How does specificity for a receptor correlate with pockets forming in a specific state? From this, how does one design an agonist vs. an antagonist vs. an inverse agonist? Does breaking a specific network select a function of the drug? How would another group follow up on this work, for example in a virtual screening campaign?

    We agree with the reviewers that a more comprehensive discussion of the points they mention is highly relevant to the study. Firstly, we have rephrased the last sentence of the abstract as "This information can be useful to ease the design of allosteric modulators for A1R" to ensure the significance of our results is not overstated. However, to address the reviewer's feedback more thoroughly, we conducted additional simulations with a positive allosteric modulator (MIPS51) and added an additional sub-section in the results, which includes three new figures (Figure 7-8 and S12):

    "ADO and MIPS51 PAM have a significant impact on the energy networks. In order to establish a connection between the energy networks and the mode of action of allosteric modulators, we focus on exploring the effect of MIPS521 positive allosteric modulator (PAM) and ADO agonist as a proof of concept. Experimental assays and Gaussian accelerated MD determined that MIPS521 PAM increases the binding affinity of ADO in the orthosteric site.[19] Thus, PB and PD must be allosterically coupled. Among MIPS521 PAM pocket residues, only L2426.43, L2456.46, S2466.47 and G2797.44 were experimentally found to affect the PAM cooperativity (Figure S11). Interestingly, the PEN obtained in presence of ADO captures these key residues along activation, including TM6 (L2426.43 and L2456.46) in the intermediate, L2426.43 and S2466.47 in the pre-active and L2426.43 and TM7(G2797.44) allosteric residues in the fully-active ensemble (Figure 7). Indeed, G2797.44 becomes a key node in the PEN of the fully-active ensemble. This evidence suggests that although both PD and PB are open in all conformational states, their energy coupling is particularly stronger during the receptor activation.

    This prompted us to investigate whether the binding of ADO and the MIPS521 PAM can affect the allosteric communication between PB and PD sites. To that end, we performed cMD of the heterotrimeric Gi2 protein ADO-A1R-Gi2 complex in presence of the PAM (PAM-ADO-A1R-Gi2 complex) and in absence of adenosine (A1R-Gi2 complex) in order to compute their conformational landscape and energy networks following the same protocol for the ADO-A1R-Gi2 complex (Figure 8 and S12). The analysis of the PEN of A1R-Gi2 complex reveals that in the absence ADO, the receptor displays a reduced allosteric communication between PB and functional regions of the receptor, such as the extracellular allosteric center, TM6 and PD allosteric site. As expected, the presence of ADO restores the allosteric coupling between PB and TM6, which could explain the increase of receptor activity associated with agonist binding. Additionally, our analysis of the PAM-ADO-A1R-Gi2 complex shows that the PAM reinforces the TM7-ECL3-ECL2 allosteric pathway that couple PD with PB, and ECL2 now communicates to the intracellular region through TM5 (Figure 8). Notably, a recently published study reported that the orhosteric pocket contracts after ADO binding, as demonstrated by shortened distances of the so-called vestibular lid (defined as the sum of length of the triangle perimeters formed by E17045.51-Y2717.36-E17245.53 interacting residues) and the E17245.53-K26567 salt bridge.[48] Remarkably, the TM7-ECL3-ECL2 enhanced pathway by PAM effect contains the vestibular lid and the E17245.53-K26567 salt bridge residues (Figure 8). This suggest that PAM promotes the contraction of PB, leading to the stabilization of the ADO-bound state. Thus, the enhanced energy coupling between PB and PD may be responsible for the increase in the binding affinity of ADO in presence of the PAM, as observed experimentally.[19] This data indicates that allosteric modulators are able to enhance and redistribute the energy networks, which is likely attributed for their effects on the receptor activity."

    The new insights gained from our additional simulations have significantly enriched the discussion on how the protein energy network analysis can contribute to the rational design of specific modulators with desired modes of action. In light of these finding, the last paragraph of the discussion has been rephrased as:

    "As a proof of concept, we focus on the PD, which corresponds to the binding site of MIPS52, a positive allosteric modulator (PAM) that increases Adenosine binding affinity in the orthosteric site (PB). Although PD is open in all conformational states the communication between PB and PD is enhanced along activation capturing the allosteric residues that were found to affect its PAM from the intermediate to the fully active. Based on this observation, we hypothesize that drugs that bind pockets and interact with PEN residues, which progresses towards regions of the receptor where function can be altered, may potentially affect the receptor activity through allosteric effects. Additionally, the pocket where the drug binds must be open at least in the conformation state that is targeted. As a practical aspect, virtual screening campaigns could use this information during the design procedure by selecting drug candidates that perform stronger interactions with the PEN contained in the pockets.

    To further support this hypothesis, we explored the allosteric effects of ADO and MIPS52 PAM on the PEN. Interestingly, we observed that ADO is crucial for the formation of the extracellular center and its connection with TM6 pathway. Furthermore, MIPS52 PAM reinforces the pathway that connects PB and PD pockets and redistribute other connections. These alterations in the PEN can be related with their mode of action. ADO may increase the activity of the receptor through its communication with TM6 and the PAM may increase ADO binding affinity though stronger energy coupling between PD and PB pockets. These findings imply that the mode of action of allosteric drugs could be predicted depending on how they redistribute the PEN."

    Accordingly, the last paragraph of the conclusions has been rephrased as:

    "As a proof of concept, Adenosine and a previously experimentally determined positive allosteric modulator were found to enhance and redistribute the energy networks of the receptor in a manner that is consistent with their respective functions. The prediction of drug effects depending on how they redistribute the protein energy networks presents a promising avenue for drug discovery. All these system-specific structural dynamics understanding provides useful information to advance the design of A1R allosteric modulators on the basis of structure-based drug design. This computational approach can be also transferable to other GPCRs and related receptors, which is of interest for the design of novel allosteric drugs."

    1. Free energy calculations:

    a. A proof of convergence of the free energy calculations is missing. The authors argue that obtaining landscapes that do not change over time is proof of convergence, but this is incorrect in well-tempered metadynamics. The fact that the heights of the Gaussians decrease over time guarantees that the landscape will be stable over time, and the way to check convergence is to show that the collective variables become diffusive after convergence. In addition, to validate that the choice of collective variables (CV) is actually appropriate, they should check that CVs that were not biased are also diffusive. This would be best studied by looking at the behavior of microswitches that were not considered, such as ones describing the PIF motif, the NPxxY motif, the ligand binding pose, etc.

    The goal of the well-tempered approach [Phys. Rev. Lett. 2008, 100, 020603] is to improve the convergence of the energy landscapes. This is achieved by gradually decreasing the height of gaussians over simulation. In this fashion, the height of the Gaussians is proportional to a decaying exponential function of the potential deposited in the currently visited point of the CV space. This technique has the added benefit of constraining reconstruction to the region of interest, reducing the risk of irreversible movement towards physically irrelevant regions of the CV space. As noted in the plumed tutorial (https://www.plumed.org/doc-v2.7/user-doc/html/master-_i_s_d_d-2.html), the fact that the Gaussian height is gradually decreasing should not be used as a measure of convergence. Rather, convergence can be assessed by monitoring the energy differences between chosen regions of the energy landscape over time (e.g. the inactive, intermediate and pre-active local energy minima used in our work). If this energy differences do not change significantly as a function of time, this can be taken as an indicator of convergence. We also would like to emphasize that our aim is to recover the major conformational states involved in the pathway of receptor activation rather than the study of subtle energy barriers and relative stability differences of the energy minima upon system perturbations. This objective has been made clear in the text.

    We agree with the reviewers' comment that our measure of the convergence could be strengthened by additional analysis that verify the computed conformational pathway of receptor activation.

    As suggested by the reviewers, we have plotted the CV1 values over time to asses convergence of our simulations (Figure S4). However, it should be noted that in this study, we have employed the walkers approach [J. Phys. Chem. B 2006, 110, 3533], which utilizes 10 replicas (walkers) to parallelize free energy reconstruction. Each walker simultaneously reconstructs the energy landscape by reading the Gaussian potentials deposited by the other walkers. Consequently, the correct time to stop the metadynamics simulation using the walkers approach becomes more problematic. To facilitate efficient sampling of the CV space and achieve convergence more rapidly, we utilized a sampling strategy that involved starting the simulation with walkers that spanned the entire CV space of interest. In this context, the fact that the walkers do not become trapped in the initial CV space and are able explore and cross into regions occupied by other walkers may serve as a useful indicator of convergence. This assessment of the convergence has been implemented before for Tryptophan Synthase, in which the resulting energy landscapes were consistent with experimental data. [J. Am. Chem. Soc. 2019, 141, 13049] and [ACS Catal. 2021, 11, 13733]

    We have added the following paragraph in the convergence section including Figure S4:

    "We have also assessed convergence by analyzing the CV1 values over simulation time. Figure S4A shows that during the first 100ns, walkers primary oscillates around their initial CV1 values. Subsequently, at around 200 ns walkers exhibit a higher frequency of crossing into regions occupied by other walkers. This is further supported by the exploration of W1 and W10, as shown in Figure S4B. These two walkers initially start the landscape reconstruction at the opposite extremes of the CV space. At 120 ns, they are able to escape from their respective basins and approach each other, sampling similar CV values (at approximately 240 ns). At this point of the simulation, only these two walkers have covered the entire conformational space of activation. Subsequently, they tend to return to previously sampled CV space. The observation that walkers do not become trapped in their initial CVs region, but instead explore and cross into other regions suggests that our sampling strategy, which involved starting the simulations with walkers that spanned the entire CV space of interest, has facilitated the exploration of the relevant conformational space. Although we cannot guarantee full convergence of the free energy landscape under these conditions, we successfully reconstructed the major conformational states of the receptor activation at 250 ns."

    Accordingly, in the results section we have replaced "After 250 ns of accumulated time the FEL was considered to be converged (see Figure S3 and S4)" by "After 250 ns of accumulated time, we successfully reconstructed the major conformational states of the FEL (see convergence assessment in Figure S3 and S4)."

    To further verify the accuracy of the collected conformational landscape, we have conducted additional analysis, which include the following:

    "As a complementary analysis, we conducted the reweighting of the metadynamics simulations[28] to determine the free energy as a function of previously identified A1R micro-switches (ionic-lock, PIF motif, water-lock and toggle switch). The fact that we capture the distinct energy barriers associated with unbiased micro-switches highlights the accuracy of the metadynamics simulations in reproducing the pathway of activation and provides useful information to guide the selection of collective variables for future GPCR landscape calculations (Figure S5, S6 and S7)."

    b. The authors should characterize the uncertainties/statistical errors on the measured free energy profiles to better evaluate the significance of change (e.g. for inspiration: https://www.plumed.org/doc-v2.7/user-doc/html/masterclass-21-2.html).

    In response to the reviewers' comment, we have included a new sub-section in materials and methods together with an additional figure in the Supplementary Information (Figure S13), as follows:

    "Error: We estimated the error on the 2D free energy landscape of the first collective variable (CV1), which is the TM3-TM6 intracellular ends distance (Figure S13) using the block averaging technique, as described in the PLUMED tutorial on calculating error bars (https://www.plumed.org/doc-v2.8/user-doc/html/lugano-4.html). We calculated the weights using the metadynamics bias potential obtained at the end of the simulation, and assuming a constant bias during the entire course of the simulation.[28] Specifically, we calculate the error using blocks of histograms of 25 ns each, covering the entire 250 ns simulation time."

    c. In the cMD trajectories, a large part of phase space is sampled, which does not appear consistent with what one would expect based on the free energy landscapes. For instance, it does not seem reasonable to cover an almost complete conformational transition in 500ns when the barrier of the system is on the order of 5-8kcal/mol. The definition of CVs may have led to an overestimation of the free energy barrier. Hence an independent validation of the free energy barrier height is needed, by e.g. changing the CV definition.

    We agree with the reviewers' that the cMD simulations cover a large part of the phase space. This is in part due to running simulations staring from both the inactive and active structures. For the latest, we removed the G-proteins from the receptor. This situation increases the flexibility of the receptor and induces a population shift respect to the starting point. However, as expected, we did not observe any complete transitions in our cMD simulations either staring from the inactive or active structures (see Figure 1C and S2).

    Regarding the activation energy barrier, we would like to clarify that the aim of our study is not to compare subtle differences in energy barriers after system perturbations or compare them with experimental data. As far as we know, there is not NMR data available that confirms the exact time-scale of activation for A1R receptor, suggesting that A1R could be highly flexible. Notably, we report a rather low activation energy barrier of approximately 4 kcal/mol derived from the metadynamics simulations of A1R in the presence of adenosine. This is consistent with other computational studies of A1R where a complete transition from active to inactive [PNAS 2022, 119, E2203702119] and from inactive to pre-active [Nature 2021 597, 571] states is sampled in the course of nanosecond-scale Gaussian accelerated MD simulations. In addition, similar activation barrier values have been computed for the b2 adrenergic GPCR in the presence of adrenaline using different collective variables, as reported in [PLoS Comput. Biol. 2011, 7, e1002193] and [eLife 2021, 10, e60715]."

    After careful consideration, we found relevant to validate our path of conformations by reweighing of our free energy into other collective variables. As previously stated, we have reweighted the original free energy into the ionic-lock, PIF motif, water-lock and toggle micro-switches (refer to the last paragraph of Essential revision 2.1 and Figure S5 and S6). Our analysis reveals that our original CV2 (TM6 torsion), the PIF motif, and the toggle display modest energy barriers, while our original CV1 (TM3-TM6 distance) presents a rather high energy barrier. Moreover, the ionic-lock and the water-lock exhibit the highest energy barriers. Thus, if we had chosen to use our initial CV1 and the PIF motif, we would have obtained a similar energy barrier, whereas choosing our initial CV1 and the water-lock would have resulted in a higher energy barrier, as predicted by our reweighting calculations (see Figure S7).

    As mentioned in the manuscript, this data highlights the ability of our simulations to reproduce the activation pathway and provides interesting insights that can guide the selection of collective variables for future GPCR landscape calculations.

    1. Configurations extracted from both conventional MD and wt-metadynamics are mixed in the analyses of the allosteric networks and the pockets. A more accurate way to integrate these datasets would be to modulate the weights of the configurations by their statistical weights, which can be retrieved from the metadynamics simulations.

    We thank the reviewers for this suggestion and we will consider to use configurations by their statistical weights in future work. For this case, we aimed to include as much as configurations as possible from each conformational state. We then included all configurations from the inactive, intermediate and pre-active derived from the metadynamics and for the fully-active we applied a stride value in order to collect a similar number of structures.

    1. Related to Figure S6, it is essential to compare the dynamics for all of the key class A activation motifs including the Na binding site, PIF motif, and NPxxY.

    Based on the reviewers' comments, we have generated histograms for the relevant micro-switches corresponding to the inactive, intermediate, and pre-active states (see Figure S10). This analysis provides further support to the activation pathway derived from the metadynamics simulations. Notably, the population distributions of these micro-switches in the inactive, intermediate, and pre-active states exhibit a correlated progression that mirrors the receptor's activation pathway.

    1. Please provide clarification on why 500 ns was chosen as the time-scale of the MD simulations and inclusion of the time course for the three independent MD simulations for each of the key structural features (e.g. TM6 torsion angles and TM3-TM6 distances).

    To ensure adequate initial sampling of receptor activation, we performed three independent MD replicas of 500 ns each, starting from both the inactive and active structures. The purpose of these MD simulations is to provide an initial sampling of the receptor instead of describing the complete activation pathway. Based on this initial sampling, we selected a path of 10 conformations as starting points for the walker metadynamics simulations. We have added this information to the text for clarification.

    "For each starting point we computed 3 replicas of 500ns, which is a reasonable simulation time to provide an initial sampling of the receptor activation."

    1. The validation of the results in the form of previously published mutagenesis results does not appear completely convincing. Large parts of the protein are included in the allosteric network, making it likely that mutations in some of these residues will have an effect if mutated. In addition, the fact that mutations in ECL2 and ECL3 affect allostery is expected and does not constitute a good validation of the results. If no other results are included, we recommend that the language be toned down so as not to overstate the significance of the results.

    Following the reviewers' comment, we have removed the following sentences from the results:"Among the PEN residues, we successfully captured most of the allosteric residues previously identified by mutagenesis studies, which highlights the reliability of the allosteric networks computed" and "The high predictive power of the PEN to identify allosteric residues highlights the reliability of the characterized allosteric pathways"

    1. What is the justification for using an energy-based scoring for network analysis, given that a conventionally correlation-based approach has been used successfully in the field? The concern with an energy-based approach is that the interaction energy calculations do not consider the dielectric effect, i.e., if water molecules interfere with two interacting residues. Since the dynamic network is one of the critical aspects of this study, we believe the authors need to explore other tools such as the one implemented in VMD (https://www.ks.uiuc.edu/Research/vmd/plugins/networkview/) and compare the results.

    We decided to perform protein energy networks analysis because our aim was to investigate how the networks change in different conformational states along activation. We selected this approach because energy networks can provide a more detailed insight into how communication within the protein changes during activation, as compared to cross-correlation networks, which are more suited to characterizing communication through correlated motions in the global ensemble. In order to compare both protein energy networks and correlation networks we performed the cross-correlation analysis in the global ensemble. Note that both approaches yield some similarities and provides complementary information.

    We also would like to thank the reviewers for raising concerns about the methodology employed in our work.We acknowledge that gRINN, which we used to generate the pairwise residue mean interaction energy matrix, does not include water molecules and ligands during the matrix generation process. As a result, we are unable to capture communication pathways that involve water-mediated connections or interactions between ligands and residues. For example, in the pre-active ensemble, Y2005.58 and Y2887.53 from the NPxxY motif are connected by a hydrogen bond facilitated by a bridging water molecule (the son-called water-lock). However, such communication is not captured in our analysis due to the absence of water molecules (Figure 3). We highlight this major limitation in the text. Another example can be observed in the simulation of the PAM-ADO-A1R-Gi2 system, where communication between L242 and S246, two residues involved in the Positive Allosteric Modulator (PAM) binding site, is missing in the PEN (Figure 8). Since these residues must be connected through the PAM, our methodology cannot detect their communication. A promising tool to considered in future studies is webPSN v2.0 [Nucleic Acids Res. 2020, 48, W95], a protein structure network analysis that includes nucleic acids and more than 30,000 biologically relevant molecules and ions, which is highly advantageous to study the effect of drugs on the protein communication.

    1. Provide generic residue numbers such as GPCRdb or Ballesteros Weinstein numbering for all mentioned residues in text and figures, as is standard for structural papers.

    As the reviewers suggested, we have renumbered all residues mentioned in the text and figures according to the Ballesteros Weinstein numbering scheme.

    Additional suggestions for the authors to consider:

    1. For the PEN analysis it would be useful to digest these communication networks with respect to the established structural activation motifs of class A GPCRs (Na binding site, PIF, and NPxxY) that are present at the A1R.

    In order to make the PEN analysis more digestive, we have revised the second paragraph of the "Energy Networks captures the dynamic allosteric pathways along A1R activation" results section. Specifically, we have highlighted the most relevant micro-switches captured in the PEN, with a particular focus on the ionic and water-lock switches, which are the most prominent for the protein communication.

    1. It is unclear why the authors chose two largely correlated CVs (See comment 2c). In addition, the choice of CV is likely contributing to the distortion of S6, as displayed in Figure 1E. It has been shown that choosing a different CV set that describes the motion between states in a more distributed way is more likely to lead to a converged conformational ensemble. We suggest repeating the metadynamics simulations with a more distributed CV set that encompasses all of the microswitches in the receptor.

    Regarding the concern raised by the reviewers about the distortion observed in the TM6 end, we want to clarify that it is not attributed to the selection of collective variables (CVs) since it is already explored in the initial conventional MD simulations (see Figure S2B, W5 structure). We selected two CVs that may seem largely correlated, but they actually describe different aspects of the TM6 inward-to-outward transition. The first CV, which measures the distance between the center of mass of TM3 and TM6, is more related to the dynamics of the ionic lock in the intracellular region. On the other hand, the second CV (TM6 torsion) is related to the forces sensed by the upper region of TM6, including the dynamics of the W2476.48 toggle switch. Therefore, we believe that the combination of these two CVs provides a comprehensive description of the TM6 transition.

    It is also worth mentioning that a more distributed set of CVs may be beneficial to better reproduce the activation energy barriers of the receptor. In fact, as shown in the reweighting calculations of the metadynamic bias potential (see Figures S6 and S7), using the TM3-TM6 COM end distances and the Y2005.58-Y2887.53 distance (water-lock) as CVs appears to be a good choice for this purpose.

    1. To support the vision on how the analysis of activation pathway, energy networks and transient pockets could be used "to ease the design of allosteric modulators for A1R" (last sentence of the abstract), it might be necessary to show that the combination of these methods can indeed be predictive for the binding and effect of known ligands. This might provide a first step towards establishing that molecules that bind to pockets "near allosteric networks" is a promising avenue for drug discovery.

    This highly relevant point has been addressed in Essential revisions 1.

    1. The specific TM3-TM6 residues should be specified in figures and text. Commonly used TM3-TM6 comparisons include the measured maximum distance between 2x46 to 6x37, which could be used here also (e.g. see https://docs.gpcrdb.org/structures.html#structure-descriptors).

    The specific residues of TM3 and TM6 that were used in the analysis have been clearly specified in both the Materials and Methods section and the figure captions of Figure 1 and 4.

    1. Even though the "A1R in complex with PSB36 (PDB 5N2S)" is an inactive structure, PSB36 is an agonist. Hence, the authors should consider using the DU172 antagonist-bound structure for comparison (PPDB 5UEN)

    According to literature, PSB36 is selective antagonist for A1R. In fact, experimental data showed that PSB36 exhibit low inverse agonist activity [Chem. Med. Chem. 2006, 1, 891]. Although PDB 5N2S and PPDB 5UEN structures are almost identical, the bulkier DU172 ligand causes a displacement of TM2 in the extracellular region. Therefore, we chose to use the 5N2S structure in our study. However, we will consider using the 5UEN structure in future studies.

    1. How does adenosine and MIPS521 binding impact the different conformational states and PEN.

    This highly relevant point has been addressed in Essential Revisions 1. For a more detailed analysis, refer to Figure 8 and S12, which shows the impact of MIPS521 on the PEN and the conformational landscape, respectively.

    1. It would be interesting to note how the findings from this study compare/contrast to a very recently published report by Li et al, PNAS, 2022 "The full activation mechanism of the adenosine A1 receptor revealed by GaMD and Su-GaMD simulations". Similarly with regards to the determination of allosteric binding pockets in this recent publication: "The pocketome of G-protein-coupled receptors reveals previously untargeted allosteric sites" (https://doi.org/10.1038/s41467-022-29609-6)

    According to the reviewers' suggestion, we have compared our findings to those of a recently published study [PNAS 2022, 119, E2203702119], which explored the full activation mechanism of A1R using both su-GaMD and GaMD.

    This published work serves to further confirm the combined activation mechanism that we observed for A1R in our study, which entails the formation of a pre-active state in the presence of Adenosine and the stabilization of a fully-active state in the presence of both Adenosine and G-proteins. Moreover, their study reports the pre-activation of the receptor from the inactive state within 150 ns of GaMD, indicating a rather low activation energy barrier of the receptor in presence of Adenosine. This is consistent with the approximately 4 kcal/mol activation energy barrier we calculated for A1R-ADO in our 250 ns metadynamic simulations.

    We would also like to highlight an additional noteworthy point we have included in the results as: "Notably, a recently published study reported that the orthosteric pocket contracts after ADO binding, as demonstrated by shortened distances of the so-called vestibular lid (defined as the sum of length of the triangle perimeters formed by E17045.51-Y2717.36-E17245.53 interacting residues) and the E17245.53-K26567 salt bridge.[48] Remarkably, the TM7-ECL3-ECL2 enhanced pathway by PAM effect contains the vestibular lid and the E17245.53-K26567 salt bridge residues (Figure 8). This suggest that PAM promotes the contraction of PB, leading to the stabilization of the ADO-bound state. Thus, the enhanced energy coupling between PB and PD may be responsible for the increase in the binding affinity of ADO in presence of the PAM, as observed experimentally.[19]"

    1. A major advantage of allosteric drugs is the potential to achieve higher selectivity. Expansion of this study to include other adenosine receptor subtypes or linking to other types of molecular pharmacology (e.g. biased signalling, subtype selectivity, etc.) would be a major benefit to the field.

    We are grateful to the reviewers for recognizing the potential impact of our work in various applications beyond our initial scope. We will consider to incorporate their valuable suggestions in our future research endeavors.

    1. Consider including an explanation of the physiological and pharmacological relevance of A1AR in the introduction.

    According to this suggestion, we have incorporated a new sentence in the introduction section as follows:

    "The adenosine A1 receptor (A1R) is a member of the class A G protein-coupled receptor (GPCR) family that preferentially couples with Gi/o proteins. It is widely distributed in multiple organs mediating a variety of physiological processes, including those in the brain and the heart. Thus, A1R has significant therapeutic potential in the treatment of numerous diseases and disorders.[18]"

    1. Even if not entirely necessary for the results, it would be more consistent if the study would include metadynamics of the G-protein bound state.

    Performing a metadynamics calculation of the G-protein bound state is challenging as it requires careful consideration of the G-protein binding process. As a reference work, Giulio Mattedi et al. successfully implemented this calculation for the glucagon receptor [Proc. Natl. Acad. Sci USA, 2020, 117, 15414]. However, in our study the effect of the G-proteins in the activation landscape is a minor remark. Our study places a greater emphasis on sampling the most stable conformations associated with the fully-active conformational state to compute the protein energy networks.

    1. Methods: "In other words, once the free energy surface does not change significantly during a relatively long period of time in the last part of the simulation". What is "relatively long period of time" and "change significantly". The convergence, should be stated as a quantitative description of the observed energy differences.

    We have addressed this technical issue in Essential revisions 2. It is worth noting that "the relatively long period of time" required for convergence may vary depending on the specific system under study. Nonetheless, we believe that observing a stable energy surface over the course of 50-100 ns, while the system explores different relevant regions of the CV space, provides a good criterion for convergence."

    1. The authors should strongly consider making their analysis code and simulation data publicly available (e.g. on GitHub or Zenodo) so that others can replicate and build upon this work.

    We thank the reviewers for this suggestion. We will make all the output files generated from the get Residue Interaction eNergies and Networks (gRINN) calculation available to the public.By doing so, users will be able to visualize the results in the gRINN visual interface and perform customized network analysis using the pairwise residue mean interaction energy matrices. Additionally, we will provide all Pymol sessions that include the protein energy networks as well as the Isosurface representation of the frequency maps of the transient pockets. We believe that these materials will provide better visualization compared to the current figures presented in the manuscript and supporting information, which will be helpful in guiding future structure-based drug design campaigns.

    (This is a response to peer review conducted by Biophysics Colab on version 3 of this preprint.)

  9. Biophysics Colab

    Consolidated peer review report (28 November 2022)

    GENERAL ASSESSMENT

    The objectives of the study:

    This paper aims to characterize the dynamics that drive allostery of the adenosine A1 receptor (A1R) via computational analysis of its activation free energy landscape and measurements of the appropriate geometrical parameters. This is done by focusing on the allosteric signaling pathways in different activation states, from inactive to active states via intermediate and pre-active ones, as well as the characterization of putative drug-binding pockets. The long-term objectives are to eventually be able to aid drug discovery efforts for this therapeutically important GPCR.

    Key findings and major conclusions:

    Conventional MD does not enable the sampling of the complete conformational landscape of receptor activation. Instead, enhanced sampling MD simulations are required to achieve this. Using metadynamics, the authors decipher the activation pathway of A1R, decode the allosteric networks and identify transient pockets. The protein energy networks computed throughout the inactive, intermediate active, pre-active and active conformational states unravel the extra and intracellular allosteric centers and the communication pathways that couple them, whereby the pathways are reinforced in the activated state. These conformations primarily differ in the dynamics of the ionic lock motif that couples TM3 to TM6 in the inactive conformation and reveal that G-proteins are required to fully stabilize the active conformation. Support for these findings comes from prior mutagenesis work on the A1R that identified key allosteric residues that in many cases map to identified communication nodes. Finally, the authors identified allosteric pockets throughout the A1R in four different conformational states that support prior experimental and MD studies on the mechanism of the positive allosteric modulator MIPS521 and which could be targeted for the design of modulators. Overall, these findings provide complementary support to a structure-based mechanism of activation and allosteric modulation of A1R, and extend the findings to incorporate dynamics across the full activation pathway.

    The perceived strengths and weaknesses:

    This preprint employs a combination of computational techniques to successfully reconstruct and analyze the conformational ensemble of the A1R activation. The metadynamics simulations supported the aim of the study, the results are clearly presented and the work is very well written. The authors could improve the discussion of how the protein energy network analysis could further advance rational design of specific modulators with a desired mode of action. The computational approach needs to be refined to be robust, with a focus on characterizing the convergence of the free energy landscapes. Overall, A1R is a good choice as the target for this study as there is existing structural and pharmacological data to support preliminary findings. Moreover, the framework presented herein could be adapted and scaled to other GPCRs with structural templates, which might enable comparison of allosteric pathways across families and classes.

    RECOMMENDATIONS

    Revisions essential for endorsement:

    1. The paper could better demonstrate how the insights gained herein will or could lead to progress in the rational design of specific modulators with a desired effect. The authors should outline and discuss how they envision the modeling pipeline they have designed will be used towards this goal and tone-down or explain why “this information is essential to ease the design of allosteric modulators for A1R.”. A recent study on FFAR1, where the authors targeted a specific dynamic pocket could be helpful in this respect (https://www.pnas.org/doi/full/10.1073/pnas.1811066116). Specifically this might entail: How does specificity for a receptor correlate with pockets forming in a specific state? From this, how does one design an agonist vs. an antagonist vs. an inverse agonist? Does breaking a specific network select a function of the drug? How would another group follow up on this work, for example in a virtual screening campaign?

    2. Free energy calculations:

    a. A proof of convergence of the free energy calculations is missing. The authors argue that obtaining landscapes that do not change over time is proof of convergence, but this is incorrect in well-tempered metadynamics. The fact that the heights of the Gaussians decrease over time guarantees that the landscape will be stable over time, and the way to check convergence is to show that the collective variables become diffusive after convergence. In addition, to validate that the choice of collective variables (CV) is actually appropriate, they should check that CVs that were not biased are also diffusive. This would be best studied by looking at the behavior of microswitches that were not considered, such as ones describing the PIF motif, the NPxxY motif, the ligand binding pose, etc.

    b. The authors should characterize the uncertainties/statistical errors on the measured free energy profiles to better evaluate the significance of change (e.g. for inspiration: https://www.plumed.org/doc-v2.7/user-doc/html/masterclass-21-2.html).

    c. In the cMD trajectories, a large part of phase space is sampled, which does not appear consistent with what one would expect based on the free energy landscapes. For instance, it does not seem reasonable to cover an almost complete conformational transition in 500ns when the barrier of the system is on the order of 5-8kcal/mol. The definition of CVs may have led to an overestimation of the free energy barrier. Hence an independent validation of the free energy barrier height is needed, by e.g. changing the CV definition.

    1. Configurations extracted from both conventional MD and wt-metadynamics are mixed in the analyses of the allosteric networks and the pockets. A more accurate way to integrate these datasets would be to modulate the weights of the configurations by their statistical weights, which can be retrieved from the metadynamics simulations.

    2. Related to Figure S6, it is essential to compare the dynamics for all of the key class A activation motifs including the Na binding site, PIF motif, and NPxxY.

    3. Please provide clarification on why 500 ns was chosen as the time-scale of the MD simulations and inclusion of the time course for the three independent MD simulations for each of the key structural features (e.g. TM6 torsion angles and TM3-TM6 distances).

    4. The validation of the results in the form of previously published mutagenesis results does not appear completely convincing. Large parts of the protein are included in the allosteric network, making it likely that mutations in some of these residues will have an effect if mutated. In addition, the fact that mutations in ECL2 and ECL3 affect allostery is expected and does not constitute a good validation of the results. If no other results are included, we recommend that the language be toned down so as not to overstate the significance of the results.

    5. What is the justification for using an energy-based scoring for network analysis, given that a conventionally correlation-based approach has been used successfully in the field? The concern with an energy-based approach is that the interaction energy calculations do not consider the dielectric effect, i.e., if water molecules interfere with two interacting residues. Since the dynamic network is one of the critical aspects of this study, we believe the authors need to explore other tools such as the one implemented in VMD (https://www.ks.uiuc.edu/Research/vmd/plugins/networkview/) and compare the results.

    6. Provide generic residue numbers such as GPCRdb or Ballesteros Weinstein numbering for all mentioned residues in text and figures, as is standard for structural papers.

    Additional suggestions for the authors to consider:

    1. For the PEN analysis it would be useful to digest these communication networks with respect to the established structural activation motifs of class A GPCRs (Na binding site, PIF, and NPxxY) that are present at the A1R.

    2. It is unclear why the authors chose two largely correlated CVs (See comment 2c). In addition, the choice of CV is likely contributing to the distortion of S6, as displayed in Figure 1E. It has been shown that choosing a different CV set that describes the motion between states in a more distributed way is more likely to lead to a converged conformational ensemble. We suggest repeating the metadynamics simulations with a more distributed CV set that encompasses all of the microswitches in the receptor.

    3. To support the vision on how the analysis of activation pathway, energy networks and transient pockets could be used “to ease the design of allosteric modulators for A1R” (last sentence of the abstract), it might be necessary to show that the combination of these methods can indeed be predictive for the binding and effect of known ligands. This might provide a first step towards establishing that molecules that bind to pockets “near allosteric networks” is a promising avenue for drug discovery.

    4. The specific TM3-TM6 residues should be specified in figures and text. Commonly used TM3-TM6 comparisons include the measured maximum distance between 2x46 to 6x37, which could be used here also (e.g. see https://docs.gpcrdb.org/structures.html#structure-descriptors).

    5. Even though the "A1R in complex with PSB36 (PDB 5N2S)" is an inactive structure, PSB36 is an agonist. Hence, the authors should consider using the DU172 antagonist-bound structure for comparison (PPDB 5UEN)

    6. How does adenosine and MIPS521 binding impact the different conformational states and PEN.

    7. It would be interesting to note how the findings from this study compare/contrast to a very recently published report by Li et al, PNAS, 2022 “The full activation mechanism of the adenosine A1 receptor revealed by GaMD and Su-GaMD simulations”. Similarly with regards to the determination of allosteric binding pockets in this recent publication: “The pocketome of G-protein-coupled receptors reveals previously untargeted allosteric sites” (https://doi.org/10.1038/s41467-022-29609-6)

    8. A major advantage of allosteric drugs is the potential to achieve higher selectivity. Expansion of this study to include other adenosine receptor subtypes or linking to other types of molecular pharmacology (e.g. biased signalling, subtype selectivity, etc.) would be a major benefit to the field.

    9. Consider including an explanation of the physiological and pharmacological relevance of A1AR in the introduction.

    10. Even if not entirely necessary for the results, it would be more consistent if the study would include metadynamics of the G-protein bound state.

    11. Methods: "In other words, once the free energy surface does not change significantly during a relatively long period of time in the last part of the simulation". What is “relatively long period of time” and “change significantly”. The convergence, should be stated as a quantitative description of the observed energy differences.

    12. The authors should strongly consider making their analysis code and simulation data publicly available (e.g. on GitHub or Zenodo) so that others can replicate and build upon this work

    REVIEWING TEAM

    Reviewed by:

    Antonios Kolocouris, Professor, Department of Medicinal Chemistry Faculty of Pharmacy National and Kapodistrian University of Athens, Greece:

    CADD and computational biophysics on adenosine receptors

    David Thal, Senior Research Officer, Monash University, Australia:

    structural biology and molecular pharmacology of allosteric mechanisms underlying Class A GPCRs

    SciLifeLab Journal Club, Stockholm, Sweden (see Appendix for members)

    Curated by:

    Alexander S. Hauser, Associate Professor, University of Copenhagen, Denmark

    APPENDIX

    SciLifeLab Journal Club:

    Feedback was generated in a meeting of the journal club involving:

    Lucie Delemotte (Journal Club oversight), Associate Professor of Biophysics, KTH Royal Institute of Technology, Sweden: modeling and enhanced sampling of GPCRs and other membrane proteins.

    Olivia Andén, PhD student, Stockholm University: cryo-EM and functional characterization of membrane proteins.

    Cathrine Bergh, PhD student, KTH Royal Institute of Technology: enhanced sampling simulations of membrane proteins.

    Koushik Choudhury, PhD student, KTH Royal Institute of Technology: membrane protein modeling, enhanced sampling simulations.

    John Cowgill, postdoctoral scholar, Stockholm University: cryo-EM and simulations of membrane proteins.

    Chen Fan, postdoctoral scholar, Stockholm University: cryo-EM and simulations of membrane proteins.

    Nandan Haloi, postdoctoral scholar, KTH Royal Institute of Technology: membrane protein modeling, free energy calculations, structure refinement in cryo-EM maps.

    Rebecca J Howard, researcher, Stockholm University: membrane protein structure-function, allosteric modulation.

    Marie Lycksell, PhD student, Stockholm University: structure and simulations of membrane proteins.

    Antoni Marciniak, PhD student, KTH Royal Institute of Technology: enhanced sampling simulations of GPCRs and other membrane proteins.

    Darko Mitrovic, PhD student, KTH Royal Institute of Technology: membrane protein modeling, enhanced sampling, machine learning.

    Alex Payne, PhD student, Memorial Sloan Kettering Center for Cancer Research: membrane protein modeling, cryo-EM structure determination, drug discovery.

    Urška Rovšnik, PhD student, Stockholm University: cryo-EM and functional characterization of membrane proteins.

    Akshay Sridhar, postdoctoral scholar, KTH Royal Institute of Technology: membrane protein modeling, enhanced sampling simulations.

    Amanda Dyrholm Stange, PhD student, Aarhus University: membrane protein modeling, enhanced sampling simulations.

    (This consolidated report is a result of peer review conducted by Biophysics Colab on version 3 of this preprint. Minor corrections and presentational issues have been omitted for brevity.)

  10. Version published to 10.1101/2022.08.08.502889v5 on bioRxiv
  11. Version published to 10.1101/2022.08.08.502889v4 on bioRxiv
  12. Version published to 10.1101/2022.08.08.502889v3 on bioRxiv
  13. Version published to 10.1101/2022.08.08.502889v2 on bioRxiv
  14. Version published to 10.1101/2022.08.08.502889v1 on bioRxiv