Coevolution-Driven Method for Efficiently Simulating Conformational Changes in Proteins Reveals Molecular Details of Ligand Effects in the β2AR Receptor

With the advent of AI-powered structure prediction, the scientific community is inching closer to solving protein folding. An unresolved enigma, however, is to accurately, reliably, and deterministically predict alternative conformational states that are crucial for the function of, e.g., transporters, receptors, or ion channels where conformational cycling is innately coupled to protein function. Accurately discovering and exploring all conformational states of membrane proteins has been challenging due to the need to retain atomistic detail while enhancing the sampling along interesting degrees of freedom. The challenges include but are not limited to finding which degrees of freedom are relevant, how to accelerate the sampling along them, and then quantifying the populations of each micro- and macrostate. In this work, we present a methodology that finds relevant degrees of freedom by combining evolution and physics through machine learning and apply it to the conformational sampling of the β2 adrenergic receptor. In addition to predicting new conformations that are beyond the training set, we have computed free energy surfaces associated with the protein’s conformational landscape. We then show that the methodology is able to quantitatively predict the effect of an array of ligands on the β2 adrenergic receptor activation through the discovery of new metastable states not present in the training set. Lastly, we also stake out the structural determinants of activation and inactivation pathway signaling through different ligands and compare them to functional experiments to validate our methodology and potentially gain further insights into the activation mechanism of the β2 adrenergic receptor.


■ INTRODUCTION
From receptors to transporters, cycling between conformational states is how many membrane proteins conduct their function.−7 This leaves many potentially functionally relevant conformational states undiscovered, and insights into functional experiments may be misguided.−10 Though successful in predicting folds, the question of which conformational state is predicted is not clear.There is an inherent problem with the lack of balanced data sets for different conformational states in many families, which means that models would be discouraged to predict potentially correct alternative states in favor of the ones that are being trained on. 11,12In addition, transition states and the kinetically accessible paths between alternative states are not captured experimentally and are thus impossible to train on thus far.
To probe the short-time scale (<10 ns) dynamics of conformational change effectively, one can employ molecular dynamics (MD) simulations, where one parametrizes each atomic interaction with a force field, and evaluates Newton's equations of motion in the time domain to simulate a protein structure.Due to the high computational complexity, large systems over long time scales are difficult to simulate in equilibrium MD, a regime in which many important processes such as receptor activation, membrane transport, or ion conduction fall.One approach to prolong the time scales of MD simulations is to enhance the sampling by perturbing the Hamiltonian adaptively based on historical sampling to encourage sampling of new regions and discourage sampling of already visited regions. 13In this work, we utilize the accelerated weight histogram (AWH) algorithm, 14 which adaptively changes a bias potential along a predefined collective variable (CV), which in turn is a function of the atomic coordinates in the system (see the Methods section).Since the added bias of the AWH method is known at all times, it is possible to calculate the convolved free energy surface of a bias-free simulation based on the measured probability distribution.This feature makes it possible to evaluate the energetic stability of the explored conformational space.
Finding an appropriate set of collective variables to sample is crucial for this algorithm to converge and properly sample alternative conformational states.Learning the CV while having no observation of the motion is difficult, and our previous efforts have been mostly centered on a supervised learning approach. 15,16There are however some general principles for designing a good CV: (i) a good CV should describe the conformational change such that all relevant states are separable and uniquely correspond to different points along the CV; 17 (ii) a good CV for AWH or other adaptive biasing methods should contain the degrees of freedom that are responsible for the highest energetic or kinetic barriers. 13n essential hurdle to overcome is that the space of degrees of freedom is very large (3N) and searching all possible configurations would be intangible for very large systems such as proteins or other biomolecules.Thus, it is clear that we need a way to reduce the space of possible configurations by means other than structural information deposited in the PDB.We previously introduced a supervised coevolution-driven methodology to predict conformational states across entire families; however, this required at least one solved structure in each conformational state. 15In this work, we utilize coevolution, how residues that are in contact tend to mutate together to retain the function of the protein, and only one reference structure.By choosing a reference structure and pinpointing which residue pairs are far away in 3D space yet have a high coevolutionary signal.Such residue pairs are in this work referred to as false positive coevolving pairs and could be due to either contacts forming during the folding procedure, longrange interactions, or contacts formed in other functional states.Since evolutionary information is agnostic of whether structural information is available and tightly coupled to function, we wondered whether exploring the space of possible contacts from coevolving residue pairs could aid in discovering new states.We reasoned that if we restricted all possible conformations based on kinetic accessibility and thermodynamics as defined by physical force fields in simulations while enhancing the sampling along CVs containing evolutionary information, we could find new, functionally relevant conformational states.Our expectation of correlation with function is based on coevolving residues being tightly coupled to the function of the protein since an evolutionary pressure must be present to produce statistically significant coevolution signal.By analyzing, the unsatisfied coevolving residue pairs thus provide a convenient dimensionality reduction transformation, whereas the relevant hypervolume to be explored can be determined by an energy-based sampling algorithm such as MD simulations.
To benchmark this methodology, we chose a family of proteins where it is possible to extract coevolutionary information from the multiple sequence alignments and where the conformational landscape is well defined by experimental structures.One of the largest and most pharmaceutically relevant families of membrane proteins are G-protein-coupled receptors (GPCRs). 18,19The energetics of these receptors are finely tuned to respond to a multitude of external ligands, from druglike molecules to lipids or peptides.Thanks to their importance in drug development and signaling pathways, there are a plethora of structures available in various states, making them the perfect candidate on which to validate our methodology.We chose the β2 adrenergic receptor (β2AR) as our focus because of the abundance of both structural and functional information in the form of the deposited PDB structures and activity experiments.As this system has been a subject of numerous computational studies, we considered it an appropriate target for benchmarking and for comparison of the time scales of convergence.We further validated our methodology by simulating the β2AR binding to nine different ligands with experimentally determined activities and structurally resolved binding modes.The resulting free energy landscapes evaluated that our method is able to sample functional conformational states and accurately predict the efficacies of β2AR ligands, providing insights into the molecular mechanism of activation and inactivation and potentially shedding light on the structural basis of biased agonism.

