G protein-coupled receptors (GPCRs) are membrane proteins which activate cellular signaling in response to extracellular stimuli. This process is controlled by extracellular ligands such as hormones and neurotransmitters, and the binding of these increases the probability of activating intracellular partners. GPCRs are vital in many physiological processes and constitute the most common class of drug targets (Hauser et al., 2017).

Much of the current understanding of GPCR signaling at the molecular level can be attributed to the progress in GPCR structure determination during the last decade (Cherezov et al., 2007; Hanson et al., 2008; Masureel et al., 2018; Rasmussen et al., 2011a; Ring et al., 2013; Wacker et al., 2010). GPCRs interconvert between inactive (R) and active (R*) states, which control G protein binding to a conserved intracellular domain via conformational rearrangements among the seven transmembrane (TM) helices (Figure 1a; Manglik and Kruse, 2017). In the absence of a bound agonist, this process is called basal activity. Ligands can bind to the orthosteric site in the receptor’s extracellular domain and thereby control conformational rearrangements. Orthosteric ligands are traditionally classified as either agonists, which promote activation, antagonists, which bind to the orthosteric site but do not alter basal activity, or inverse agonists that also reduce basal activity. However, this classical view of ligand efficacy is complicated by the fact that GPCRs can signal via several intracellular partners; for example, G-proteins or β-arrestins. Most agonists will activate several signaling pathways, but agonists with the ability to activate one specific intracellular partner have also been identified, a phenomenon referred to as biased signaling. Based on spectroscopy and structure determination studies, conformational changes in the receptor govern the activation of signaling pathways (Frei et al., 2020; Liu et al., 2012; Masureel et al., 2018), although the underlying molecular mechanisms remain elusive. Characterization of the allosteric process guiding interactions with intracellular partners is a major challenge and can only be fully understood by using a combination of different methodologies.

Figure 1

Download asset Open asset

The term microswitch, or molecular switch, describes local structural changes in the receptor that contribute to controlling activation and can, for example, involve side chain rotamers, movement of two domains relative to each other, or a helix twist. Two microswitches implicated in the activation of class A GPCRs are an outward displacement of the transmembrane helix 6 (TM6) and twist of the highly conserved N(7.49)P(7.50)xxY(7.53) motif (superscripts notation according to Ballesteros–Weinstein numbering; Ballesteros and Weinstein, 1995), which take part in the formation of the intracellular binding site. In the orthosteric site, microswitches are typically less conserved and depend on the type of ligand recognized by the receptor (Manglik and Kruse, 2017).

Considering the high dimensionality of a protein with over 300 interacting residues, it is difficult to identify relevant microswitches from the sequence or static experimental structures. Historically, sequence analysis and mutagenesis experiments (Gregorio et al., 2017; Lamichhane et al., 2020, Lamichhane et al., 2015; Manglik et al., 2015; Picard et al., 2019) have been used to characterize motifs important for signaling, but this approach may overlook the role of less conserved residues, water, and ions in ligand recognition and receptor activation (Chen et al., 2020). Molecular dynamics (MD) simulations can generate trajectories from experimental starting structures, capture the dynamics of all microswitches and allow us to derive the free energy landscapes governing the equilibrium between protein states. MD simulations have indeed been used extensively over the last decade to study GPCR activation (Bhattacharya and Vaidehi, 2010; Dror et al., 2011; Hu et al., 2019; Kohlhoff et al., 2014; Li et al., 2013; Miao and McCammon, 2016; Niesen et al., 2011; Shan et al., 2012; Tikhonova et al., 2013). Due to the computational cost of brute-force MD simulations, it is nearly impossible to obtain converged results without enhanced sampling methods, although the use of special purpose hardware (Dror et al., 2011) has pushed the boundaries of what is achievable by conventional MD simulations. Enhanced sampling strategies have emerged as an alternative, where exploration of the conformational landscape is promoted by the introduction of a bias in the simulations (Harpole and Delemotte, 2018). In a post-processing step, the bias can be adjusted for, and it is thus possible to derive theoretically exact results at a fraction of the cost of unbiased simulations. However, the ever-increasing size of simulation data and diverse conformations sampled by MD simulations makes it difficult to identify and determine the importance of microswitches by simple visualization of the conformational ensemble. Data-driven and machine learning approaches can help to condense the data and reduce human bias in the interpretation of the results (Fleetwood et al., 2020c; Hu et al., 2019).

In this study, we focus on the prototypical β₂ adrenergic receptor (β₂AR), which interacts with Gs proteins to trigger a cyclic adenosine monophosphate (cAMP) response, and arrestins, which control endocytosis and kinase activation (Jean-Charles et al., 2017). Both pathways are physiologically relevant and are modulated by therapeutic drugs. The β₂AR is a drug target for bronchoconstriction medication and was the first receptor crystallized in complex with a G protein (Rasmussen et al., 2011a; Rasmussen et al., 2011b). Experimental studies, including crystallography (Masureel et al., 2018; Rasmussen et al., 2011b; Ring et al., 2013), spectroscopy methods (Gregorio et al., 2017; Imai et al., 2020; Kofuku et al., 2012; Lamichhane et al., 2020; Liu et al., 2012), and computational methods (Provasi et al., 2011; Chen et al., 2020; Dror et al., 2011; Kohlhoff et al., 2014; Tikhonova et al., 2013), have investigated the activation mechanism of the β₂AR. Agonists bound in experimental structures show a key interaction with S207(5.46) (Chan et al., 2016) and an inward bulge of TM5 in the active state. In the TM domain between the orthosteric site and G protein-binding site, the connector region (Weis and Kobilka, 2018), partially overlapping with the P(5.50)I(3.40)F(6.44) motif, undergoes a rotameric change and thereby influences the hydrated cavity surrounding the conserved D79(2.50) (Imai et al., 2020), which in turn interacts with the conserved NPxxY motif in TM7 and reorients Y326(7.53) to form water-mediated interaction with Y219(5.58) (the Y–Y motif) (Latorraca et al., 2017). The combination of several microswitches leads to conformational changes that promote an outward movement of TM6 and binding of an intracellular binding partner, such as a G protein or arrestin. Understanding how ligands modulate individual microswitches could aid the development of biased agonists.

Enhanced sampling techniques have been used to characterize the activation mechanism of β₂AR, from early coarse-grained protocols (Bhattacharya and Vaidehi, 2010; Niesen et al., 2011) to more refined methodologies involving Gaussian accelerated MD (Tikhonova et al., 2013), metadynamics using path collective variables (CVs) derived from adiabatic biased MD simulations (Provasi et al., 2011), or adaptive sampling on cloud-based computing resources (Kohlhoff et al., 2014). Following in these footsteps, we recently introduced a version of the string with swarms of trajectories method designed to capture the activation pathway and the free energy landscapes along various microswitches (Fleetwood et al., 2020b).

In this study, thanks to our cost-effective computational approach, we have derived the activation free energy and characterized the details of the active-like state of the β₂AR (Figure 1a) in its ligand-free state and bound to six ligands with different efficacy profiles (Figure 1b), all of which were resolved bound to the β₂AR (Figure 1b) and several of which are clinically approved drugs (Woo and Robinson, 2015). The free energy landscapes revealed a stabilization of active-like states for the receptor bound to agonists and a shift toward inactive-like states for the receptor bound to antagonists or inverse agonists. Remarkably, we obtained a strong quantitative correlation between experimentally measured intracellular cAMP responses and the expectation values of the upper and transmembrane microswitches, highlighting the predictive power of our approach. In a second step, we introduce an adaptive sampling protocol developed to quantitatively sample the most stabilized states kinetically accessible from the activated starting structure (which we will refer to as the active-like state). Using dimensionality reduction techniques, we find that all ligands stabilize distinct receptor states and that ligands with similar pharmacological properties cluster together. Several of the microswitches considered to be of significance for GPCR activation, such as the NPxxY motif and the extracellular end of TM5, were automatically identified as important with our protocol. Combined with the activation free energies, our results show how ligands control the population of states. They modulate the conformational equilibrium by tuning important allosteric microswitches, in particular near the G protein-binding site. By inspecting the inter-residue contacts formed for different ligands, we identified an allosteric pathway between the two binding sites and a large heterogeneity of TM7 states. Our results thus build on the earlier use of enhanced sampling methods and demonstrate how such protocols combined with today’s computational capacities and availability of high-resolution structures in various states can provide insights into the structural basis of allosteric communication and ligand efficacy profiles, and potentially find use in the design of novel GPCR drug candidates.

Following the progress in GPCR research, it has become evident that a simple two-state model of activation is an oversimplification with considerable limitations. To explain biased and partial agonism, there is a need for a more comprehensive model. Many ligands have been characterized as full, partial, or biased agonists (van der Westhuizen et al., 2014). However, a systematic characterization of the molecular mechanisms which transmit this allosteric communication across the cell membrane remains elusive. Researchers have successfully managed to discriminate between arrestin and G protein bias using spectroscopic probes (Lamichhane et al., 2020; Liu et al., 2012), but the conformational changes and dynamics of many microswitches cannot be captured in a single measurement. MD simulations have the potential to provide additional insights thanks to the atomistic-level description they enable (Lamim Ribeiro and Filizola, 2019). Whereas the free energy landscapes clearly show that ligands influence several microswitches, a direct comparison between free energy profiles may be misleading if the other orthogonal microswitches are ignored. For example, the local minima in the connector ΔRMSD landscapes are located at similar positions for all ligands, but our analysis clearly shows that the overall conformational ensembles differ (Figure 2a and Figure 2—figure supplement 7b). Such projections onto single variables can thus obscure major differences in other microswitches. To address this limitation, we have used data-driven analysis methods, which are better suited for handling high-dimensional data than mere visual inspection, and found that there were indeed significant conformational differences in the states stabilized by the different ligands. Remarkably, this protocol automatically identified both receptor-specific and conserved motifs considered to be of significance for GPCR activation, such as the TM5 bulge and the NPxxY motif, as important. The different dimensionality reduction techniques found similar partitioning of the data, which strongly indicates that the results are not due to fortuitous parameter tuning or method choice. One of the remaining enigmas in GPCR research is to understand how the same overall activation mechanism can be conserved in spite of the fact that very different ligands are recognized by the family. Our machine learning-inspired data analysis protocol provides an unbiased approach to identify key features of activation for different receptor types.

The derivation of free energy landscapes and the corresponding microswitch expectation values provide a tool for estimating the stability of activation states, and thus also the relative efficacy of different ligands. Given the high correlations between microswitch expectation values and experimental data, we anticipate that microswitches located in the orthosteric and connector regions can be used for future predictions of ligand efficacy. Additionally, an advantage of the methods used in this study is that results were derived from several simulation replicas, which reduces the statistical error related to the stochastic nature of MD simulations on short time scales. Since we allow the strings to diffuse around the converged equilibrium pathway for many iterations, the statistical error in the microswitch expectation values and energy landscapes is small, although it may be somewhat underestimated since swarm trajectories launched from the same point are correlated to each other. The systematic error is likely bigger for reasons related to the choice of force field (Guvench and MacKerell, 2008). An important limitation is that the string method can only identify one out of multiple activation pathways. Our control simulation starting from the inactive crystal structure, which converged to a different active state, indeed demonstrated this issue and highlights the importance of starting from a well-chosen input pathway.

Understanding ligand control of receptor activation is of great interest in drug development. The conformational selection of intracellular binding partners enables the construction of molecules or nanobodies that have a high binding affinity only in combination with a specific ligand (Sencanski et al., 2019), which has been utilized in structure determination (Masureel et al., 2018; Rasmussen et al., 2011a; Ring et al., 2013; Staus et al., 2016). Similarly, a well-designed nanobody or allosteric modulator could enhance or block the binding of a specific compound (Staus et al., 2016). However, in this study, the ligands were bound to the receptor throughout the simulations and no intracellular binding partner was considered. GPCRs only assume their fully active state in the presence of an intracellular binding partner (Gregorio et al., 2017; Manglik et al., 2015; Nygaard et al., 2013). The β₂AR undergoes a basal activity, where it fluctuates between active-like and inactive states (Gregorio et al., 2017; Lamichhane et al., 2015), which can be inferred from the relatively low free energy difference between the basins in the energy landscapes. While extracellular and TM microswitches only require an agonist to maintain their active state conformations, intracellular microswitches interact with an intracellular binding partner directly and are more affected by its absence, which is likely why correlations between microswitch expectation value and Emax are worse for the intracellular microswitches.

Other important aspects of receptor–ligand interactions, such as identifying the binding pose of novel compounds or estimating a ligand’s binding free energy, cannot be estimated with this approach. There are other complementary techniques such as free energy perturbation methods (Cournia et al., 2017; Matricon et al., 2021; Matricon et al., 2017), which can be used to rigorously estimate binding affinity. Combined with enhanced sampling MD, we move toward having a complete toolkit for in silico drug design for development of high-affinity GPCR ligands with a specific efficacy.

Despite the increasing number of β₂AR structures, the receptor has not been solved in a form bound to arrestin. However, conformations observed in our simulations share properties with other arrestin-bound states observed for other receptors. For example, the supervised learning methods identified the C-terminal and H8 below the NPxxY motif as relatively important (Figure 4a, Figure 4—figure supplements 1–4). This region is known to be stabilized in ligand-dependent states for the angiotensin II type 1 receptor (Suomivuori et al., 2020; Wingler et al., 2019). While the reorganization of H8 may be a secondary effect due to modulation of the NPxxY motif, this region could be important for arrestin recruitment (Lally et al., 2017; Staus et al., 2018). Salmeterol’s distinct NPxxY state only formed in combination with a lost interaction between salmeterol and S207(5.46) and S203(5.43), which is remarkable considering that the two binding sites are believed to be only loosely coupled (Dror et al., 2011; Fleetwood et al., 2020b). A similar phenomenon has been reported in a recent structure of β₁AR (Lee et al., 2020), where the corresponding serines in the orthosteric site experienced weakened interactions to the biased agonist formoterol. The fact that Y(7.53) forms contacts with arrestin for β₁AR (Lee et al., 2020) suggests that our derived β₂AR conformation may have biased signaling properties.

Our results show that the activation pathways as well as the stabilized states are significantly altered upon ligand binding, and that ligands with shared efficacy profiles generate similar, albeit not identical, ensembles of states. It therefore cannot be taken for granted that two ligands which lead to a similar downstream response necessarily stabilize the same receptor conformations. As we considered several compounds in this study, similarities and differences between different compound classes emerged. The results in this study provide a good starting point for further analysis and allowed us to catch a glimpse of the complexity underlying receptor signaling. A thorough quantification of biased and partial agonism will require studies of ligands that stimulate various signaling pathways to different extents.