Coevolution Analysis And False Positive Detection.
First, the multiple sequence alignment (MSA) for the GPCR class A family was downloaded from the PFAM database, 20 PFAM ID: PF00001, name: 7TM_1, since this was the family to which the target protein belonged.As to not divert the model's attention from the most conserved regions, we modified the MSA and the training parameters accordingly.Specifically, we filtered positions with more than 20% gaps and assigned an effective weight to each sequence equal to 1/N 0.9 , where N 0.9 is the number of sequences with over 90% identity.Then, the MSA was used as a training set for direct coupling analysis (DCA), with the goal to fit a Potts model in the pseudo-maximum likelihood sense (eq 1) where the conditional probability encodes the position-wise (i, j) information from the entire MSA given the model parameters.The model parameters correspond to column (position)-specific fields v and column pair-specific couplings w.The actual loss function that is optimized is then the natural logarithm of this expression, summed over all sequences n and positions i, and weighted sequence-wise by the N 1 0.9 metric.The partition function Z is a normalization constant such that the probability over position i sums to 1.The loss function was optimized using the Adam optimizer 21 for 300 iterations as implemented in TensorFlow, 22 with a learning rate of 0.0001 and parameters initialized randomly by drawing from a normal distribution of mean 0 and standard deviation 1.The code was adapted in the described manner from the original GREMLIN implementation. 23he pairwise coupling parameter w is in the form of a 20 × 20 matrix, which is then summed over all axes to yield a scalar value of the coevolution strength for a certain (i,j) pair.To The Journal of Physical Chemistry B reduce the row-wise noise, an average product correction (APC) correction is applied so that comparable scores are attained.
The coevolution scores are then overlaid with the physical contact map of the target protein's starting structure, and the false positives are identified and clustered based on coevolution strength (see the Exploration Simulation section for details).We have made this script available on the attached Zenodo project which can be found by following this link 10.5281/ zenodo.8164731.
MD Simulations.All Molecular Dynamics (MD) simulations were carried out in GROMACS2021.The simulation systems containing the β2AR in the active state (PDB ID 3P0G) bound to different ligands were prepared using the CHARMM-GUI membrane builder, 24 and the ligands were parametrized using the built-in cgenff 25 procedure available in CHARMM-GUI.Two histidines (H172 4.64 and H178 4.70 ) were protonated at their epsilon positions, and the rest of the histidines were left protonated at the delta position.Furthermore, a Na + ion was placed at the sodium binding site 2.50 and left unprotonated, as done in earlier work. 26The systems contained the protein embedded in a POPC bilayer plunged in a 0.1 M KCl solution.The initial periodic boundary condition (PBC) box was 85 × 85 × 94 Å 3 , ensuring at least 12.5 Å of water molecules between the protein and the PBC box end at least 10 lipid molecules between each PBC copy of the protein.The force field used was CHARMM36m 27 for protein and lipids, and TIP3P for water.The models were equilibrated using the default CHARMM-GUI scheme with one minimization step and 6 100 ps restraint cycles with gradually released restraints in the NPT ensemble, followed by a production simulation of 10 ns.The simulations were carried out by using a 2 fs time step.The target temperature and pressure were set to 303.15 K and 1 bar, respectively, and maintained by a Nose-Hoover thermostat (coupling separately protein, lipids, and solvent) and the recently implemented Crescale barostat with semi-isotropic coupling (p = 5.0, compressibility 4.5 × 10 −5 ).Hydrogen bonds were constrained by using the linear constraint solver (LINCS), and long-range electrostatics were accounted for by using the particle mesh Ewald (PME) method beyond the 12 Å electrostatic cutoff.A neighbor-list cutoff was used for vdw interactions with r vdw = 12 Å and a switching function starting at 10 Å.
Exploration Simulations.The nonequilibrium exploration simulations for the apo system used the same parameters as described above, but the Hamiltonian of the system was further modified after equilibration.The false positives identified by aligning the coevolution map and the contact map calculated from the structure were clustered into three clusters based on coevolution strength with the K-means algorithm, and the minimum-distance atom pairs were used to define individual distances by the pull code in GROMACS.Then, three transformation pull coordinates corresponding to the three clusters (determined by the mean shift algorithm) expressed the mean of all intracluster distances, respectively, for each of the three clusters.Then, umbrella potentials were applied to each transformation pull coordinate in four independent replicas, in order to push the system toward an alternative state.The force constant used was 200 kJ/mol nm 2 per contact, and the target was 0.2 nm.After 10 ns of pulling, the added potential was switched off and an additional 10 ns of unrestrained MD simulations was run.The resulting trajectories were then used to train an SVM classifier with a linear kernel to distinguish the ensembles from the equilibrium simulations based on the reference structure 3P0G and end state of the nonequilibrium pulling simulations.The top 20 contacts were then used to define two CVs through a weighted sum: one with the positive and one with the negative coefficients.Ultimately, these were then used to define the CV space for the AWH simulations.The exploration simulations were not repeated for each ligand system, but rather the CV generated for the apo system was used for all ligand-bound systems, ensuring that a common projection in a set of low-dimensional projections.The necessary files to reproduce these simulations can be found in the associated Zenodo project (see the following link: 10.5281/zenodo.8164731).
AWH Simulations.The AWH simulations were run with 4 walkers, two starting from the equilibrated reference structure based on 3P0G and two starting from the end state as explored by the nonequilibrium pulling simulations.The walkers move independently but share a bias and update frequency to the same histogram.To increase the overlap of walkers in the early stages of the simulation, a point is considered covered only once a walker has sampled all points within the cover diameter (0.4 nm) of a point.Additionally, the Wang−Landau algorithm 28 in the final stage is only initiated once the total histogram has been equilibrated to reduce wildly different responses from each walker to the same walker.All necessary input files including tpr files are available at the Zenodo project (10.5281/zenodo.8164731)As the explored state may be far from equilibrium, we set the growth factor parameter to 2.0, which helped prolong the number of initial sweeps required for the simulations to initiate the Wang−Landau algorithm.As mentioned in our previous work, we have experienced that free energy estimates stemming from this modified parameter produce long-term reliability and flatter distributions without walkers getting stuck.Furthermore, while the cover diameter, energy cutoff, number of steps between each AWH sample, and number of samples between updates were, respectively, set to 0.4 nm, 120 kJ/mol, 10, and 10 for all simulations, parameters for individual simulations were fine-tuned to allow for convergence on faster time scales (see Table 1).
The CV spaces were all defined in the ranges of 0.4−1.3 in CV 1 , and 0.3−1.0 in CV 2 .
Network Analysis.Collective variables need to reduce the dimensionality of the measured system for interpretability and, in this case, computational efficiency in converging a free energy landscape in the same space.This brings with it the risk that important degrees of freedom for moving between states The Journal of Physical Chemistry B may not be included in the CV description.Thus, after convergence, we prolonged the simulations (with a nearly static bias) another 100 ns to achieve sufficient sampling for degrees of freedom not directly incorporated in the CV.How to exactly determine what "sufficient" sampling means is difficult to estimate even a posteriori, but we estimated that we would probably need as many coverings as would be required for convergence (ca.20−50).Based on this sampling, we estimated the frame-wise free energy as a function of the collective variable (eq 2).
= w t u i t X ( ) e / ( , , ) where U(X(t)) is the point-wise free energy estimate, RT is the thermodynamic constant at 298 K, and u(i,t,X) reflects the binning procedure which, when summed over, gives the number of frames in the bin that is visited at time t.Moreover, X(t) represents the value of the reaction coordinate at time t, simply taken from the simulation frame.Given a new reaction coordinate Y, we could calculate the bin-wise energies according to eq 3.
While any observable Y i could be used in this formulation, we projected the free energy surface onto every minimum distance between residue pairs in the protein−ligand system.We then extracted the difference between the lowest free energy basin and highest free energy barrier to assess if there was tight or loose coupling between any alternative conformational state of each residue pair.We then used this measure as a measure of energetic coupling, which we used to construct a network via an adjacency matrix.We identified a substrate-coordinating residue on the extracellular side of the receptor, F193 as a source, and determined the most efficient coupling pathways to three sinks (P330, R131, and E268) on the intracellular side using Dijkstra's algorithm 29 as implemented in the python package Networkx. 30To get comparable estimates of a residue's importance in this network, we calculated the betweenness centrality, the number of shortest paths passing through each residue.All scripts utilized in the analysis are available in the associated Zenodo project (10.5281/zenodo.8164731).