The data necessary to reproduce the findings presented in this paper can be found on OSF (DOI 10.17605/OSF.IO/6XPYV). The code used to run and analyze simulations has been deposited on GitHub (https://github.com/delemottelab/demystifying, https://github.com/delemottelab/gpcr-string-method-2019 and https://github.com/delemottelab/state-sampling; copies archived at https://archive.softwareheritage.org/swh:1:rev:e8527b52d5fbe0570cd391921ecda5aefceb797a/, https://archive.softwareheritage.org/swh:1:rev:bc3b7ce2e74e5ac95644d57a1b24f717a7ec74a4/, https://archive.softwareheritage.org/swh:1:rev:f0f56430ce581f0338771c126da212ecc2f218a0/).

1. 1. Abraham MJ
   2. Murtola T
   3. Schulz R
   4. Páll S
   5. Smith JC
   6. Hess B
   7. Lindahl E

   (2015) GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers

   SoftwareX 1-2:19–25.

   https://doi.org/10.1016/j.softx.2015.06.001

   • Google Scholar
2. 1. Ballesteros JA
   2. Weinstein H

   (1995) [19] Integrated methods for the construction of three-dimensional models and computational probing of structure-function relations in G protein-coupled receptors, in: receptor molecular biology, methods in neurosciences

   Methods in Neurosciences 25:366–428.

   https://doi.org/10.1016/S1043-9471(05)80049-7

   • Google Scholar
3. 1. Bhattacharya S
   2. Vaidehi N

   (2010) Computational mapping of the conformational transitions in agonist selective pathways of a G-protein coupled receptor

   Journal of the American Chemical Society 132:5205–5214.

   https://doi.org/10.1021/ja910700y

   • PubMed
   • Google Scholar
4. Book

   1. Borg H
   2. Groenen P

   (2005)

   Modern Multidimensional Scaling

   Springer .

   • Google Scholar
5. Book

   1. Breiman L
   2. Friedman JH
   3. Olshen RA
   4. Stone CJ

   (1984)

   Classification and Regression Trees

   Wadsworth & Brooks/Cole Advanced Books & Software.

   • Google Scholar
6. 1. Chan HC
   2. Filipek S
   3. Yuan S

   (2016) The principles of ligand specificity on beta-2-adrenergic receptor

   Scientific Reports 6:34736.

   https://doi.org/10.1038/srep34736

   • PubMed
   • Google Scholar
7. 1. Chen KM
   2. Keri D
   3. Barth P

   (2020) Computational design of G Protein-Coupled receptor allosteric signal transductions

   Nature Chemical Biology 16:77–86.

   https://doi.org/10.1038/s41589-019-0407-2

   • PubMed
   • Google Scholar
8. 1. Cherezov V
   2. Rosenbaum DM
   3. Hanson MA
   4. Rasmussen SG
   5. Thian FS
   6. Kobilka TS
   7. Choi HJ
   8. Kuhn P
   9. Weis WI
   10. Kobilka BK
   11. Stevens RC

   (2007) High-resolution crystal structure of an engineered human beta2-adrenergic G protein-coupled receptor

   Science 318:1258–1265.

   https://doi.org/10.1126/science.1150577

   • PubMed
   • Google Scholar
9. 1. Chung FZ
   2. Wang CD
   3. Potter PC
   4. Venter JC
   5. Fraser CM

   (1988) Site-directed mutagenesis and continuous expression of human beta-adrenergic receptors. identification of a conserved aspartate residue involved in agonist binding and receptor activation

   Journal of Biological Chemistry 263:4052–4055.

   https://doi.org/10.1016/S0021-9258(18)68888-X

   • PubMed
   • Google Scholar
10. 1. Cournia Z
    2. Allen B
    3. Sherman W

    (2017) Relative binding free energy calculations in drug discovery: recent advances and practical considerations

    Journal of Chemical Information and Modeling 57:2911–2937.

    https://doi.org/10.1021/acs.jcim.7b00564

    • PubMed
    • Google Scholar
11. 1. Dallongeville J
    2. Helbecque N
    3. Cottel D
    4. Amouyel P
    5. Meirhaeghe A

    (2003) The Gly16-->Arg16 and Gln27-->Glu27 polymorphisms of beta2-adrenergic receptor are associated with metabolic syndrome in men

    The Journal of Clinical Endocrinology and Metabolism 88:4862–4866.

    https://doi.org/10.1210/jc.2003-030173

    • PubMed
    • Google Scholar
12. 1. Dror RO
    2. Arlow DH
    3. Borhani DW
    4. Jensen MØ
    5. Piana S
    6. Shaw DE

    (2009) Identification of two distinct inactive conformations of the beta2-adrenergic receptor reconciles structural and biochemical observations

    PNAS 106:4689–4694.

    https://doi.org/10.1073/pnas.0811065106

    • PubMed
    • Google Scholar
13. 1. Dror RO
    2. Arlow DH
    3. Maragakis P
    4. Mildorf TJ
    5. Pan AC
    6. Xu H
    7. Borhani DW
    8. Shaw DE

    (2011) Activation mechanism of the 2-adrenergic receptor

    PNAS 108:18684–18689.

    https://doi.org/10.1073/pnas.1110499108

    • Google Scholar
14. Software

    1. Fleetwood O

    (2019a) Delemottelab/gpcr-String-Method-2019 Initial-Submission

    Delemottelab/gpcr-String-Method-2019 Initial-Submission.

    https://zenodo.org/record/2662201/export/json#.YBQHO1kzZ80
15. Software

    1. Fleetwood O

    (2019b) Delemottelab/demystifying: Updated Benchmarks, Zenodo

    Delemottelab/demystifying: Updated Benchmarks, Zenodo.

    https://zenodo.org/record/3269704#.YBQHV1kzZ80
16. Software

    1. Fleetwood O

    (2020a) Delemottelab/state-Sampling Initial-Release, Zenodo

    Delemottelab/state-Sampling Initial-Release, Zenodo.

    https://zenodo.org/record/3904591#.YBQHtlkzZ80
17. 1. Fleetwood O
    2. Matricon P
    3. Carlsson J
    4. Delemotte L

    (2020b) Energy landscapes reveal agonist control of G Protein-Coupled receptor activation via microswitches

    Biochemistry 59:880–891.

    https://doi.org/10.1021/acs.biochem.9b00842

    • PubMed
    • Google Scholar
18. 1. Fleetwood O
    2. Kasimova MA
    3. Westerlund AM
    4. Delemotte L

    (2020c) Molecular insights from conformational ensembles via machine learning

    Biophysical Journal 118:765–780.

    https://doi.org/10.1016/j.bpj.2019.12.016

    • PubMed
    • Google Scholar
19. 1. Flood E
    2. Boiteux C
    3. Lev B
    4. Vorobyov I
    5. Allen TW

    (2019) Atomistic simulations of membrane ion channel conduction, gating, and modulation

    Chemical Reviews 119:7737–7832.

    https://doi.org/10.1021/acs.chemrev.8b00630

    • PubMed
    • Google Scholar
20. 1. Frei JN
    2. Broadhurst RW
    3. Bostock MJ
    4. Solt A
    5. Jones AJY
    6. Gabriel F
    7. Tandale A
    8. Shrestha B
    9. Nietlispach D

    (2020) Conformational plasticity of ligand-bound and ternary GPCR complexes studied by 19F NMR of the β₁-adrenergic receptor

    Nature Communications 11:669.

    https://doi.org/10.1038/s41467-020-14526-3

    • PubMed
    • Google Scholar
21. 1. Gregorio GG
    2. Masureel M
    3. Hilger D
    4. Terry DS
    5. Juette M
    6. Zhao H
    7. Zhou Z
    8. Perez-Aguilar JM
    9. Hauge M
    10. Mathiasen S
    11. Javitch JA
    12. Weinstein H
    13. Kobilka BK
    14. Blanchard SC

    (2017) Single-molecule analysis of ligand efficacy in β2AR–G-protein activation

    Nature 547:68–73.

    https://doi.org/10.1038/nature22354

    • Google Scholar
22. 1. Guvench O
    2. MacKerell AD

    (2008) Comparison of protein force fields for molecular dynamics simulations

    Methods in Molecular Biology 443:63–88.

    https://doi.org/10.1007/978-1-59745-177-2_4

    • PubMed
    • Google Scholar
23. 1. Hanson MA
    2. Cherezov V
    3. Griffith MT
    4. Roth CB
    5. Jaakola VP
    6. Chien EY
    7. Velasquez J
    8. Kuhn P
    9. Stevens RC

    (2008) A specific cholesterol binding site is established by the 2.8 A structure of the human beta2-adrenergic receptor

    Structure 16:897–905.

    https://doi.org/10.1016/j.str.2008.05.001

    • PubMed
    • Google Scholar
24. 1. Harpole TJ
    2. Delemotte L

    (2018) Conformational landscapes of membrane proteins delineated by enhanced sampling molecular dynamics simulations

    Biochimica Et Biophysica Acta (BBA) - Biomembranes 1860:909–926.

    https://doi.org/10.1016/j.bbamem.2017.10.033

    • Google Scholar
25. 1. Harrigan MP
    2. Sultan MM
    3. Hernández CX
    4. Husic BE
    5. Eastman P
    6. Schwantes CR
    7. Beauchamp KA
    8. McGibbon RT
    9. Pande VS

    (2017) MSMBuilder: statistical models for biomolecular dynamics

    Biophysical Journal 112:10–15.

    https://doi.org/10.1016/j.bpj.2016.10.042

    • PubMed
    • Google Scholar
26. 1. Hauser AS
    2. Attwood MM
    3. Rask-Andersen M
    4. Schiöth HB
    5. Gloriam DE

    (2017) Trends in GPCR drug discovery: new agents, targets and indications

    Nature Reviews Drug Discovery 16:829–842.

    https://doi.org/10.1038/nrd.2017.178

    • PubMed
    • Google Scholar
27. Conference

    1. Ho TK

    (1995)

    Random decision forests

    Proceedings of 3rd International Conference on Document Analysis and Recognition. pp. 278–282.

    • Google Scholar
28. 1. Hu X
    2. Wang Y
    3. Hunkele A
    4. Provasi D
    5. Pasternak GW
    6. Filizola M

    (2019) Kinetic and thermodynamic insights into sodium ion translocation through the μ-opioid receptor from molecular dynamics and machine learning analysis

    PLOS Computational Biology 15:e1006689.

    https://doi.org/10.1371/journal.pcbi.1006689

    • PubMed
    • Google Scholar
29. 1. Huang J
    2. Rauscher S
    3. Nawrocki G
    4. Ran T
    5. Feig M
    6. de Groot BL
    7. Grubmüller H
    8. MacKerell AD

    (2017) CHARMM36m: an improved force field for folded and intrinsically disordered proteins

    Nature Methods 14:71–73.

    https://doi.org/10.1038/nmeth.4067

    • PubMed
    • Google Scholar
30. 1. Imai S
    2. Yokomizo T
    3. Kofuku Y
    4. Shiraishi Y
    5. Ueda T
    6. Shimada I

    (2020) Structural equilibrium underlying ligand-dependent activation of β₂-adrenoreceptor

    Nature Chemical Biology 16:430–439.

    https://doi.org/10.1038/s41589-019-0457-5

    • PubMed
    • Google Scholar
31. 1. Isberg V
    2. Vroling B
    3. van der Kant R
    4. Li K
    5. Vriend G
    6. Gloriam D

    (2014) GPCRDB: an information system for G protein-coupled receptors

    Nucleic Acids Research 42:D422–D425.

    https://doi.org/10.1093/nar/gkt1255

    • PubMed
    • Google Scholar
32. 1. Jean-Charles PY
    2. Kaur S
    3. Shenoy SK

    (2017) G Protein-Coupled receptor signaling through β-Arrestin-Dependent mechanisms

    Journal of Cardiovascular Pharmacology 70:142–158.

    https://doi.org/10.1097/FJC.0000000000000482

    • PubMed
    • Google Scholar
33. 1. Jorgensen WL
    2. Chandrasekhar J
    3. Madura JD
    4. Impey RW
    5. Klein ML

    (1983) Comparison of simple potential functions for simulating liquid water

    The Journal of Chemical Physics 79:926–935.

    https://doi.org/10.1063/1.445869

    • Google Scholar
34. 1. Klauda JB
    2. Venable RM
    3. Freites JA
    4. O'Connor JW
    5. Tobias DJ
    6. Mondragon-Ramirez C
    7. Vorobyov I
    8. MacKerell AD
    9. Pastor RW

    (2010) Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types

    The Journal of Physical Chemistry B 114:7830–7843.

    https://doi.org/10.1021/jp101759q

    • PubMed
    • Google Scholar
35. 1. Kofuku Y
    2. Ueda T
    3. Okude J
    4. Shiraishi Y
    5. Kondo K
    6. Maeda M
    7. Tsujishita H
    8. Shimada I

    (2012) Efficacy of the β₂-adrenergic receptor is determined by conformational equilibrium in the transmembrane region

    Nature Communications 3:1045.

    https://doi.org/10.1038/ncomms2046

    • PubMed
    • Google Scholar
36. 1. Kohlhoff KJ
    2. Shukla D
    3. Lawrenz M
    4. Bowman GR
    5. Konerding DE
    6. Belov D
    7. Altman RB
    8. Pande VS

    (2014) Cloud-based simulations on google exacycle reveal ligand modulation of GPCR activation pathways

    Nature Chemistry 6:15–21.

    https://doi.org/10.1038/nchem.1821

    • PubMed
    • Google Scholar
37. 1. Kullback S
    2. Leibler RA

    (1951) On information and sufficiency

    The Annals of Mathematical Statistics 22:79–86.

    https://doi.org/10.1214/aoms/1177729694

    • Google Scholar
38. 1. Lally CC
    2. Bauer B
    3. Selent J
    4. Sommer ME

    (2017) C-edge loops of arrestin function as a membrane anchor

    Nature Communications 8:14258.

    https://doi.org/10.1038/ncomms14258

    • PubMed
    • Google Scholar
39. 1. Lamichhane R
    2. Liu JJ
    3. Pljevaljcic G
    4. White KL
    5. van der Schans E
    6. Katritch V
    7. Stevens RC
    8. Wüthrich K
    9. Millar DP

    (2015) Single-molecule view of basal activity and activation mechanisms of the G protein-coupled receptor β2ar

    PNAS 112:14254–14259.

    https://doi.org/10.1073/pnas.1519626112

    • PubMed
    • Google Scholar
40. 1. Lamichhane R
    2. Liu JJ
    3. White KL
    4. Katritch V
    5. Stevens RC
    6. Wüthrich K
    7. Millar DP

    (2020) Biased signaling of the G-Protein-Coupled receptor β₂AR Is Governed by Conformational Exchange Kinetics

    Structure 28:371–377.

    https://doi.org/10.1016/j.str.2020.01.001

    • PubMed
    • Google Scholar
41. 1. Lamim Ribeiro JM
    2. Filizola M

    (2019) Allostery in G protein-coupled receptors investigated by molecular dynamics simulations

    Current Opinion in Structural Biology 55:121–128.

    https://doi.org/10.1016/j.sbi.2019.03.016

    • PubMed
    • Google Scholar
42. 1. Latorraca NR
    2. Venkatakrishnan AJ
    3. Dror RO

    (2017) GPCR dynamics: structures in motion

    Chemical Reviews 117:139–155.

    https://doi.org/10.1021/acs.chemrev.6b00177

    • PubMed
    • Google Scholar
43. 1. Lee J
    2. Cheng X
    3. Swails JM
    4. Yeom MS
    5. Eastman PK
    6. Lemkul JA
    7. Wei S
    8. Buckner J
    9. Jeong JC
    10. Qi Y
    11. Jo S
    12. Pande VS
    13. Case DA
    14. Brooks CL
    15. MacKerell AD
    16. Klauda JB
    17. Im W

    (2016) CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field

    Journal of Chemical Theory and Computation 12:405–413.

    https://doi.org/10.1021/acs.jctc.5b00935

    • PubMed
    • Google Scholar
44. Preprint

    1. Lee Y
    2. Warne T
    3. Nehmé R
    4. Pandey S
    5. Dwivedi-Agnihotri H
    6. Edwards PC
    7. García-Nafría J
    8. Leslie AGW
    9. Shukla AK
    10. Tate CG

    (2020) Molecular determinants of β-arrestin coupling to formoterol-bound β₁ -adrenoceptor

    bioRxiv.

    https://doi.org/10.1101/2020.03.27.011585

    • Google Scholar
45. 1. Lev B
    2. Murail S
    3. Poitevin F
    4. Cromer BA
    5. Baaden M
    6. Delarue M
    7. Allen TW

    (2017) String method solution of the gating pathways for a pentameric ligand-gated ion channel

    PNAS 114:E4158–E4167.

    https://doi.org/10.1073/pnas.1617567114

    • PubMed
    • Google Scholar
46. 1. Li J
    2. Jonsson AL
    3. Beuming T
    4. Shelley JC
    5. Voth GA

    (2013) Ligand-Dependent Activation and Deactivation of the Human Adenosine A _2AReceptor

    Journal of the American Chemical Society 135:8749–8759.

    https://doi.org/10.1021/ja404391q

    • Google Scholar
47. 1. Liu JJ
    2. Horst R
    3. Katritch V
    4. Stevens RC
    5. Wuthrich K

    (2012) Biased Signaling Pathways in 2-Adrenergic Receptor Characterized by 19F-NMR

    Science 335:1106–1110.

    https://doi.org/10.1126/science.1215802

    • Google Scholar
48. 1. Liu X
    2. Kaindl J
    3. Korczynska M
    4. Stößel A
    5. Dengler D
    6. Stanek M
    7. Hübner H
    8. Clark MJ
    9. Mahoney J
    10. Matt RA
    11. Xu X
    12. Hirata K
    13. Shoichet BK
    14. Sunahara RK
    15. Kobilka BK
    16. Gmeiner P

    (2020) An allosteric modulator binds to a conformational hub in the β2 adrenergic receptor

    Nature Chemical Biology 16:749–755.

    https://doi.org/10.1038/s41589-020-0549-2

    • PubMed
    • Google Scholar
49. 1. Ma X
    2. Hu Y
    3. Batebi H
    4. Heng J
    5. Xu J
    6. Liu X
    7. Niu X
    8. Li H
    9. Hildebrand PW
    10. Jin C
    11. Kobilka BK

    (2020) Analysis of β2AR-Gs and β2AR-Gi complex formation by NMR spectroscopy

    PNAS 117:23096–23105.

    https://doi.org/10.1073/pnas.2009786117

    • Google Scholar
50. 1. Manglik A
    2. Kim TH
    3. Masureel M
    4. Altenbach C
    5. Yang Z
    6. Hilger D
    7. Lerch MT
    8. Kobilka TS
    9. Thian FS
    10. Hubbell WL
    11. Prosser RS
    12. Kobilka BK

    (2015) Structural insights into the dynamic process of β2-Adrenergic receptor signaling

    Cell 161:1101–1111.

    https://doi.org/10.1016/j.cell.2015.04.043

    • PubMed
    • Google Scholar
51. 1. Manglik A
    2. Kruse AC

    (2017) Structural basis for G Protein-Coupled receptor activation

    Biochemistry 56:5628–5634.

    https://doi.org/10.1021/acs.biochem.7b00747

    • PubMed
    • Google Scholar
52. 1. Masureel M
    2. Zou Y
    3. Picard LP
    4. van der Westhuizen E
    5. Mahoney JP
    6. Rodrigues J
    7. Mildorf TJ
    8. Dror RO
    9. Shaw DE
    10. Bouvier M
    11. Pardon E
    12. Steyaert J
    13. Sunahara RK
    14. Weis WI
    15. Zhang C
    16. Kobilka BK

    (2018) Structural insights into binding specificity, efficacy and Bias of a β₂AR partial agonist

    Nature Chemical Biology 14:1059–1066.

    https://doi.org/10.1038/s41589-018-0145-x

    • PubMed
    • Google Scholar
53. 1. Matricon P
    2. Ranganathan A
    3. Warnick E
    4. Gao ZG
    5. Rudling A
    6. Lambertucci C
    7. Marucci G
    8. Ezzati A
    9. Jaiteh M
    10. Dal Ben D
    11. Jacobson KA
    12. Carlsson J

    (2017) Fragment optimization for GPCRs by molecular dynamics free energy calculations: probing druggable subpockets of the A 2A adenosine receptor binding site

    Scientific Reports 7:6398.

    https://doi.org/10.1038/s41598-017-04905-0

    • PubMed
    • Google Scholar
54. 1. Matricon P
    2. Suresh RR
    3. Gao Z-G
    4. Panel N
    5. Jacobson KA
    6. Carlsson J

    (2021) Ligand design by targeting a binding site water

    Chemical Science 12:960–968.

    https://doi.org/10.1039/D0SC04938G

    • Google Scholar
55. 1. Miao Y
    2. McCammon JA

    (2016) Graded activation and free energy landscapes of a muscarinic G-protein-coupled receptor

    PNAS 113:12162–12167.

    https://doi.org/10.1073/pnas.1614538113

    • PubMed
    • Google Scholar
56. 1. Millman KJ
    2. Aivazis M

    (2011) Python for scientists and engineers

    Computing in Science & Engineering 13:9–12.

    https://doi.org/10.1109/MCSE.2011.36

    • Google Scholar
57. 1. Montavon G
    2. Samek W
    3. Müller K-R

    (2018) Methods for interpreting and understanding deep neural networks

    Digital Signal Processing 73:1–15.

    https://doi.org/10.1016/j.dsp.2017.10.011

    • Google Scholar
58. 1. Niesen MJM
    2. Bhattacharya S
    3. Vaidehi N

    (2011) The Role of Conformational Ensembles in Ligand Recognition in G-Protein Coupled Receptors

    Journal of the American Chemical Society 133:13197–13204.

    https://doi.org/10.1021/ja205313h

    • Google Scholar
59. 1. Nygaard R
    2. Zou Y
    3. Dror RO
    4. Mildorf TJ
    5. Arlow DH
    6. Manglik A
    7. Pan AC
    8. Liu CW
    9. Fung JJ
    10. Bokoch MP
    11. Thian FS
    12. Kobilka TS
    13. Shaw DE
    14. Mueller L
    15. Prosser RS
    16. Kobilka BK

    (2013) The Dynamic Process of β2-Adrenergic Receptor Activation

    Cell 152:532–542.

    https://doi.org/10.1016/j.cell.2013.01.008

    • Google Scholar
60. 1. Pan AC
    2. Sezer D
    3. Roux B

    (2008) Finding Transition Pathways Using the String Method with Swarms of Trajectories

    The Journal of Physical Chemistry B 112:3432–3440.

    https://doi.org/10.1021/jp0777059

    • Google Scholar
61. Preprint

    1. Pedregosa F
    2. Varoquaux G
    3. Gramfort A

    (2011) Scikit-learn: machine learning in Python

    arXiv.

    https://arxiv.org/abs/1201.0490

    • Google Scholar
62. 1. Picard L-P
    2. Schonegge A-M
    3. Bouvier M

    (2019) Structural Insight into G Protein-Coupled Receptor Signaling Efficacy and Bias between Gs and β-Arrestin

    ACS Pharmacology & Translational Science 2:148–154.

    https://doi.org/10.1021/acsptsci.9b00012

    • Google Scholar
63. 1. Provasi D
    2. Artacho MC
    3. Negri A
    4. Mobarec JC
    5. Filizola M

    (2011) Ligand-Induced Modulation of the Free-Energy Landscape of G Protein-Coupled Receptors Explored by Adaptive Biasing Techniques

    PLOS Computational Biology 7:e1002193.

    https://doi.org/10.1371/journal.pcbi.1002193

    • Google Scholar
64. 1. Rasmussen SGF
    2. Choi H-J
    3. Fung JJ
    4. Pardon E
    5. Casarosa P
    6. Chae PS
    7. DeVree BT
    8. Rosenbaum DM
    9. Thian FS
    10. Kobilka TS
    11. Schnapp A
    12. Konetzki I
    13. Sunahara RK
    14. Gellman SH
    15. Pautsch A
    16. Steyaert J
    17. Weis WI
    18. Kobilka BK

    (2011a) Structure of a nanobody-stabilized active state of the β2 adrenoceptor

    Nature 469:175–180.

    https://doi.org/10.1038/nature09648

    • Google Scholar
65. 1. Rasmussen SGF
    2. DeVree BT
    3. Zou Y
    4. Kruse AC
    5. Chung KY
    6. Kobilka TS
    7. Thian FS
    8. Chae PS
    9. Pardon E
    10. Calinski D
    11. Mathiesen JM
    12. Shah STA
    13. Lyons JA
    14. Caffrey M
    15. Gellman SH
    16. Steyaert J
    17. Skiniotis G
    18. Weis WI
    19. Sunahara RK
    20. Kobilka BK

    (2011b) Crystal structure of the β2 adrenergic receptor–Gs protein complex

    Nature 477:549–555.

    https://doi.org/10.1038/nature10361

    • Google Scholar
66. 1. Ring AM
    2. Manglik A
    3. Kruse AC
    4. Enos MD
    5. Weis WI
    6. Garcia KC
    7. Kobilka BK

    (2013) Adrenaline-activated structure of β2-adrenoceptor stabilized by an engineered nanobody

    Nature 502:575–579.

    https://doi.org/10.1038/nature12572

    • Google Scholar
67. Software

    1. Schubert E

    (2017)

    interpretation - Clustering on the output of t-SNE

    interpretation - Clustering on the output of t-SNE.
68. 1. Sencanski M
    2. Glisic S
    3. Šnajder M
    4. Veljkovic N
    5. Poklar Ulrih N
    6. Mavri J
    7. Vrecl M

    (2019) Computational design and characterization of nanobody-derived peptides that stabilize the active conformation of the β2-adrenergic receptor (β2-AR)

    Scientific Reports 9:16555.

    https://doi.org/10.1038/s41598-019-52934-8

    • Google Scholar
69. 1. Shan J
    2. Khelashvili G
    3. Mondal S
    4. Mehler EL
    5. Weinstein H

    (2012) Ligand-Dependent Conformations and Dynamics of the Serotonin 5-HT2A Receptor Determine Its Activation and Membrane-Driven Oligomerization Properties

    PLOS Computational Biology 8:e1002473.

    https://doi.org/10.1371/journal.pcbi.1002473

    • Google Scholar
70. Conference

    1. Smolensky P

    (1986)

    Information processing in dynamical systems: foundations of harmony theory

    Parallel Distributed Processing: Explorations in the Microstructure of Cognition: Foundations. pp. 194–281.

    • Google Scholar
71. 1. Staus DP
    2. Strachan RT
    3. Manglik A
    4. Pani B
    5. Kahsai AW
    6. Kim TH
    7. Wingler LM
    8. Ahn S
    9. Chatterjee A
    10. Masoudi A
    11. Kruse AC
    12. Pardon E
    13. Steyaert J
    14. Weis WI
    15. Prosser RS
    16. Kobilka BK
    17. Costa T
    18. Lefkowitz RJ

    (2016) Allosteric nanobodies reveal the dynamic range and diverse mechanisms of G-protein-coupled receptor activation

    Nature 535:448–452.

    https://doi.org/10.1038/nature18636

    • Google Scholar
72. 1. Staus DP
    2. Wingler LM
    3. Choi M
    4. Pani B
    5. Manglik A
    6. Kruse AC
    7. Lefkowitz RJ

    (2018) Sortase ligation enables homogeneous GPCR phosphorylation to reveal diversity in β-arrestin coupling

    PNAS 115:3834–3839.

    https://doi.org/10.1073/pnas.1722336115

    • Google Scholar
73. 1. Suomivuori C-M
    2. Latorraca NR
    3. Wingler LM
    4. Eismann S
    5. King MC
    6. Kleinhenz ALW
    7. Skiba MA
    8. Staus DP
    9. Kruse AC
    10. Lefkowitz RJ
    11. Dror RO

    (2020) Molecular mechanism of biased signaling in a prototypical G protein–coupled receptor

    Science 367:881–887.

    https://doi.org/10.1126/science.aaz0326

    • Google Scholar
74. 1. Tikhonova IG
    2. Selvam B
    3. Ivetac A
    4. Wereszczynski J
    5. McCammon JA

    (2013) Simulations of Biased Agonists in the β ₂ Adrenergic Receptor with Accelerated Molecular Dynamics

    Biochemistry 52:5593–5603.

    https://doi.org/10.1021/bi400499n

    • Google Scholar
75. 1. Tipping ME
    2. Bishop CM

    (1999) Probabilistic principal component analysis

    Journal of the Royal Statistical Society: Series B 61:611–622.

    https://doi.org/10.1111/1467-9868.00196

    • Google Scholar
76. 1. Tribello GA
    2. Bonomi M
    3. Branduardi D
    4. Camilloni C
    5. Bussi G

    (2014) PLUMED 2: new feathers for an old bird

    Computer Physics Communications 185:604–613.

    https://doi.org/10.1016/j.cpc.2013.09.018

    • Google Scholar
77. 1. Ueda T
    2. Kofuku Y
    3. Okude J
    4. Imai S
    5. Shiraishi Y
    6. Shimada I

    (2019) Function-related conformational dynamics of G protein-coupled receptors revealed by NMR

    Biophysical Reviews 11:409–418.

    https://doi.org/10.1007/s12551-019-00539-w

    • PubMed
    • Google Scholar
78. 1. van der Maaten L
    2. Hinton G

    (2008)

    Visualizing data using t-SNE

    Journal of Machine Learning Research : JMLR 9:2579–2605.

    • Google Scholar
79. 1. van der Westhuizen ET
    2. Breton B
    3. Christopoulos A
    4. Bouvier M

    (2014) Quantification of ligand Bias for clinically relevant β2-adrenergic receptor ligands: implications for drug taxonomy

    Molecular Pharmacology 85:492–509.

    https://doi.org/10.1124/mol.113.088880

    • PubMed
    • Google Scholar
80. 1. Venkatakrishnan AJ
    2. Ma AK
    3. Fonseca R
    4. Latorraca NR
    5. Kelly B
    6. Betz RM
    7. Asawa C
    8. Kobilka BK
    9. Dror RO

    (2019) Diverse GPCRs exhibit conserved water networks for stabilization and activation

    PNAS 116:3288–3293.

    https://doi.org/10.1073/pnas.1809251116

    • PubMed
    • Google Scholar
81. 1. Wacker D
    2. Fenalti G
    3. Brown MA
    4. Katritch V
    5. Abagyan R
    6. Cherezov V
    7. Stevens RC

    (2010) Conserved binding mode of human beta2 adrenergic receptor inverse agonists and antagonist revealed by X-ray crystallography

    Journal of the American Chemical Society 132:11443–11445.

    https://doi.org/10.1021/ja105108q

    • PubMed
    • Google Scholar
82. 1. Weis WI
    2. Kobilka BK

    (2018) The molecular basis of G Protein-Coupled receptor activation

    Annual Review of Biochemistry 87:897–919.

    https://doi.org/10.1146/annurev-biochem-060614-033910

    • PubMed
    • Google Scholar
83. 1. Wingler LM
    2. Elgeti M
    3. Hilger D
    4. Latorraca NR
    5. Lerch MT
    6. Staus DP
    7. Dror RO
    8. Kobilka BK
    9. Hubbell WL
    10. Lefkowitz RJ

    (2019) Angiotensin analogs with divergent Bias stabilize distinct receptor conformations

    Cell 176:468–478.

    https://doi.org/10.1016/j.cell.2018.12.005

    • PubMed
    • Google Scholar
84. Book

    1. Woo TM
    2. Robinson MV

    (2015)

    Pharmacotherapeutics for Advanced Practice Nurse Prescribers

    F. A. Davis Company.

    • Google Scholar

Acceptance summary:

This manuscript reports advanced molecular dynamics simulations to describe

β₂ adrenergic receptor modulation, expanding significantly on existing knowledge. The study has made use of an atomistic string method to measure the effects of agonists, antagonists and inverse agonists and to understand how ligands affect GPCR activity. The authors have presented sufficient analysis to demonstrate statistical significance of the data and have made connections to experimental measurements that are well described in the revised manuscript. Overall, this is a high-level computational study of biological significance that will be impressive to many eLife readers.

Decision letter after peer review:

Thank you for submitting your article "Identification of ligand-specific G protein-coupled receptor states and prediction of efficacy via data-driven modeling" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by José Faraldo-Gómez as the Senior Editor. The reviewers have opted to remain anonymous.

Based on the reviews received, which are enclosed below, as well as a follow-up discussion amongst reviewers and editors, we regret to inform you that we cannot publish your manuscript in its current form. However, the reviewers have agreed to re-evaluate a revised version of the manuscript that convincingly addresses the concerns raised. While we encourage you to submit a revision, please note that a positive outcome is not guaranteed.

As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

The reviewers have appreciated many aspects of the string method solution for the allosteric modulation of the β₂ adrenergic receptor, but were concerned about the reliability of the simulations, their connections to experiments and the absence of discussion on the relationship to past GPCR activation studies. Concern is expressed about the use of Emax as the sole proxy for activated conformation, as well as a lack of connection of free energy calculations to experiments. A critical point to address is the absence of a description of the statistical reliability and error estimates from the string method solutions, as well as the dependence of the converged solutions on the starting conformation. This necessitates at least one demonstration that reliable results emerge with a very different initial structure. A significantly revised manuscript would need to adequately address the above concerns to the satisfaction of the reviewers.

Essential revisions:

The full reviews have been included below. Essential revisions include:

* Proof of independence of string solutions on starting structure.

* Proof of statistical reliability and error estimates.

* Improved connections to experiments that do not rely solely on Emax values.

* Improved relationship to past studies.

Reviewer #1:

This manuscript reports advanced MD simulation to describe the allosteric modulation of the β₂ adrenergic receptor, expanding on the authors' recent simulations (Biochemistry 2020, 59:880-891) that explored communication of ligand binding to the G protein-binding site via "microswitches" within the protein. This study uses the same string method to examine the effects of full and biased agonists, antagonists and inverse agonists using previously solved Xray structures. The study also calls on a series of data-based analysis methods to seek answers to the question of how different ligands communicate their changes and affect GPCR activity. Observation of stabilization of active states for agonists and inactive states for other ligands is an important achievement for MD simulation, even if reproducing known experimental results. Also, characterizing how the distribution of kinetically stabilized (active-like) states are controlled by ligands, and correlations between microswitch expectation values and cAMP response experiments are important. However, the manuscript is highly technical and jargon-based such that many readers will not follow the analysis. The statistical reliability and dependence of the converged string solutions on starting conformation do not appear to have been tested. Finally, relationship between this and the previous manuscript (Fleetwood et al., 2020), as well as how results distinguish themselves from knowledge in the field could be made more clear.

String methods are notorious for depending on starting structure as well as the number and choice of collective variables, and may converge slowly on the optimal path, if at all. The authors do not show plots of convergence (how each PMF changes over #iterations), and do not appear to undertake tests for reproducibility. i.e. Starting from very different initial structures to show that the same final pathways emerge, leading to consistent free energy profiles. For example, starting again with an initial structure for an inverse agonist instead of one for an agonist and show the same results are computed. I note that in the Discussion it is said that several simulation replicas were used to reduce error. But to what does this apply (apparently not to the swarms)? Also, collective variables (CVs) are those defined in the previous study (Fleetwood et al., 2020) but what if the set of CVs is changed or the number of CVs is reduced/increased? While I do not suggest redoing with a different set, how does the reader judge the robustness of the CV-reliant results?

The method for free energy calculation uses a transition matrix on a CV grid with stationary solution. In the subsection “Free energy estimation”, the authors cite Pan, Sezer and Roux, 2008 (and Fleetwood et al., 2020), but it is not clear that Pan, Sezer and Roux, 2008, did this, and in fact I think the first to do this from a swarms of trajectories solution was Lev et al., 2017 (see also Flood, et al., 2019), not cited here.

Having computed free energy projections along collective variables, receptor stability in relation to efficacy instead relies on estimation values for variables and not free energies, and the justifications for this could be improved. The authors write that the width and depth of the basins determine the most probable state projected into that space. This sentence is not clear, but its meaning is important, because ultimately the relative values of free energies are discussed. If the 2D maps are valid projections of the full configurational space (albeit with sampling guided/biased by the string), then Boltzmann integrals over each basin should yield a valid equilibrium constant. I presume the concern is that projections along different pairs of CVs can lead to different apparent free energies, because different CVs map out different proportions of phase space, and potentially envelop multiple states of the system, but does this mean a Boltzmann integral over a site is not a true thermodynamic quantity? In the Discussion, the authors return to add that examination of one projection can overlook what is happening in other coordinates/switches, but more discussion/justification is needed. Related, the authors state that they have "accurately captured the relative stability of states", but this is confusing given the relative stability of states from the maps appears to have been discredited.

Instead of using free energies, the authors correlate "expectation values" of their CVs to experiments in the bottom of Figure 2. I assume a plot against free energy was abandoned as it was not working as planned? Ideally one would want to see free energies plotted against experimental efficacy, because estimation values of variables such as TM5…, may correlate to efficacy, but not uniquely map to shifting state equilibrium.

Many analysis tasks were completed by existing packages for which the meanings are not obvious. e.g. Demystifying, Scikit-learn… Materials and methods, even if accepted, deserve a sentence to explain and motivate. e.g. CVs were short inverse inter-residue distances used to train a restructured Boltzmann machine (the principles of this machine could be explained simply). The subsection “Data driven analysis reveals that ligands stabilize unique states” is highly technical and not well explained. While many will know about PCA, less will be known about MDS and T-SNE. The ability to identify signaling hotspots in Figure 4 is impressive, leading to a model for allosteric communication. But the machine learning approaches used come across as black box and require better physical interpretation.

Reviewer #2:

This manuscript describes the use of enhanced sampling molecular dynamics to calculate free energy landscapes for the β₂ adrenergic receptor. A particular focus is the conformational dynamics and thermodynamics of microswitches that change conformation upon receptor activation. A variety of ligands are investigated to identify shared and divergent features of receptor conformational modulation. Unbiased methods are shown to identify known conformational switches as key regulators of receptor activation.

Overall the manuscript is interesting and clearly written. Figures are very clear and well presented. The subject matter has been very extensively studied in the GPCR family as a whole and in the β₂ receptor in particular, and the results presented here largely align with existing understanding of GPCR activation and conformational regulation by ligands. A major concern is the question of whether the results presented here truly enhance our understanding of GPCR activation, or simply confirm things that are already known. Many groups have published detailed analysis of similar questions to those presented here, including the Dror, Nagarajan, and McCammon groups, among others. There is very little discussion of this prior work, which makes it difficult to see how the present manuscript fits within the broader context of GPCR activation molecular dynamics analysis.

A significant technical concern is the reliance on cAMP Emax values as the sole experimental validation. A prospective experimental test of hypotheses generated here would significantly enhance the manuscript. Barring this, a more extensive comparison of computational results with experimental data is essential in my view. The cAMP pathway is highly amplified, which may confound analysis linking Emax values directly to conformational equilibria. Manglik et al., 2015, presents an actual biophysical analysis of conformational equilibria (albeit with fewer ligands) and comparison to this would be helpful. Do relative energy predictions match these experimental observations?

A smaller issue is that the utility of the results presented here appears to be a bit overstated. For example, it is claimed (Introduction) that the results provide insight into how ligands with specific efficacy profiles can be designed. To support this statement, the authors should present examples of compounds they have designed based on their results, together with experimental data confirming these molecules have the intended efficacy profiles.

Reviewer #3:

This is an interesting paper, and I like to attempt to connect analysis of allostery to function. However, I'm extremely concerned about statistical uncertainty - it's not really discussed, and it would be easy to chalk all of the results up to limited sampling. It will be important for the authors to demonstrate this isn't the case.

Basically, my concern is that there's essentially no mention of statistical uncertainty or convergence anywhere in the document. One of the major claims is that the different agonists each populate a distinct substate - this is an incredibly important and interesting observation if true, but could also be explained by saying each simulation wandered in its own space and didn't have time to explore anywhere else. If it were run again, it might wander into a totally different place. I don't have a great feel for how rapidly swarms explore, but I do know the total sampling time of 1.4 µs per ligand sounds awfully small. There are examples in the literature showing that conventional simulations several times this length are not converged as far the configurations in the ligand binding pocket go (for example, Leioatts et al., Biophysical Journal, 2015, or some of the work from Dror's group).

I have a couple of ideas for how the authors could convince me this isn't just a sampling artifact. The best would probably be to pick one ligand and redo the whole calculation 4 more times, start to finish - looking at the variation between those replicates would be a decent estimate of the uncertainty. Ideally, they'd do this for all of the ligands, but I recognize that's not a reasonable request, which is why I suggest doing it for 1 ligand.

A weaker test would be to break the individual ligand calculations into blocks, somehow - I'm not totally sure how to do it - and show that the blocks are self-similar. If they're all wandering through the whole of the blob in Figure 3, you're more likely to be ok. If each block populates a discrete chunk of the blob (and the overall swarm data doesn't revisit), then there's a big problem.

I'd want to see error estimates on the free energies in Figure 2, plus a redo of the dimension reduction in Figure 3 to see if the ligands still separate more than the replicates of 1 ligand.

On a more minor note, the use of p-values in the subsection “Ligands control efficacy by reshaping microswitches’ probability distributions”, isn't really correct. The point the authors are trying to make is that the computed value doesn't predict the experiment (for good reason), and the correlation coefficients make that point for you.

Essential revisions:

The full reviews have been included below. Essential revisions include:

* Proof of independence of string solutions on starting structure.

* Proof of statistical reliability and error estimates.

* Improved connections to experiments that do not rely solely on Emax values.

* Improved relationship to past studies.

1. Proof of statistical reliability and error estimates.

We have performed a number of steps to better demonstrate the statistical reliability of our results:

String simulations

The free energy estimates constitute the most important result derived from our string simulations (Fig.2). To better assess their convergence, we analyzed the convergence of the underlying transition matrix, T, of the swarm trajectories’ transitions between microstates. Since many factors go into constructing T, the impact of all the statistical uncertainties on the final free energy landscape cannot be directly derived. Instead, we applied a Bayesian approach and used Metropolis Markov chain Monte Carlo (MCMC) to sample over the posterior distribution of transition matrices (Harrigan et al., 2017), and thereby obtained a distribution of free energy profiles. From these we could compute the mean, standard deviation, whiskers and outliers etc., which are shown in Figure 2–figure supplement 7A and B. Although this convergence assessment was already included in the original manuscript, we acknowledge that we did not properly inform the reader of this analysis. This has been corrected in this revised manuscript.