Quantification of Populations.
To quantify the fraction of active populations, the populations of each basin were estimated using thermodynamic weights estimated for each bin of the free energy surfaces.The basins were defined as centered in the minima and continuing up to the border of inflection points, as determined by the InfleCS methodology. 31Additionally, to reduce the noise from microstates, basins composed of less than 3 bins were discarded.Moreover, the functional state of the basins (resting or active) was assigned through the following procedure: (1) The frames within each basin but after convergence was reached were extracted.(2)  The snapshots were aligned with a reference active structure with PDB ID 4LDE.(3) The RMSD between the TM6 helix of each frame and the active reference structure of TM6 was calculated, where TM6 was defined as residues 269−298.The classification was such that if the average RMSD was below 3.0 Å, the basin was labeled as active, and vice versa for the resting label (see Table 2) In addition to simple resting/active labels, we also labeled basins based on their position in the CV space, with generally

The Journal of Physical Chemistry B
the upper right being more resting and lower left more active.With the observation of multiple potentially active basins, we also introduced the A 1 and A 2 states in the lower left and lower right corners of the CV space, respectively.For the calculation of the fraction of A 1 -like population in the salmeterol and salbutamol simulations, we defined the A 1 state as defined by the borders of the A 1 basin in the adrenaline-bound simulations.
Error Analysis.Due to the nature of AWH simulations with an adaptive bias that changes for every iteration based on prior sampling, it is difficult to apply statistical analysis to gauge the point-wise error in the free energy estimations.Our strategy for quantifying the error was to measure the deviation from a flat coordinate distribution when the bias converges.This approach is sound because the negative bias function − U(X(r)) should approximate the underlying Hamiltonian H(X(r)), producing an effective Hamiltonian that corresponds to a flat coordinate distribution (eq 4).
where the effective coordinate distribution at equilibrium would be P eff .After the histogram equilibrated (meaning that the mean of the sampled distribution should be within 80% of the target distribution), we calculated the transition matrices between neighboring bins, which tell us about the transition imbalance in different regions in the free energy landscape.For an optimally converged system, all transitions should be roughly equally probable, and the higher the imbalance the higher the free energy estimate error is in this region.This can The Journal of Physical Chemistry B in turn be translated into free energy deviations from a flat distribution.
When the free energy estimation is based on several different walkers, it is important to check whether the sampling is consistent between them.We used the overlap metric, which is defined as the number of walkers that have a coordinate distribution of more than 10% of the mean coordinate distribution.This is calculated for all points and gives an estimation of how many walkers sample a flat distribution according to the shared bias function.
■ RESULTS AND DISCUSSION Computational Pipeline for Structure Prediction and Free Energy Landscape Computation.For the methodology to be applicable, we need a starting structure (experimentally resolved or obtained thanks to structure prediction algorithms) and a deep MSA.We initiated this work from the active structure of the β2AR (PDB ID 3P0G) using an MSA of the GPCR class A family as a basis (Figure 1A) and analyzed residue pair couplings using direct coupling analysis (DCA). 32We then compared the coevolution coupling scores inferred from the MSA with the contacts found in the structure (Figure 1B), where the coevolving residue pairs not in contact in the structure (false positives) are clustered according to their coupling scores using hierarchical clustering, with the number of clusters determined using the mean shift algorithm (Figure 1C).The resulting contact clusters were then aggregated into collective variables, as weighted sums where the weights correspond to the coupling scores (and where a CV value of zero corresponds to all contacts being formed) and then used as coordinates for steered molecular dynamics to enhance the sampling along the chosen dimensions.The ensemble explored with steered MD simulations (Figure 1D) is then analyzed using a machine learning model to extract collective variables describing the conformational ensemble obtained thus far.A good CV to describe the conformational landscape should encompass the degrees of freedom with the highest intrinsic energy barriers.We thus attempt to find these degrees of freedom by finding those that are maximally separated in the maximum-margin sense between two states.Maximum-margin separation would result in a steep potential function and is thus a good approximation for finding these degrees of freedom.In this work, we used a support vector machine (SVM) with a linear kernel. 33The data set is based on the tracked coevolving residue pair distances, both false and true positives, corresponding to contacts that form and break during the conformational change, respectively (Figure 1E).This CV was then divided into positive and negative coefficients and used for AWH simulations in order to ensure a monotonic response to increasing distances in the simulations, which we have observed produces faster diffusion rates across the CV space.In principle, should the exploration not be sufficient (identified by, e.g., the lack of other basins or structural features missing), the already explored degrees of freedom could be downweighed and another round of exploration could be made until sufficient conformational exploration is seen.In this work, one iteration was sufficient.
In order to probe the conformational range of the β2AR receptor, a CV space was trained, as described in Figure 1.Starting from the active state based on the PDB structure 3P0G, the discovered end state after steered MD simulations (Figure 1D) features the hallmark TM6 inward movement as is expected for any GPCR deactivation. 1 This does not, however, guarantee that the intricate movement of all of the conserved microswitches is captured.To assess whether all multidimensional aspects of the inactivation mechanism were modeled, free energy surfaces were calculated in the aforementioned final CV space using the accelerated weight histogram (AWH) method 14 (Figure 3), natively implemented in GROMACS 2022, running with 4 parallel bias-sharing walkers. 34Using the test of the flat probability distribution described in more detail in the Methods section, we found that the free energy surface converged in relatively fast time scales (∼50 ns, Figure S1).In most cases, a rough indication of whether the active or inactive states are favored is evident as early as 10−20 ns per walker.As discussed in our previous work on the sugar transporter family, 15 it seems that combining coevolving residues with contacts observed in structural ensembles of different states indeed effectively captures the degrees of freedom with the highest free energy barriers, and thus enables to quickly sample the conformational landscape with adaptive biasing.
To summarize, the methodology requires a starting structure and a deep enough MSA to technically enable conformational exploration.Additionally, assessing the number of cycles required for sufficient exploration can be done in several ways, from completely system-agnostic by finding new basins in the free energy landscape or completely system-specific by measuring functionally important measures.It should be noted that assessing the functional relevance of explored conformations merely by measuring a selected subset of structural features is not enough, and alternative state conformations that are more stable than the original starting structure should be seen to rule out hysteresis effects (see Energetical analysis).To effectively pursue such a strategy, one usually requires additional information about state-stabilizing conditions, such as pH, bound ligands, etc.
Conformational Analysis.To assess and validate the accuracy of the methodology, we first aimed to characterize the accuracy of the generated structural ensemble.We did this by extracting frames from MD trajectories corresponding to the main energetic basins.The basins that appeared on the free energy surfaces were labeled according to the TM6 position (residues 266−298) of the configurations making up the free energy minima.The latter was defined by a C-α RMSD metric to the reference (starting) position of the TM6 helix.Thus, we could classify the basins into inactive (above 3 Å RMSD of the reference TM6 position), and active (below 3 Å RMSD of the reference TM6 position) states (Table 2).Additionally, two separate basins that could both be classified as active were separately labeled as A 1 and A 2 based on their positions in CV space (see Methods section).Since the generative procedure started with a reference structure in the active state, we assessed the accuracy of the exploration by comparing the C-α RMSD of the ensembles to two other reference structures (PDB ID 4LDE 35 as an active reference, and PDB ID 2RH1 36 as an inactive reference).The ensembles investigated were (1) equilibrium simulations of the starting structure based on 3P0G, (2) equilibrium simulations after nonequilibrium pulling, and (3) the ensemble of snapshots taken from the resting basin after the AWH simulations had converged (Figure 2A).The RMSD toward the active state is increased by the end of the exploration simulation, but only slightly decreased toward the inactive reference; however, the exploration is enough for an appropriate CV to be learned such that the correct inactive ensemble (<1.5 Å RMSD) is The Journal of Physical Chemistry B indeed sampled during the AWH simulations.As one might expect, the nonequilibrium nature of the exploration simulations gives rise to unphysical artifacts such as helices breaking or poor side chain stacking, but this is consistently resolved on the time scales of the AWH simulations.An important result is that the basin states from the AWH simulations overlapped with the experimentally determined resting states of the receptor (Figure 2A), confirming that we can indeed find alternative states in the conformational cycle of a receptor by using this method.
As the entropy of a dissociated TM6 helix in the active conformation is higher, we divided the active population into A1 and A2 states based on their position within the CV space (Figure 3 and Table 2).Upon visual inspection of the conformational ensembles in these basins, it became evident that they corresponded to alternative conformations of helix TM6 adopting different radial positions at similar TM6-TM3 distances (Figure 2B).We found that the angle between the TM6 helix conformations was on average 16°radially and that the hinge position of this movement was located at G280 6.42 where the TM6 helix also exhibited a slight bulge.The movement of TM6 was also observed to be accompanied by a similar hinge-driven motion in TM4 that starts with a twist at T164.Additionally, different ligands stabilized the A1 and A2 states differentially (Figure 3), which may hint at the remnants of a molecular mechanism responsible for biased agonism.Interestingly, both active states are separable in this CV space even though the reference structure (4LDE) more closely resembled the A1 basin, showing that alternate local states can also be explored and sampled with this methodology.
Energetical Analysis.Even though the correct conformational space may have been searched, as indicated previously (Figure 2A), we aimed to evaluate the quality of our free energy landscapes by comparing the fraction of the β2AR population that resides in the active state.Thanks to the availability of a multitude of functional experiments measuring the efficacy of various ligands, 37−43 we could correlate our values to the %Emax values in relation to adrenaline, which was an available data point in all studies used.Before we could analyze the free energy surfaces, we ensured that convergence had been achieved.Thanks to the monotonic growth of the reference weight histogram in AWH simulations, the free energy update size is guaranteed to decrease for every iteration, which we have demonstrated occurs in our simulations (Figure S1).To confirm that the correct solution is found, we would expect that the effective Hamiltonian is close to uniform, which in turn should give a flat coordinate distribution and an exchange of walkers (Figure S2).The expectation of a flat coordinate distribution also enables analysis of the approximate deviation in kJ/mol.For details of the utilized metrics and algorithms, see the Methods section.
First, we needed a global definition for assigning a functional state to each basin.Even though the CVs acted on the same atoms in each condition, the surfaces were not completely superimposable in terms of where the basins fell in the CV space.Importantly, the resting and active states are separable under each of the simulation conditions even though they may be shifted.Interestingly, this suggests that there are slight differential stabilizations of molecular conformations that still correspond to similar functional states classified by the TM6 position.This may be due to shifts in the atomic coordinates as a response to each ligand, since chemically similar ligands such as salmeterol and salbutamol produced superimposable results.This result, however, meant that any analysis of the  The Journal of Physical Chemistry B populations had to be based on basin population and could not be based on specific CV values.Furthermore, the basins needed to be labeled as either resting or active, after which the populations could be summed over all active states.We based this on the RMSD of the TM6 position toward the previously used reference active structure with PDB ID 4LDE (Table 2).
Visually one can clearly distinguish antagonists and inverse agonists (apo, carazolol, alprenolol, and timolol) from full agonists (formoterol, adrenaline, salbutamol, salmeterol, and BI-167107), where the most energetically favorable state is either resting or active, respectively.Without any further quantitative analysis, it is possible to distinguish ligands such as alprenolol and mirabegron, which cover a conformational space similar to the apo simulation, and carazolol and timolol, whose basins are shifted away from any active state.Thus, it is possible to qualitatively distinguish the behavior of the receptors bound to these different ligands.Our free energy landscapes directly show the existence of multiple distinct conformations being favored or disfavored, contrary to a simple two-state model.This agrees with modern observations of the existence of many functionally relevant and accessible conformations in the form of metastable states. 44owever, to test the hypothesis that the CVs generated by the unsupervised conformational exploration were able to capture the intricacies of the β2AR receptor activation mechanism, we quantified the fraction of the β2AR population in all active states for each condition (Figure 4), with the different activation states being defined as all bins within the inflection border of each local minimum.Additionally, we plotted these active fractions against experimentally measured relative Emax values (using adrenaline Emax levels as a reference).This ensured experiments could be compared to one another.
Additionally, in the cases where no experimental values could be found in the literature (BI-167107 and Carazolol), we used the mean value obtained from our previous work on this receptor using string simulations. 45As mentioned before, salbutamol and salmeterol exhibited uniquely different stabilization patterns of the A1 and A2 states, which we hypothesized could be related to other molecular mechanisms, such as biased agonism.The existence of an orthogonal active state unique when bound to at least salbutamol has previously been reported in accelerated MD simulations by Tikhonova et al. 46 To probe this hypothesis, we compared the fraction of active population based on all active states (green squares, Figure 4) and only the A1 state, as defined by the adrenaline simulation (purple pentagons, Figure 4).
The quantification analysis shows that, in addition to displaying differentially stable resting and active states comparing full and inverse agonists, the relative populations of said basins are consistent with relative Emax values, including those determined for partial agonists (Figure 4).In contrast to previous work on similar systems, we were able to estimate the populations directly, instead of inferring dynamics from structural observables, such as the TM6 angle or microswitch positions.This suggests that the CV trained on coevolution-driven conformational exploration can indeed capture the intricacies of the microswitches that govern the ligand response from the receptor.Thanks to the correlation based on the fraction of active populations, we can conclude that ligands such as salmeterol completely activate the receptor by skewing the entire population to the active state.This  4. Correlation between experimental (Emax values relative to adrenaline along the x-axis) 37−43 and computational activities (calculated fractions of the active population basins along the y-axis).The correlation line was fitted with least-squares, and the R-value was calculated based on this fitted line.Due to differing efficiency values in assays measuring G-protein and β-arrestin signaling of salmeterol and salbutamol, we calculated both the fraction of active population of all active states (green, top right) and only the A 1 state as defined by the basin seen in the adrenaline simulations (purple pentagons, middle).The error bars along x are taken from the references, and the error bars along y are calculated based on the population difference that would arise from the calculated free energy landscape errors.* = Fraction of active population calculated based on the population of the A1 basin, where the A 1 basin is defined by where it is located in the adrenaline-bound simulations.These points were not included in the calculation of the correlation factor R ** = Assumed computational values averaged between Ichikawa et al. (2016) 42 and Fleetwood et al. (2020). 45These points were not included in the calculation of the correlation factor R.

The Journal of Physical Chemistry B
suggests that the salmeterol activation of the β2AR is the theoretically highest activation response one could obtain.
In 19 F NMR spectroscopy experiments, Kim et al. suggested the β 2 AR enthalpy changes could stem from the basal activity of β2AR 47 due to weaker interactions stabilizing the resting state than in the active state. 48A similar conclusion was reached using orthogonal 19 F NMR experiments as well. 49The existence of basal activity is probably due to a small fraction of the receptor residing in the active state, which according to our analysis should be around 38% (Figure 4).In the same study, Kim et al. found that "two inactive states accounted for roughly 60% of the total spectral intensity for the apo and inverse agonist saturated samples", 48 which quantifiably agrees with our estimate of 100−38 = 62%.
Most experiments utilized the cAMP Emax value as a percentage of adrenaline activation of the G-protein pathway, while others additionally probed the β-arrestin pathway and quantified biased signaling. 38Interestingly, while most agonists stabilized active-like A 1 and A 2 basins, salmeterol and salbutamol exhibited the strongest shift toward one of these two basins.Curiously, these ligands have been reported as biased agonists in different experiments, where in cAMP-based assays, the activity is higher than that of adrenaline, whereas when probing the β-arrestin pathway 43 or in FRET-based experiments 41 they are both classified as partial agonists.Since salbutamol and salmeterol both exhibited similar and unique behavior in that the A 2 basin was significantly more stabilized than the A 1 basin, we wondered whether the unique behavior could be related to the biased signaling observed experimentally.We thus set out to separately quantify the relative A 1 population for these two ligands to see if we could relate the two different basins (Figure 2B) with biased agonism.While the relative population of the A 1 basin matches the experimental Emax values obtained when probing the βarrestin pathway (Figure 4, purple pentagons), it should be mentioned that this is highly dependent on the definition of the functionally relevant A 1 state.Nevertheless, we can still qualitatively identify differentially stabilized basins and (B) The β2AR apo system is colored by the per-residue relative betweenness centrality, with three main allosteric paths from the ligand-binding pocket (F193 45.52 ) to the G-protein binding site (E268 6.30 , R131 3.50 , P330 8.48 ) highlighted as sticks.Analogously plotted are the β2AR carazolol system (C), the β2AR timolol system (D), the β2AR adrenaline system (E), and the β2AR salmeterol system (F).The shortest paths shown in this figure have been made available for all ligand-bound systems in Figure S4.