As shown in Figure 2–figure supplement 7A and B, the error bars are typically small near the basins, and thus the standard error of the microswitch expectation values will also be small. The error bars may therefore be difficult to discern for some values in the correlation figures in Fig. 2. This use of bootstrapping analysis methods to assess convergence is common practice in the analysis of Markov State Modelling (Prinz et al., 2011). Note that if we have many transitions across free energy barriers, this approach is prone to underestimate the total error since systematic errors due to force field are not taken into account.We also monitored the convergence of the pathway for every simulation in CV space. Author response image 1 shows that after approximately 100 iterations the string does not move on average. Instead, it fluctuates around an equilibrium value.

Author response image 1

Download asset Open asset

State sampling simulations

The state sampling simulations form the basis for Fig. 3-5. It is an iterative sampling approach where we track the center point in CV space of multiple walkers. By iteratively terminating walkers which drift far away from the center (since they may be considered as statistical outliers), we monitor how the center point converges to a fixed value. Figure 3 – figure supplement 3, which shows the distance between center points, indicates convergence for all systems within 8 iterations. This figure has now been updated to also show the standard error of the walkers’ distance to the center point. The distance between center points is smaller than the error and the walkers’ center is fluctuating around their equilibrium position. The distribution of walkers around this position is what we refer to as the pre-active state.