The Journal of Physical Chemistry B
reasonably explain biased agonism.Interestingly, this correlation leads to interesting mechanistic implications of the A 1 and A 2 structural snapshots of Figure 2B, suggesting that biased signaling is dependent on the local rearrangement of the TM4, TM1, and TM6 helices.Our observation is consistent with the model proposed by Liu et al. based on evidence for ligands favoring one of two distinct conformations. 50Additionally, Zhao et al. observed that different conformers of the receptor are necessary for stable binding of Gs or β -arrestin, which have different binding modes. 51In particular, Zhao et al. observed that ICL3 was less stabilized during β2AR-β-arrestin 1 simulations than in β2AR-Gs simulations due to a salt bridge between R239 5.77 of β2AR and D343 of Gs.In our simulations we observed that in the A 1 basin, the TM4 and ICL3 indeed adopt a more closed conformation than in the A 2 basin (Figure 2B), suggesting that population dynamics of the receptor directly impairs not only signaling but also binding of downstream pathway proteins.Our methodology could potentially help in identifying alternative conformational states not necessarily on the main functional pathway.However, to firmly establish the functional relevance of the A 1 and A 2 states, further studies on the preferential binding of the G-protein and β-arrestin will be required.
Allosteric Signaling Pathway Analysis.Given that the CV was able to capture the activation mechanism of the β2AR, we aimed to uncover the details of the molecular mechanism and which molecular determinants were responsible for the differential responses by different ligands.Therefore, we analyzed the structural characteristics of the entire conformational ensembles under each condition by computing interresidue minimum distances, including interactions with any ligand along the 4 × 200−400 ns long trajectories.Then, we projected the free energy surfaces of Figure 3 onto each interresidue distance (see the Methods section for details).The methodology yielded the relative energetic coupling (as defined by the largest barrier between local minima along the distance) between all residue pairs under each condition.By performing a principal component analysis (PCA) on the energetical couplings for each condition, we were able to cluster the ligand-specific molecular mechanisms of activation (Figure 5A).The first two components explained 74% of the variance in the data set and were able to cluster ligands with different effects together.Naturally, the apo simulations were unique, given that no ligand was present.Nevertheless, the closest resemblances to the coupling patterns seen in the apo simulations were the inverse agonist and antagonist simulations.Two distinct clusters were found for agonist pathways, containing both full and partial agonists, suggesting two modes of signaling through the receptor.
In order to further analyze exactly how these coupling patterns were different, we then constructed weighted graphs representing each residue as a node and each edge as the energetic coupling score mentioned above.We then analyzed these graphs by evaluating per-residue betweenness centrality and computing the shortest path from the extracellular to the intracellular faces of the graphs.Combining the paths and the PCA clusters, we were able to identify four main avenues of allosteric communication (Figure 5B−E).The unmodulated energetics and pathways of the receptor can be seen in the apo state.The resulting shortest pathways of the ligandwise network analyses of all simulated systems can be seen in the Supporting Information (Figure S4).
Evidently, for a signal to be transduced from the orthosteric binding pocket (defined as F193 45.52 ) to the intracellular binding pocket of the receptor (defined as three sinks at P330 8.48 , R131 3.50 , and E268 6.30 ), several different pathways are exploited.While the TM6 pathway is activated in almost all cases for communication to E268 6.30 , there are alternative pathways that characteristically correspond to agonistic or antagonistic behavior.In the absence of ligand, anchoring residues around the binding site such as N293 6.55 , N312 7.38 , and Y316 7.42 are responsible for passing the signal along to other helices in the receptor (Figure 5B).The vertical signal also passes through D79 2.50 which has a significantly higher betweenness centrality, along with Y316 7.42 and Y326 7.53 .Indeed, D79 2.50 acts as a switch to bridge the signal between the two tyrosines, which sheds light on the mechanism of deactivation once the ligand has left the activated receptor.This unique mechanism of deactivation in the apo state, as seen from the isolated position in PCA space (Figure 5A) suggests that the receptor activates and deactivates asymmetrically.
Unsurprisingly, the closest behavior to the apo system is found in cluster 1 containing the partial inverse agonist carazolol and antagonist alprenolol.While D79 2.50 still has a high betweenness centrality score, the communication with Y326 7.42 occurs via TM6 as opposed to through TM2.The TM6 pathway is conserved, and the rigid communication through the backbone seems to be inherent to favoring the resting state, which also agrees with the less mobile conformation of TM6 in this state.In accordance with the lock-and-key model of receptor activation, the ligand-binding residues of cluster 1 have the lowest betweenness centrality of all ligand-bound simulations.(Figure 5C) Instead, they seem to block the communication to TM7 and TM2 by not responding to the activation signal.Upon inactivation, communication between TM6 and TM3 is made possible by closing of the intracellular G-protein-binding cavity.The communication is mainly facilitated by M279 6.41 ̧L272 6. 34 , and L275 6.37 is consistent with the observation from Dror et al. when simulating the inactivation process, suggesting that a similar movement and arrangement of side chains is achieved in our simulations as well. 38n both clusters 2 and 3, the signal is preferentially transduced via the ligand, indicating that the ligand has a larger energetic effect compared to the unmodulated free energy landscape.Interestingly, BI-167107 and adrenaline, two agonists, fall in the same cluster as the most extreme inverse agonist timolol (Figure 5A).This could suggest that the energetics are communicated through similar pathways but the energy balance is shifted to favor different states.When analyzing the details of the communication pathways in cluster 2, it becomes clear that there is a preferential signaling via TM3, in particular through a connector region centered around I121 3.40 and L124 3.43 .Additionally, S207 5.46 becomes a central ligand-binding site residue that transfers the signal to different pathways in all agonist-bound systems, a residue that has been indicated in the related β1AR as a determinant of agonism 52 (Figure 5D,E).These regions with strongly statedependent interaction patterns have been previously observed in equilibrium MD simulations by Dror et al., where the hydrophobic connector region of I121 3.40 and S207 5.461 comprised two spatially distinct clusters of residues that moved together. 38In both the adrenaline and BI-167107 cases, the pathways split at the S207 5.461 residue and are transduced The Journal of Physical Chemistry B within TM3 and TM7.Additional dynamics between N322 7.49 , S120 3.39 , and Y326 7.53 facilitates signaling through this pathway rather than utilizing only the TM3 pathway as in the timololbound system (Figure 5E).This interaction pattern corresponds to the hallmark NPxxY motif, which has been shown to be a crucial microswitch, as well as functionally important and mutationally sensitive. 1,53urthermore, a clear difference between clusters 2 and 1 is the absence of the TM6 pathway.In an NMR study where, among others, M279 6.41 was probed for conformational change in the presence of an array of ligands, Kofuku et al. measured a clear shift in carazolol, but no signal in the presence of agonists. 49,54This is consistent with our observation that the pathway is utilized only in cluster 1 and apo simulations and is passed over in other simulations.
Finally, the third cluster consisting of salmeterol, salbutamol, mirabegron, and formoterol exhibits more similar allosteric behavior to each other than the aforementioned clusters.In particular, there are generally two points of contact with the F193 45.52 residue that transduce the signal through TM6 and TM3 to TM7.Additionally, Y326 7.53 , which has been shown by deep mutational scanning to be particularly sensitive to mutations, 55 is transducing the signal via the aforementioned connector region centered around L124 3.43 to I127 3.46 .This showcases how single residues and motifs could take on multiple roles depending on the external stimuli.From an evolutionary perspective, single residues taking on multiple roles could aid multiple interactions to take place to evolve fine-tuned receptor mechanisms and provide mutational stability as a sort of safety net.
From a mechanistic perspective, the strongest interactions under a given condition would dominate the communication in the allosteric network.As expected, we find that most interactions that have been shown to stabilize the inactive state are favored as pathways under inactivating conditions (when bound to antagonists) and vice versa under activating conditions (when bound to agonists).To summarize our findings, we see that ionic locks R131 3.50 and E268 6.3019 are utilized primarily in the antagonist-and inverse agonist-bound simulations.Moreover, internal signaling to the ionic lock region through TM3 via D130 3.49 is preferred when bound to full agonists, which is also consistent with the local change in side chain rearrangements in both active and G-protein bound structures. 19A second highly conserved motif, the NPxxY is expected to form interactions through Y326 7.53 across the Gprotein binding pocket in the inactive state, chiefly the backbone of L124 3.44 and the side chain of Y219 5.58 . 19These cross-pocket interactions are only seen in the antagonist cases and notably absent in the case of the inverse agonist timolol, which we have already described disrupted these interactions allosterically.Continuing up the receptor, water-mediated network has been suggested to involve D79 2.50 to transfer the allosteric signal from the orthosteric site to the G-protein binding pocket in the inactive state, which is consistent with our findings only in the apo state, while other inactive statestabilizing conditions do not exhibit the same pathway, even though the residue lights up as it has a high betweenness centrality.Finally, at the orthosteric site, the most common ligand-binding residues such as S207 5.46 , N312 7.39 , I121 3.40 , and F193 45.52 are all found as stabilizing interactions in a multitude of structures, suggesting that the interaction networks start from these residues. 19CONCLUSIONS Predicting alternative conformational states of membrane proteins is inherently coupled to comparing conformers that are energetically fine-tuned to exist under physiological conditions.This is further exacerbated by the fact that even a light external bias from an experimental structure determination may completely shift the delicate balance of populations, leading to complications in relying on the PDB database. 11In this work, we present a computational methodology to simultaneously explore and enhance sampling of the conformational change of the β2AR receptor.By accelerating the sampling with the AWH algorithm, we also enabled an on-the-fly free energy landscape estimation, providing the crucial thermodynamic validation of newly discovered states.We found excellent structural agreement between the discovered conformations and available experimental structures.The fluctuations were comparable with typical equilibrium simulations of crystal structures (∼1.5 Å RMSD), suggesting that we indeed reached the resting state by simulations initiated in the active state.
Upon quantifying the fraction of active population under binding of different ligands, the stunning agreement between experimental and computationally determined percentages of activation relative to adrenaline showed that we did not only describe the appropriate global rearrangements but also achieve a detailed exploration of individual microswitch conformations.We found further evidence to support the model of a set of loosely coupled switches and an ensemble of pathways that fine-tune the behavior of the receptor, rather than a single cascade that shifts the population of the receptor. 38dditionally, the unsupervised nature of the conformational exploration enabled not only the discovery of alternative conformational states but also revealed structures that plausibly relate to biased agonism identified in the salbutamol and salmeterol bound simulations.Through this observation, we could assign the functional significance of each active basin, bridging the gap between the structure and function in activated ensembles.This observation stresses the importance of analyzing full ensembles of potentially undiscovered areas of the conformational space and possibly sheds light on the discrepancy between experiments in the cases of salbutamol and salmeterol.
By combining the recent surge of computationally predicted structures in a single conformational state (using algorithms such as Alphafold2, Omegafold, or trRosetta), 8−10 with coevolution-driven conformational exploration, it would be possible to shed light on previously unvisited corners of the protein universe.Coevolving pairs are not only prevalent in alternative conformational states but also may reduce the necessary search space for protein−protein interactions or folding pathway studies using expanded ensemble MD simulations.In principle, one could repeat the exploration many times to explore and seed new states.One may use the resulting free energy landscapes to assess the energetics of the resulting conformations.Should no new states arise, then one should continue the cycle.In this paper, only the CV estimation after the first iteration of conformational exploration by nonequilibrium pulling was enough for both quick and accurate convergence.
To conclude, the fast convergence time of this work (ca.50−100 ns per walker, i.e., 0.2−0.4μs in total simulation time The Journal of Physical Chemistry B per system) is around 10 times faster than some previous attempts at simulating the conformational change in similar GPCR systems, 56−58 and indeed this exact system (4.3 μs total simulation time per system with the string method). 45This study thus provides a powerful proof of concept of combining evolutionary data and physics through machine learning to accelerate the exploration and sampling of the conformational space to discover new conformational states.Since the methodology is not tied to this specific system or the GPCR family, this could in principle be applied to any GPCR or membrane protein family of interest.The fact that only one conformational state is needed could provide a useful way to probe the conformational landscape of proteins that are difficult to solve experimentally or are less well studied in structural biology.Combining the applicability of the methodology with the fast time scales could have wide implications of using the methodology for medium-throughput drug screening or discovery.
Free energy estimation convergence of each AWH simulation run (Figure S1); covering and overlap of different walkers throughout the AWH simulations (Figure S2); quantification of the error of the free energy surfaces (Figure S3); and all contact networks and shortest path, a subset of which are selected and presented in Figure 5 (Figure S4) (PDF) ■