Following Reviewer #3’s suggestion, we also performed cross validation with a subset of the walkers and showed that the results are self-similar. In other words, by only considering one third of the walkers at the time, we arrived at very similar dimensionality reduction figures. We thus found that our claims regarding the similarities and differences between agonists and nonagonists still held after these convergence assessments. These results have been included in the revised Figure 3 – figure supplement 3.

Investigation of independence of string solutions on starting structure.

Our original starting method calculations were initiated from the active structure with PDB code 3P0G. To address the editor and reviewers’ comment, we initiated a new string simulation from the inactive state carazolol-bound structure with PDB code 2RH1. We note however that a fully active conformation is observed in the presence of intracellular binding (e.g. G protein), which we do not have in our system and that thus driving the system from the inactivated state to the activated one is expected to be energetically unfavorable and thus more challenging than the reverse. Apart from the difference in conformation, the two receptor structures were identical in terms of side chain protonation and treatment of loops and termini. The initial pathway was different too, as we applied targeted MD on the 2RH1-initiated structure toward the active structure, using the same CVs as before. After convergence, we obtained a similar pathway: calculating the expectation values of the microswitches for the ensemble derived from the system initiated from the inactive state structure, preserves the overall trend for carazolol to cluster away from agonists for transmembrane microswitches (Figure 2 – figures supplement 10).

Interestingly, the pathway is not entirely identical (Figure 2 – figures supplement 9) and some states in the activated region of the landscape, in particular, appear different. The states in the inactive region, on the other hand are equivalent. Since the fully activated state is only stabilized in the presence of an intracellular binding partner, we expect the activated state region to not be sampled fully adequately starting from the inactive state. We would thus recommend initiating simulations from active structures if such structures are available. In addition to this new simulation, we have previously shown that the swarm trajectories near the active state basin were similar to the G protein bound structure with PDB code 3SN6 (Fig. S3 in (Fleetwood et al., 2020), even more so than to the nanobody-bound structure 3P0G, which was used to initialize the system. This suggests that our results are not heavily dependent on the choice of starting structure when it is appropriately chosen.

2. Improved connections to experiments that do not rely solely on Emax values.

Our ambition was to present a quantitative comparison of microswitch simulations to experimental results. Quantitative comparison is notoriously difficult since it typically involves a comparison of physical properties across widely different length- and time-scales. In this study, we thus considered it an achievement that were able to correlate predicted expectation values for microswitches to experimental measurements of downstream signaling. Nevertheless, following the reviewers’ comments, we have complement these findings with a comparison to other types of experiments, especially to biophysical studies which, combined with our results, can provide mechanistic insights.

In the first part of the results section we make a qualitative comparison between the free energy landscapes for the apo, carazolol-bound, and BI-167107-bound receptor to distance distributions in spectroscopy experiments (Manglik et al., 2015). Since the “ionic lock” microswitch is frequently discussed in that and other papers, we have replaced the TM6-TM3 distance in Fig.2 and Figure 2 – figure supplement 9, 11, with a measurement that better reflects the ‘ionic lock distance’. This does not influence the overall results of the paper—the only difference is that we compute the distance from R131 to E268 instead of L266—but it may simplify the interpretation of the results for readers well acquainted with the GPCR literature. We also investigate the change in hydration of TM6’s intracellular part, which could explain an observed change in probe environment.

We also discuss the heterogeneity of states identified by these experiments, the existence of preactive states along the activation pathway and our single equilibrated states. Digging further into the molecular basis for the agonist-specific states, we identify the transmembrane region around the PIF motif to be an allosteric switch, which explains why this region it is a target of allosteric modulators (Liu et al., 2020) (see the final paragraph of the last result section). We also identify signalling hotspot around the NPxxY motif to be the major contributor to the heterogeneity of states, in agreement with spectroscopy studies (Frei et al., 2020; Lamichhane et al., 2020; Liu et al., 2012) (see also second to last result section). Finally, we have elaborated on the general result that our analysis methods tend to detect highly conserved residues in the transmembrane region, which experiments have identified to be critical for the receptor to function (see last result section). All new text is marked in blue, and we hope the reviewers will find it easy to locate these new experimental comparisons in the revised manuscript.

3. Improved relationship to past studies.

The revised manuscript includes more citations of past studies and describes them in greater detail. The introduction has been extended to include a new paragraph that describes previous state-of-the-art MD simultion studies of the β2AR. We also present methodological development from our previous paper (Fleetwood et al., 2020), while introducing changes that make this new manuscript understandable for a reader unfamiliar with our previous work. As described in the previous section, we also put our results into the context of previous experimental results throughout the text. Note however that given the number of studies of this receptor carried out so far, a thorough review of all previous studies is not possible in an original research article.

Regarding the novelty of our work:

One of the reviewer also questioned the novelty of our work. Although GPCR activation has been previously studied with MD simulations and that our work clearly builds on previous efforts, we believe that the results presented here enhance substabtially our understanding of GPCR activation and conformational regulation by ligands. This is made possible by the availability of high quality structures of β2AR which inform our computational work. We also note that this work is also enabled by the sampling procedures we have proposed, which allow the characterization of the conformational activation ensemble at a reasonable computational cost, and thus the repetition of the procedure for a set of ligands.

First, our study shows how intracellular and transmembrane microswitches quantitatively correlate with experimental downstream response. Second, we find that all ligands stabilize distinct pre-active states, but that ligands with similar pharmacological properties stabilize similar states. Although several results agree with previous experimental or computational findings, we provide novel detailed molecular insight into the mechanism and energetics of these states, thus enhancing understanding of GPCR activation. Finally, the manuscript takes a methodological leap toward quantitative data-driven analysis of datasets with comparatively many ligands, and thereby sets a new state-of-the-art for GPCR molecular dynamics analysis.

Reviewer #1:

This manuscript reports advanced MD simulation to describe the allosteric modulation of the β₂ adrenergic receptor, expanding on the authors' recent simulations (Biochemistry 2020, 59:880-891) that explored communication of ligand binding to the G-protein binding site via "microswitches" within the protein. This study uses the same string method to examine the effects of full and biased agonists, antagonists and inverse agonists using previously solved Xray structures. The study also calls on a series of data-based analysis methods to seek answers to the question of how different ligands communicate their changes and affect GPCR activity. Observation of stabilization of active states for agonists and inactive states for other ligands is an important achievement for MD simulation, even if reproducing known experimental results. Also, characterizing how the distribution of kinetically stabilized (active-like) states are controlled by ligands, and correlations between microswitch expectation values and cAMP response experiments are important. However, the manuscript is highly technical and jargon-based such that many readers will not follow the analysis. The statistical reliability and dependence of the converged string solutions on starting conformation do not appear to have been tested. Finally, relationship between this and the previous manuscript (Fleetwood et al., 2020), as well as how results distinguish themselves from knowledge in the field could be made more clear.

We thank the reviewer for this assessment. Regarding the high technical level, we have tried to make simplifications and clarifications throughout the revised manuscript. For example, we now give more background on what the software and machine learning tools do. We have also tried to better define the connection between this manuscript and the previous work. We hope the changes appear satisfactory.

String methods are notorious for depending on starting structure as well as the number and choice of collective variables, and may converge slowly on the optimal path, if at all. The authors do not show plots of convergence (how each PMF changes over #iterations), and do not appear to undertake tests for reproducibility. i.e. Starting from very different initial structures to show that the same final pathways emerge, leading to consistent free energy profiles. For example, starting again with an initial structure for an inverse agonist instead of one for an agonist and show the same results are computed. I note that in the Discussion it is said that several simulation replicas were used to reduce error. But to what does this apply (apparently not to the swarms)? Also, collective variables (CVs) are those defined in the previous study (Fleetwood et al., 2020), but what if the set of CVs is changed or the number of CVs is reduced/increased? While I do not suggest redoing with a different set, how does the reader judge the robustness of the CV-reliant results?

In addition to the arguments listed in the cover letter, we point out that the previously published apo and BI-167107-bound string simulations were performed with a different set of Collective Variables (CVs) and initiated from different pathways (Fleetwood et al., 2020) than the new ligand-receptor simulations presented in this manuscript. Despite the differences in simulation set-up, it is evident from Figure 2 that the free energy landscapes for the agonists adrenaline and salmeterol are more similar to BI-167107 than to the other ligands. The same holds for the apo free energy landscape, which is strikingly similar to the inverse agonists and antagonists. If the starting conditions and CVs alone defined the outcome, this would not be the case.

Regarding the convergence of the PMF, we previously performed cross validation over the string iterations for the TM5 bulge microswitch expectation values (Supplementary Table 2 in the original manuscript). The microswitch expectation values were consistent between the first and second half of the last 224 string iterations, which indicates that these values have converged and were sampled from an equilibrium process. Additionally, samples drawn from an unconverged distribution, would lead to a large standard deviation and/or statistical outliers in the free energy landscapes. Thus our statistical error analysis and convergence estimation provide much strong evidence of convergence.

The method for free energy calculation uses a transition matrix on a CV grid with stationary solution. In the subsection “Free energy estimation”, the authors cite Pan, Sezer and Roux, 2008 (and Fleetwood et al., 2020), but it is not clear that Pan, Sezer and Roux, 2008, did this, and in fact I think the first to do this from a swarms of trajectories solution was Lev et al., 2017 (see also Flood, et al., 2019), not cited here.

The reviewer is indeed correct here. We have updated the citations accordingly in the manuscript.

Having computed free energy projections along collective variables, receptor stability in relation to efficacy instead relies on estimation values for variables and not free energies, and the justifications for this could be improved. The authors write that the width and depth of the basins determine the most probable state projected into that space. This sentence is not clear, but its meaning is important, because ultimately the relative values of free energies are discussed. If the 2D maps are valid projections of the full configurational space (albeit with sampling guided/biased by the string), then Boltzmann integrals over each basin should yield a valid equilibrium constant. I presume the concern is that projections along different pairs of CVs can lead to different apparent free energies, because different CVs map out different proportions of phase space, and potentially envelop multiple states of the system, but does this mean a Boltzmann integral over a site is not a true thermodynamic quantity? In the Discussion, the authors return to add that examination of one projection can overlook what is happening in other coordinates/switches, but more discussion/justification is needed. Related, the authors state that they have "accurately captured the relative stability of states", but this is confusing given the relative stability of states from the maps appears to have been discredited.

We thank the reviewer for this astute set of comments, which we have enjoyed considering and addressing.

To illustrate our point, we highlight Figure 2 and Author response image 2, in which we can certainly say qualitatively that agonists stabilize the active basin. However, due to the issues with projections pointed out by the reviewer and the fact that these energy landscapes are not ideal two-state systems, an integration over state boundaries requires us to, somewhat arbitrarily, define the state boundaries. This will in turn influence any quantitative comparison to experimental values. The expectation values (Figure 2C-F), on the other hand, require no state boundaries, justifying our choice of this methodology.

Author response image 2

Download asset Open asset

That being said, we followed the reviewer’s advice and computed the relative stability of states for the different microswitches by taking the difference in free energy between the active and inactive regions of the free energy landscape (ΔG ∝-lnK, where K is the equilibrium constant). See the results in Figure 2—figure supplement 8 of the revised manuscript. State boundaries that were chosen are depicted as vertical dashed lines in Figure 2—figure supplement 7.The correlation between ΔG and cAMP signaling is strong for transmembrane microswitches, which also had a strong correlation between cAMP and their expectation values. Therefore, we see little practical difference between the two methods, but would recommend using the expectation value because it does not require defining arbitrary state boundaries.

Instead of using free energies, the authors correlate "expectation values" of their CVs to experiments in the bottom of Figure 2. I assume a plot against free energy was abandoned as it was not working as planned? Ideally one would want to see free energies plotted against experimental efficacy, because estimation values of variables such as TM5…, may correlate to efficacy, but not uniquely map to shifting state equilibrium.

Please see our comments in the previous paragraph. We chose expectation values since it does not require us to set state boundaries. This approach provides a single numerical value which, unlike having an entire probability distribution, is easier to correlate to an experimental numerical value. With that in mind, by merely taking the average and ignoring the characteristics of the free energy landscapes, such as the shape and locations of the basins, leaves out potential useful information. For this reason, we included the microswitch free energy landscapes as a basis of discussion together with the correlation plots in Figure 2.

Many analysis tasks were completed by existing packages for which the meanings are not obvious. e.g. Demystifying, Scikit-learn… Materials and methods, even if accepted, deserve a sentence to explain and motivate. e.g. CVs were short inverse inter-residue distances used to train a Restructured Boltzmann Machine (the principles of this machine could be explained simply). The subsection “Data driven analysis reveals that ligands stabilize unique states” is highly technical and not well explained. While many will know about PCA, less will be known about MDS and T-SNE. The ability to identify signaling hotspots in Figure 4 is impressive, leading to a model for allosteric communication. But the machine learning approaches used come across as black box and require better physical interpretation.

We thank the reviewer for pointing this out. It is easy to take one’s own methodological choice for granted! We have made revisions to clarify our rationale and better describe the tools used. For example, we describe what the software demystifying does when we mention it, and try to explain how we find important microswitches by computing the KL divergence in simpler terms.

Reviewer #2:

This manuscript describes the use of enhanced sampling molecular dynamics to calculate free energy landscapes for the β₂ adrenergic receptor. A particular focus is the conformational dynamics and thermodynamics of microswitches that change conformation upon receptor activation. A variety of ligands are investigated to identify shared and divergent features of receptor conformational modulation. Unbiased methods are shown to identify known conformational switches as key regulators of receptor activation.

Overall the manuscript is interesting and clearly written. Figures are very clear and well presented. The subject matter has been very extensively studied in the GPCR family as a whole and in the β₂ receptor in particular, and the results presented here largely align with existing understanding of GPCR activation and conformational regulation by ligands. A major concern is the question of whether the results presented here truly enhance our understanding of GPCR activation, or simply confirm things that are already known. Many groups have published detailed analysis of similar questions to those presented here, including the Dror, Nagarajan, and McCammon groups, among others. There is very little discussion of this prior work, which makes it difficult to see how the present manuscript fits within the broader context of GPCR activation molecular dynamics analysis.

We appreciate that the reviewer found the manuscript well presented. Following the reviewer’s comments, we have extended the discussion of prior work in the manuscript, as also described in the cover letter. Please also find arguments about the novelty of our work in the same place.

A significant technical concern is the reliance on cAMP Emax values as the sole experimental validation. A prospective experimental test of hypotheses generated here would significantly enhance the manuscript. Barring this, a more extensive comparison of computational results with experimental data is essential in my view. The cAMP pathway is highly amplified, which may confound analysis linking Emax values directly to conformational equilibria. Manglik et al., 2015, presents an actual biophysical analysis of conformational equilibria (albeit with fewer ligands) and comparison to this would be helpful. Do relative energy predictions match these experimental observations?

We believe that a comparison to cAMP measurements is relevant and interesting, since it represents a downstream cellular response and a strong correlation suggests our approach to have predictive powers. A comparison to Emax values instead of e.g. EC50 values makes sense, since our simulations have a ligand present at all times, which corresponds to experimental conditions with a very high ligand concentration. We agree that a prospective study would be ideal, but in the absence of experimental validation, we have followed the reviewer’s advice and have improved the comparison to experimental data, and in particular Manglik et al., 2015, as this was suggested by the reviewer. For a more detailed description of what has been done in this regard, refer to our comments under “Essential revisions” above.

A smaller issue is that the utility of the results presented here appears to be a bit overstated. For example, it is claimed (Introduction) that the results provide insight into how ligands with specific efficacy profiles can be designed. To support this statement, the authors should present examples of compounds they have designed based on their results, together with experimental data confirming these molecules have the intended efficacy profiles.

In the revised manuscript we are more precise in what our results allow us to do while avoiding speculative claims. We anticipate to use these methods in prospective studies going forward, but it is indeed not part of this study.

Reviewer #3:

This is an interesting paper, and I like to attempt to connect analysis of allostery to function. However, I'm extremely concerned about statistical uncertainty - it's not really discussed, and it would be easy to chalk all of the results up to limited sampling. It will be important for the authors to demonstrate this isn't the case.

We are glad that the reviewer found our paper interesting and agree that the statistical validation was not clearly presented and also lacking specific tests.

Basically, my concern is that there's essentially no mention of statistical uncertainty or convergence anywhere in the document. One of the major claims is that the different agonists each populate a distinct substate - this is an incredibly important and interesting observation if true, but could also be explained by saying each simulation wandered in its own space and didn't have time to explore anywhere else. If it were run again, it might wander into a totally different place. I don't have a great feel for how rapidly swarms explore, but I do know the total sampling time of 1.4 microseconds per ligand sounds awfully small. There are examples in the literature showing that conventional simulations several times this length are not converged as far the configurations in the ligand binding pocket go (for example, Leioatts et al., Biophysical Journal, 2015, or some of the work from Dror's group).

The reviewer is right, we had omitted important information, relying on information present in our previous paper. We have revised the manuscript such that this manuscript now stands on its own. Specifically, our assessment of statistical reliability can be found in the cover letter.

Note that the single equilibrated states are by definition converged from a starting structure using an adaptive sampling approach with multiple walkers. Sampling a large region of conformational space would require much more than 1.4 microsecond, but a small region kinetically accessible from a starting structure may be well-sampled within this time-frame. Actually, an important methodological aspect of this manuscript is that we could identify ligand-induced pharmacological properties by adaptive sampling of a small portion of conformational space.

I have a couple of ideas for how the authors could convince me this isn't just a sampling artifact. The best would probably be to pick one ligand and redo the whole calculation 4 more times, start to finish - looking at the variation between those replicates would be a decent estimate of the uncertainty. Ideally, they'd do this for all of the ligands, but I recognize that's not a reasonable request, which is why I suggest doing it for 1 ligand.

A weaker test would be to break the individual ligand calculations into blocks, somehow - I'm not totally sure how to do it - and show that the blocks are self-similar. If they're all wandering through the whole of the blob in Figure 3, you're more likely to be ok. If each block populates a discrete chunk of the blob (and the overall swarm data doesn't revisit), then there's a big problem.