Figure 1 .
Figure1.Coevolution-driven conformational exploration framework.The prerequisites include deep MSA (A) and a reference structure (PDB code used herein: 3P0G).The first step is to identify the false positive residue pairs, i.e., coevolving residues not in contact in the reference structure (B), after which these pairs are clustered according to their coevolution coupling score using the DCA algorithm as implemented in GREMLIN.Belonging to a cluster is color-coded as a circle on the coupling score plot and on the pairs shown in the molecular model.(C) The initial CVs are linear combinations of the pairs within each cluster, where the coefficients are the normalized coevolution scores.These are used in pulling simulations to achieve the required sampling (D.).Machine learning methods are then applied to reduce the space of all pairs into informative pairs that formed and broke contacts during pulling simulation (E).Steps D and E can be iterated until all slow degrees of freedom are included, after which an AWH simulation is run to converge the conformational landscape, as defined by the new CV space.

Figure 2 .
Figure 2. Structural analysis of the explored conformations extracted from the MD simulations.(A) RMSD toward an active (PDB code: 4LDE) and inactive (PDB code: 2RH1) reference.The structures compared were taken from the start and end points of the exploration simulations, as well as extracted from the resulting resting free energy basin from AWH. (B) Structural differences between the active basins were divided into A1 and A2.Note that most of the differences stem from hinges at two positions in the TM6 and TM4 helices.