We thank the reviewer for the specific suggestions! We indeed followed the reviewer’s second suggestion and performed stratified cross validation with 3 subsets of the walkers. More details can be found under Essential Revisions. We argue that this statistical test is as legitimate as rerunning the simulation for one ligand for two reasons. First, we can perform this test for all ligands, not just one, without wasting computational resources. Second, the method is based on finding a center point in CV space by tracking a number of walkers launched from the same starting structure. Considering the high number of walkers used and that our results pass a cross validation examination, we may conclude that the position of the center point is insensitive to the random seed of the individual MD trajectories. Thus, we expect a new set of walkers launched from the same starting structure to converge to the already sampled state.

I'd want to see error estimates on the free energies in Figure 2, plus a redo of the dimension reduction in Figure 3 to see if the ligands still separate more than the replicates of 1 ligand.

See our comments above.

On a more minor note, the use of p-values in the subsection “Ligands control efficacy by reshaping microswitches’ probability distributions”, isn't really correct. The point the authors are trying to make is that the computed value doesn't predict the experiment (for good reason), and the correlation coefficients make that point for you.

We have removed the use of the p-value and only refer to the correlation between cAMP response and microswitch expectation values.

Author details

1. Oliver Fleetwood

   Science for Life Laboratory, Department of Applied Physics, KTH Royal Institute of Technology, Stockholm, Sweden

   Contribution

   Conceptualization, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-4277-2661

2. Jens Carlsson

   Science for Life Laboratory, Department of Cell and Molecular Biology, Uppsala University, Uppsala, Sweden

   Contribution

   Conceptualization, Supervision, Funding acquisition, Validation, Investigation, Visualization, Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-4623-2977

3. Lucie Delemotte

   Science for Life Laboratory, Department of Applied Physics, KTH Royal Institute of Technology, Stockholm, Sweden

   Contribution

   Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Writing - review and editing

   For correspondence

   lucie.delemotte@scilifelab.se

   Competing interests

   Reviewing editor, eLife

   "This ORCID iD identifies the author of this article:" 0000-0002-0828-3899

• 3,408

  Page views

• 494

  Downloads

• 14

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.