Figure 3 .
Figure 3. Free energy surfaces of the conformational landscape of β2AR in the given CV space in the absence and presence of different ligands, ranging from inverse agonists to antagonists, partial agonists, and full agonists.The different basins or areas of the free energy landscapes are marked with R (resting), I (intermediate), A1 (Active state 1), and A2 (Active state 2).

Figure
Figure 4. Correlation between experimental (Emax values relative to adrenaline along the x-axis)37−43  and computational activities (calculated fractions of the active population basins along the y-axis).The correlation line was fitted with least-squares, and the R-value was calculated based on this fitted line.Due to differing efficiency values in assays measuring G-protein and β-arrestin signaling of salmeterol and salbutamol, we calculated both the fraction of active population of all active states (green, top right) and only the A 1 state as defined by the basin seen in the adrenaline simulations (purple pentagons, middle).The error bars along x are taken from the references, and the error bars along y are calculated based on the population difference that would arise from the calculated free energy landscape errors.* = Fraction of active population calculated based on the population of the A1 basin, where the A 1 basin is defined by where it is located in the adrenaline-bound simulations.These points were not included in the calculation of the correlation factor R ** = Assumed computational values averaged between Ichikawa et al. (2016)42 and Fleetwood et al. (2020).45These points were not included in the calculation of the correlation factor R.

Figure 5 .
Figure 5. Ligands affect the activation of different allosteric networks.(A) PCA of the per-residue energetical contributions.(B) The β2AR apo system is colored by the per-residue relative betweenness centrality, with three main allosteric paths from the ligand-binding pocket (F19345.52  ) to the G-protein binding site (E268 6.30 , R131 3.50 , P3308.48  ) highlighted as sticks.Analogously plotted are the β2AR carazolol system (C), the β2AR timolol system (D), the β2AR adrenaline system (E), and the β2AR salmeterol system (F).The shortest paths shown in this figure have been made available for all ligand-bound systems in FigureS4.

Table 1 .
AWH Parameters Used for All Simulations

Table 2 .
Classification of All of the Energy Basins for Each Simulation Condition Based on the TM6 Average RMSD (as Defined by Residues 269-298) a a The cutoff value for classification was 3.0 Å; criteria for active states are highlighted in red.