Free energy calculations of the functional selectivity of 5-HT2B G protein-coupled receptor

G Protein-Coupled Receptors (GPCRs) mediate intracellular signaling in response to extracellular ligand binding and are the target of one-third of approved drugs. Ligand binding modulates the GPCR molecular free energy landscape by preferentially stabilizing active or inactive conformations that dictate intracellular protein recruitment and downstream signaling. We perform enhanced sampling molecular dynamics simulations to recover the free energy surfaces of a thermostable mutant of the GPCR serotonin receptor 5-HT2B in the unliganded form and bound to a lysergic acid diethylamide (LSD) agonist and lisuride antagonist. LSD binding imparts a ∼110 kJ/mol driving force for conformational rearrangement into an active state. The lisuride-bound form is structurally similar to the apo form and only ∼24 kJ/mol more stable. This work quantifies ligand-induced conformational specificity and functional selectivity of 5-HT2B and presents a platform for high-throughput virtual screening of ligands and rational engineering of the ligand-bound molecular free energy landscape.


Introduction
G Protein-Coupled Receptors (GPCRs), also called seven-pass-transmembrane domain receptors, are a large family of membrane-bound eukaryotic proteins. These proteins undergo conformational change in response to the binding of extracellular drugs or ligands that results in receptor activation, intracellular G-protein recruitment, and downstream signaling [1][2][3]. GPCRs share a common structure comprising seven transmembrane helices, three extracellular and three intracellular loops, and an extracellular N-terminus and intracellular C-terminus [1]. There is broad pharmaceutical interest in GPCRs due to their importance as targets for therapeutic drugs in the treatment of conditions such as hypertension, schizophrenia, depression, obesity, and Alzheimer's disease [4]. GPCRs constitute the target of approximately 34% of Food and Drug Administration (FDA)-approved drugs with an annual market of over US $150 billion [4].
Crystallographic structures of a number of GPCRs, including rhodopsin [5,6], β 2 -adrenergic receptor (β 2 AR) [7][8][9], and serotonin receptor 2B (5-hydroxytryptamine receptor 2B, 5-HT 2B ) [10][11][12], have provided new understanding of ligand and drug binding and the change induces downstream intracellular signaling events, typically by promoting the binding of a cellular partner protein such as G-protein or arrestin [1,10]. Experimental crystal structures and structure-directed mutagenesis experiments conducted on 5-HT 2B by McCorvy et al. [10] showed that the differential binding poses of lysergic acid diethylamide (LSD)-a prototypical agonist-and lisuride-a prototypical antagonist that differs from LSD only in its stereochemistry and two extra atoms in a NH group-were responsible for the different pharmacological activity [10]. Specifically, contact of LSD with the Leu362 residue in TM7 was shown to be seemingly crucial for its agonism and resulting intracellular G-protein or arrestin recruitment [10]. LSD, commonly known as the psychadelic drug "acid" is thought to mediate hallucinogenic "acid trips" through binding to the closely related 5-HT 2A serotonin receptor [38]. Therapeutic interest in LSD have emerged from its potential in treating alcoholism, depression, and anxiety in the terminally ill [38]. 5-HT 2B and lisuride are also pharmacologically interesting for their role in cardiovascular health: transgenic mice lacking 5-HT 2B receptors have heart issues, and this effect is mimicked in normal mice dosed with lisuride [39].
At thermodynamic equilibrium the unliganded 5-HT 2B can be conceived of existing in a dynamic equilibrium wherein the accessible conformational states are populated according to the Boltzmann distribution. Binding of a ligand into the OBP or EBP can modify the relative stabilities of the conformational states of the receptor within this ensemble resulting in a rebalancing of the equilibrium population among the various microstates. Within this thermodynamic picture, agonist binding can be thought of as preferentially stabilizing the active state and inducing a thermodynamically favorable driving force for the conformational rearrangement along the activation pathway. Mechanistically, agonist binding induces changes in trigger motifs in the vicinity of the OBP or EBP that facilitate large-scale conformational rearrangements of the transmembrane helices that prime the protein for intracellular G-protein or arrestin binding [10]. Antagonist binding, in contrast, preferentially stabilizes the inactive state and provides no thermodynamic driving force for adoption of the active state.
In this work, we consider an thermostable mutant of 5-HT 2B that we term 5-HT 2B -TM, which has truncated N and C-termini and contains a thermostabilizing replacement of native residues Tyr249-Val313 with Ala1-Leu106 from the thermostabilized apocytochrome b562 RIL [11,12]. The thermostabilizing mutations restrain the conformational flexibility of the receptor and make it well suited to crystallographic structure determination as reported in Refs. [11,12]. It is the primary goal of the present work to conduct molecular modeling of the engineered 5-HT 2B -TM mutant used in crystallography studies, draw comparisons against the experimental structural results, and conduct free energy calculations to shed light on the mechanism of ligand-induced conformational specificity and functional specificity. Our results are primarily of interest in providing deeper thermodynamic and mechanistic understanding of the 5-HT 2B -TM. The transferability of our results to the wild-type receptor may be assessed by repeating our calculations for the 5-HT 2B system. Specifically, we map out changes to the 5-HT 2B -TM free energy surface in the apo (i.e., unliganded), LSD-bound, and lisuridebound states. It is our hypothesis that the LSD agonist will stabilize an active-like state of the receptor whereas lisuride will stabilize an inactive-like state, and that these differences can provide the mechanistic basis for the observed functional selectivity [19]. The free energy change DG ligand a of ligand binding to the 5-HT 2B -TM receptor is associated with the process, ligand þ 5-HT 2B -TMðapoÞ Ð ligand À 5-HT 2B -TMðrelaxedÞ complex; wherein a ligand binds to a receptor in the apo conformation and then the ligand-receptor complex relaxes to the stable ligand-bound state (Fig 2a). This thermodynamic process may be (hypothetically) broken into two steps by virtue of free energy being a state function, ligand þ 5-HT 2B -TMðapoÞ Ð ligand À 5-HT 2B -TMðapoÞ complex ð2Þ corresponding first to binding of the ligand into the binding pocket of a 5-HT 2B -TM receptor in an apo-like conformation and then a subsequent relaxation of the ligand-receptor complex to the stable conformational state by structural rearrangement of the 5-HT 2B -TM receptor from an apo-like conformation to its thermodynamically stable ligand-bound state. It is the second step that is of principal interest in this work. We term the associated free energy change DG ligand APO!ligand , where the subscript indicates the conformational change from the apo conformation to the is the free energy change associated with the process ligand + receptor (apo) Ð ligand-receptor (relaxed) complex, wherein a LSD or lisuride ligand binds to a 5-HT 2B -TM receptor in the apo conformation and then the ligand-receptor complex relaxes to the stable ligand-bound state. The ligand binding process can be split into two successive steps ligand + receptor (apo) Ð ligand-receptor (apo) complex Ð ligandreceptor (relaxed) complex. (b) The free energy of structural relaxation DG ligand APO!ligand corresponds to the free energy change of the second step, wherein the ligand-receptor complex with the 5-HT 2B -TM receptor in the apo-like conformation relaxes to its thermodynamically stable ligand-bound state. The terminal state and thermodynamic driving force for this relaxation differ for LSD and lisuride ligands.
https://doi.org/10.1371/journal.pone.0243313.g002 stable ligand-bound conformation and the superscript reminds us that this process proceeds with the ligand bound (Fig 2b). The free energy for receptor-ligand binding DG ligand a is straightforwardly computed from experimental dissociation constants. The free energy for structural relaxation of the ligand-bound receptor DG ligand APO!ligand is much more challenging to measure experimentally and is computationally estimated in the molecular simulations performed in this work. The structure of the ligand-receptor complex at the end of the thermodynamic relaxation process is of interest in determining whether the stable ligand-bound state is active-like or inactive-like, thereby providing evidence for a mechanistic basis for 5-HT 2B -TM functional selectivity. The free energy change associated with the structural relaxation provides a measure of the downhill driving force and thermodynamic stability of the relaxed ligand-bound conformation of the receptor relative to the apo state.
It is the primary objective of this work to quantify the perturbations to the conformational free energy surfaces and thermodynamic driving forces for structural relaxation of the ligand-5-HT 2B -TM complex under LSD and lisuride binding. Since free energy is a state function, the relative stabilities of various configurations are a pure function of state (i.e., independent of history or pathway) and can be measured by estimating the free energy landscape of the receptor mapping out the relative free energy of different receptor configurations. We determine the conformational free energy landscapes of 5-HT 2B -TM in (i) the unliganded (apo) form, (ii) LSD-bound form, and (iii) lisuride-bound form using allatom molecular dynamics simulations. We employ well-tempered metadynamics [32][33][34][35] to provide good exploration of the conformational free energy landscape for each of the three systems by driving two collective variables quantifying conformational changes in the OBP and resolve the metastable and stable states of each system. The global free energy minimum conformation of the receptor in the LSD-5-HT 2B -TM complex is structurally similar to the activated state, whereas the global free energy minimum conformation of the receptor in the lisuride-5-HT 2B -TM complex is structurally similar to the inactive apo state. We employ free energy calculations to drive the ligand-bound complexes between their stable and apo-like conformations to estimate thermodynamic driving forces for structural relaxation of the LSD-5-HT 2B -TM complex of DG LSD APO!LSD � À 110 kJ=mol, and of the lisuride-5-HT 2B -TM complex of DG lisuride APO!lisuride � À 24 kJ=mol. Our results provide a deeper understanding of the ligand-induced conformational specificity and functional specificity of the 5-HT 2B -TM seratonin receptor and establish a framework for highthroughput screening or molecular design of pharmaceuticals targeting this physiologically important receptor.

System preparation and relaxation
The unliganded membrane-bound 5-HT 2B -TM protein (which we henceforward refer to as the APO system for short) was prepared from the LSD-bound crystal structure (PDB: 5TVN) reported in Ref. [11] by removing the ligand from the structure. The 5TVN sequence is an engineered 5-HT 2B -TM construct that differs from the wild type in lacking a number of N and C-terminal residues and containing thermostabilizing mutations [11]. The protonation state of all residues in the receptor were specified as that corresponding to pH 7 such that the protein carried a (+1) net charge. Protonation states were specified using a classical model of pKa calculations implemented in the ACD Percepta software [40]. We terminated GLU41 with a NH 3 + group and the ARG400 with COO − group using the Gromacs software pdb2gmx and specified residue types HISE for HIS55,242, and HISH for HIS145 and HISD for HIS1064.
Both LSD and lisuride were predicted to have pKa of 7.4 for the tertiary amine nitrogen atom that was prepared in its protonated form. The LSD-bound (LSD system) and lisuride-bound (LIS system) systems were prepared in an analogous manner to the APO system from the LSD-bound (PDB: 5TVN) and lisuridebound (PDB: 6DRX) crystal structures reported in Refs. [10,11] wacker2017crystal, mccorvy 2018structural. The lisuride from 6DRX was transferred into the EBP of the 5TVN protein structure. The LSD ligand was prepared with its tertiary amine nitrogen protonated [11], and parameterized within the semi-empirical QM0 force field using the Automated Topology Builder (ATB) server [46]. The lisuride ligand differs from LSD only in its sterochemistry and an extra NH group, and was prepared in an identical manner. Since both LSD and lisuride carry a (+1) net charge, it was not necessary to add any counterions to these two systems to maintain charge neutrality.
Each of the three systems-APO, LSD, LIS-were then subjected to 1400 steps of steepest descent energy minimization to remove high-energy atomic overlaps, followed by 200 ns of molecular dynamics simulation at at T = 323 K and P = 1 bar to accelerate structural relaxation. Systems were then cooled to T = 310.15 K (body temperature) over the course 100 ns under a linear cooling schedule. Finally, 1.5 μs production runs were conducted at T = 310.15 K and P = 1 bar. Equations of motion were integrated using the leapfrog algorithm with a time step of 2 fs [47]. Temperature was maintained using a Nosé-Hoover thermostat [48] and pressure maintained using the Parrinello-Rahman barostat [49]. Lennard-Jones interactions were switched smoothly to zero at a 1.2 nm cutoff and dispersion interactions between unlike atoms specified by Lorentz-Berthelot combining rules [50]. Electrostatics treated using particle-mesh Ewald with a real space cutoff of 1.2 nm and a reciprocal space grid spacing of 0.16 nm [51]. Bond lengths were fixed to their equilibrium values using the LINCS algorithm [52]. All calculations were performed using the Gromacs 2018 simulation suite [53]. Our simulation input files are provided as S1 File and have been uploaded to GPCRmd [54] and are locatable at https://submission.gpcrmd.org/view/213/ (APO), https://submission.gpcrmd.org/view/214/ (LSD), and https://submission.gpcrmd.org/view/217/ (LIS).

Well-tempered metadynamics
Following Provasi et al. [19] provasi2011ligand, the conformational free energy surfaces of the APO, LSD, and LIS systems were estimated using well-tempered metadynamics (WTMD) [32][33][34][35] to enhance sampling of configurational space and recover an estimate of the molecular free energy landscape. Accelerated sampling is conducted in a small number k of user-defined collective variables s ¼ fs i ðRÞg 2 R k that are themselves a function of the atomic coordinates R 2 R 3N , where N is the number of atoms in the system. In brief, this approach lays down a history-dependent bias within the collective variables that drives the simulation to explore previously unvisited regions of configurational space. The metadynamics time-dependent bias potential takes the form of a time-dependent sum of Gaussians, where t 0 is an integer multiple of the deposition time τ that dictates the frequency with which Gaussians are deposited, ω t 0 controls the height of the Gaussians, and σ i controls their width.
Since we employ well-tempered metadynamics, ω t 0 is time dependent and takes the form, where ΔT is a parameter with units of temperature, k B is Boltzmann's constant, and ω is the maximum height of the deposited Gaussians. ΔT controls how aggressively the metadynamics bias is applied: in the limit ΔT ! 1 we recover standard metadynamics, whereas for ΔT = 0 we recover standard molecular dynamics. For this reason it can be convenient to express ΔT in terms of the so-called bias factor, where T is the temperature of the simulation. As t ! 1, the bias potential converges and the free energy landscape in the collective variables can be estimated-up to an arbitrary additive constant-as, The WTMD calculations were conducted using the PLUMED plugin [55] to Gromacs 2018. We adopted as our k = 2 collective variables the root mean squared deviation (RMSD) of the OBP (cf. Fig 1b) rotationally and translationally aligned to that of the LSD and LIS systems at the end of the 1.5 μs unbiased simulations. The OBP was defined as the contact residues between LSD ligand and the 5-HT 2B -TM receptor within the LSD-bound 5-HT 2B -TM crystal structure (PDB: 5TVN) [11] computed using MSDsite [56]. So defined, the OBP comprises residues ALA130, LEU132, PHE133, ASP135, VAL136, SER139, THR140, ILE186, VAL190, ILE205, CYS207, VAL208, LEU209, PHE217, MET218, GLY221, ALA224, ALA225, TRP337, PHE340, PHE341, ASN344, ILE345, VAL348, LEU362, PHE365, VAL366, GLY369, and TYR370. This pair of collective variables, (RMSD LIS , RMSD LSD ), measures the deviation of the OBP from the lisuride-bound and LSD-bound conformations, and is therefore expected to adequately span the space of configurational motions within the OBP observed in the three systems, resolve metastable states associated with the active and inactive forms of the receptor, and provide interpretable 2D projections of the the free energy surfaces. Clearly this choice of collective variables directly considers only those residues within the OBP, but since ligandinduced perturbations in the OBP are known to induce global conformational changes in the 5-HT 2B -TM we anticipated that they should also serve as good descriptors of the global protein state. In other words, we make the assumption that under good sampling induced by WTMD, an active (inactive) conformation of the OBP should correspond to an active (inactive) conformation of the entire protein. We demonstrate below that this is indeed the case by verifying that the well-tempered metadynamics calculations achieve convergence, that the free energy landscapes smoothly span the conformational space corresponding to the (meta)stable APO, LSD, and LIS structures, and that the structural changes in key motifs that are indicative of receptor activation are consistent with experimental reports. We observe in passing that more sophisticated approaches to learn good collective variables from the data, as opposed to prespecifying these in advance, could have been employed [57,58].
After some preliminary exploration, we found the following choice of parameters to perform well for our WTMD calculations, σ 1 = σ 2 = 0.01 nm, ω = 1.0 kJ/mol, and γ = 5.0. We obtained converged free energy landscapes spanning *120 kJ/mol for each of the three systems over the course of 400 ns production runs. We test for convergence of the simulations by monitoring that there are no longer Gaussians larger than 5% of ω being deposited over the course of the terminal 100 ns of the simulation. We also recovered the free energy landscapes within five 20 ns blocks spanning the terminal 100 ns portion of the simulation to verify that the landscapes are converged within this block analysis to better than 0.1 kJ/mol. We define the active structure predicted by the converged WTMD calculations to be the structure occupying the global free minimum of the free energy landscape. Where possible we present a comparison of the computational prediction with available experimental data on the active state. Our PLUMED metadynamics protocol has been uploaded to PLUMED-NEST [59] and is available at https://www.plumed-nest.org/eggs/20/026.

Umbrella sampling and WHAM
Having recovered estimates of the free energy surfaces parameterized by two collective variables RMSD LIS and RMSD LSD using WTMD, we estimated the thermodynamic driving forces for conformational relaxation of the ligand-receptor complex from an initial apo-like state of the receptor induced by LSD and lisuride binding (Fig 2b). The LIS free energy surface was sufficiently well sampled for us to simply read off the free energy difference DG LIS APO!LIS between the lisuride-bound global free energy minimum and the coordinates of the global free energy minimum of the APO system (i.e., the lisuride-receptor complex with the receptor in an apolike configuration). The LSD free energy surface, however, was sufficiently perturbed from the APO landscape that the apo-like configuration was not sufficiently sampled to allow us to read off DG LSD APO!LSD . Accordingly, we employed umbrella sampling (US) and the weighted histogram analysis method (WHAM) to construct a pathway between these states and estimate this quantity [60,61].
We first identified the (RMSD LIS , RMSD LSD ) coordinates of the global free energy minimum of the LSD system corresponding to the most stable conformation of the LSD-bound 5-HT 2B -TM receptor. This structure served as the final state for our US calculations. We then identified the (RMSD LIS , RMSD LSD ) coordinates of the global free energy minimum of the APO system corresponding to the most stable configuration of the unliganded 5-HT 2B -TM. The system configuration residing at these coordinates in the LSD free energy surface corresponds to the apo-like-configuration into which a LSD ligand has been inserted. This serves as the initial state for our US calculations. The free energy difference between the initial and final configurations corresponds to the free energy change DG LSD APO!LSD associated with the configurational relaxation of the LSD-bound system from an apo-like starting configuration (Fig 2b).
We charted a path between the (RMSD LIS , RMSD LSD ) coordinates of two end states by constructing a series of umbrella windows. As a state function, the free energy difference is independent of the path, so in regions where the LSD free energy landscape determined by WTMD is available, we construct the path to follow its low-free-energy contours for numerical stability. We employ two-dimensional harmonic restraining potentials in the k = 2 collective variables (s 1 , s 2 ) = (RMSD LIS , RMSD LSD ), where s i (R) is the value of collective variable s i corresponding to a system configuration R, κ is the strength of the restraining potential spring constant, and ξ = {ξ 1 , ξ 2 } is the location of the restraint. We used spring constants of κ = 100,000 kJ/mol.nm 2 and employed 26 umbrella windows. Each umbrella simulation was conducted for 200 ns, and the unbiased potential of mean force (PMF) along the pathway of umbrella windows was reconstructed from the biased simulation trajectories using WHAM [61] implemented in PLUMED [55]. Uncertainties in the calculated free energies are estimated by 100 rounds of bootstrap resampling. The probability distribution functions under each umbrella window rapidly converged and became time invariant on the time scale of tens of ns. An image illustrating adequate overlap between the distribution functions under the umbrella windows is presented in S1 Fig. The quantities DG LIS APO!LIS and DG LSD APO!LSD correspond to the thermodynamic driving forces for structural relaxation of the ligand-receptor complex from the apo-like state of the receptor induced by ligand binding and arise from perturbation of the underlying conformational free energy surface mediated by entry of the ligand into the binding pocket (Fig 2b). We emphasize that the initial and final states of the systems are both ligand-bound, so these free energy differences do not contain any contributions from ligand binding or release.

Results and discussion
We conducted enhanced sampling molecular dynamics simulations using well-tempered metadynamics to estimate the conformational free energy surfaces for the 5-HT 2B -TM seratonin receptor in the unliganded (APO system), LSD-bound (LSD system), and lisuride-bound (LIS system) states. Comparison of these free energy landscapes, each of which is parameterized by the RMSD of the OBP residues relative to the lisuride-bound and LSD-bound systems (RMSD LIS , RMSD LSD ), provide a quantitative measure of the perturbation to the conformational free energy landscape of the receptor induced by ligand binding, the driving force for conformational rearrangement of the ligand-receptor complex, and new molecular-level understanding for ligand-induced conformational specificity and functional selectivity.

Conformational free energy landscapes
We present in Fig 3 the conformational free energy landscapes for the APO, LSD, and LIS systems estimated using well-tempered metadynamics (WTMD) calculations.
APO system. The free energy surface for the APO system presented in Fig 3a exhibits a single-well free energy landscape with a deep global minimum at (RMSD LIS = 0.40 nm, RMSD LSD = 0.46 nm) corresponding to the unliganded and inactive form of 5-HT 2B -TM. We recall that RMSD LSD = 0.0 nm indicates that the OBP precisely matches the conformation of that in the LSD-bound system after 1.5 μs of unbiased simulation, and RMSD LIS = 0.0 nm matches the OBP state in the lisuride-bound system after an equivalent relaxation time. As anticipated, the most stable configuration of the unliganded apo protein is not coincident with either the LSD-bound or lisuride-bound structures. The structure is presented in Fig 3e (teal).
LSD system. The LSD system (Fig 3b) exhibits a qualitatively different free energy landscape from the APO system. Relative to APO, the LSD free energy surface is shifted towards lower values of RMSD LSD , spans a larger range of RMSD LIS , and exhibits two free energy minima in the surface. Our well-tempered metadynamics simulations recover the free energy landscape up to *114 kJ/mol above the global free energy minimum, and we observe essentially no overlap between the APO and LSD landscapes within this free energy range in the (RMSD LIS , RMSD LSD ) projection (Fig 3d). The implication of the observed absence of overlap is that low-free energy conformations of the unliganded receptor are high-free energy configurations of the receptor in the LSD-5-HT 2B -TM complex, and vice versa. This provides a clear indication that binding of the LSD agonist induces large perturbations in the conformational free energy landscape of the 5-HT 2B -TM receptor and mediates large changes in the distribution of the conformational ensemble of the receptor.
The receptor conformation residing at the global free energy minimum at (RMSD LIS = 0.45 nm, RMSD LSD = 0.25 nm) is structurally close to the LSD-bound crystal structure reported in Ref. [11] (PDB: 5TVN), with a RMSD Ca 5TVN ¼ 0:56 nm computed over all C α atoms. We identify this structure as our computational prediction of the active state and refer to it as the active state or stable state for brevity. The structure is presented in Fig 3e (red), and can be seen to differ substantially from the most stable conformation of the APO system (teal) most obviously by a large conformational change in the intracellular helices that move up towards the inner leaflet of the membrane. LSD is a prototypical 5-HT 2B agonist [10,12,62] and we identify the large conformational change that occurs upon LSD binding as a transition towards an activated state. This is also consistent with prior molecular simulation work that showed structural movement of the intracellular loops and helices associated with binding of the active protein with G protein or arrestin [1,13,63,64]. The conformational change is quantified by comparing the stable LSD structure to the stable APO structure over the TM5-7 regions of the receptor to compute RMSD Wacker et al. have previously identified a number of motifs involved in 5-HT 2B -TM receptor activation [11,12]. The PIF motif that forms an interface between TM3, TM5, and TM6 was observed to be partially activated in a LSD-5-HT 2B -TM crystal structure, comprising a shift inward of the TM5 Pro residue, dihedral flip of TM3 Ile, and rotation of the TM6 Phe [11]. Consistent with these findings, relative to the APO structure (teal) the stable LSD structure (red) exhibits an inward rotation of the Pro229, rotation of the Ile143, and a translational shift of the Phe333 and TM6 helix (Fig 5a). The NPxxY motif in TM7 and D(E)/RY motif in TM3 are known as "microswitches" that are implicated in conformational change and priming of the receptor for activation and G-protein binding [8,12]. The NPxxY motif was observed to be activated in the LSD-5-HT 2B -TM crystal structure, with marked rotation of the TM7 Tyr into the TM bundle [11]. We find consistency with these experimental findings, with the TM7 and Tyr380 side chain in the stable LSD structure (red) displaced and rotated relative to the APO (teal) (Fig 5b). In the case of the D(E)/RY motif, we observe relatively little change in the TM3 location or relative position of the side chains between the stable LSD (red) and APO (teal) structures (Fig 5c). This is consistent with experimental reports that this motif shows relatively little activation in the LSD-5-HT 2B -TM crystal structure [11] and that the D152-R153 salt bridge is preserved in its inactive conformation in the ergotamine-5-HT 2B -TM crystal structure, where ergotamine is the chemical precursor of LSD [12]. We performed these comparisons for the active structures residing at the bottom of the global free energy minima in the APO and LSD conformational free energy landscapes. In S1-S3 Videos we present portions of the WTMD trajectory occupying the global free energy minima to show good structural homogeneity in the PIF, NPxxY, and D(E)/RY motifs and demonstrating that the global free energy minimum structure is representative of the structural ensemble in the vicinity of the global free energy well.
The local free energy minimum at (RMSD LIS = 0.65 nm, RMSD LSD = 0.43 nm) lies ΔG � 25 kJ/mol higher in free energy than the global minimum and is substantially closer to the stable APO structure in the TM5-7 regions with a RMSD TM5À 7 APO ¼ 0:38 nm. Accordingly, we identify this structure as a more inactive-like LSD-bound structure (Fig 3e (pink)) that is metastable relative to the active-like global minimum (Fig 3e (red)).
Finally, we note that the global free energy minimum of the LSD system lies at RMSD LSD = 0.25 nm, and that conformations at RMSD LSD = 0.0 nm lie more than 100 kJ/mol higher in free energy than the global minimum. Why should configurations of the LSD system that lie  closer to RMSD LSD = 0.0 nm be less stable? We recall our definition of our reference state for computation of RMSD LSD as the terminal configuration of the 1.5 μs unbiased simulation of the LSD-bound system. Although this terminal state appeared to be a relaxed stable configuration, the free energy landscapes recovered by WTMD suggest that it was in fact a kinetically trapped configuration that was unable to relax into the global free energy minimum on the time scale of our unbiased calculations. A similar scenario was observed for the LIS system (vide infra). These observations are not surprising given the slow time scales associated with GPCR dynamics: passage times of several microseconds for the active to inactive transitions of the β 2 AR GPCR have been reported, and passage times for the inactive to active transition occur on still longer time scales and have not been observed in unbiased molecular dynamics simulations [14]. Despite this deficiency in our reference states, (RMSD LIS , RMSD LSD ) still serve as a good discriminatory collective variable in which to conduct the enhanced sampling calculations and parameterize our free energy surfaces. These observations do caution, however, against over-interpreting the structural results of unbiased calculations due to the very slow relaxation times.
LIS system. The LIS system (Fig 3c) exhibits a qualitatively similar free energy landscape to the APO system. There is substantial overlap in the APO and LIS landscapes indicating in intersection of the low-free energy structural ensembles (Fig 3d), indicating the binding of the lisuride antagonist causes a relatively mild perturbation to the 5-HT 2B -TM receptor free energy landscape within the lisuride-5-HT 2B -TM complex compared to LSD. The global free energy minimum of the LIS system resides at at (RMSD LIS = 0.31 nm, RMSD LSD = 0.41 nm), which is proximate to the APO minimum at (RMSD LIS = 0.40 nm, RMSD LSD = 0.46 nm). Structurally, these two conformations are quite close, with the LIS native structure differing from the APO native structure by only RMSD Ca APO ¼ 0:66 nm over the C α atoms. Comparison over the TM5-7 domains returns a relatively large RMSD TM5À 7 APO ¼ 0:38 nm, the magnitude of which is *85% of that under LSD binding. Comparison of the APO (teal), LSD (red), and LIS (black) stable conformations in Fig 3e shows that the TM5-7 helices in the LIS stable structure continue to point down in a similar pose to the APO structure, such that lisuride binding does not induce TM5-7 to move up and towards the membrane as observed in response to LSDmediated receptor activation. Lisuride is a prototypical antagonist for 5-HT 2B -TM and the structure of the LIS stable state supports that it stabilizes an inactive-like state of the receptor. As was observed for LSD, and consistent with McCorvy et al. [10], the OBP undergoes only a small change in response to lisuride binding with RMSD OBP APO ¼ 0:28 nm. The OBP region of the APO and LIS structures are overlaid in Fig 4b. McCorvy et al. determined a lisuride-bound crystal structure and in comparisons to the LSD-bound structure identified similar binding poses of LSD and lisuride within the OBP but different poses in the EBP [10]. This led the authors to propose that the pharmacologically distinct activities of LSD (agonist) and lisuride (antagonist) are mediated through the EBP, and in particular the Leu362 residue on TM7 [10]. Our computational predictions are consistent with these experimental findings, with the stable LSD structure (red) exhibiting a large displacement of the Leu362 residue relative to the APO (teal) compared to a much smaller relative displacement in the stable LIS structure (black) (Fig 5d). We observe similar, and even accentuated, structural changes in the PIF and NPxxY motifs in the stable LIS structures (black) systems. c) The D(E)/RY motif on TM3 shows relatively little change in terms of the TM3 location and Asp152, Arg153, and Tyr154 side chain positions in the LSD (red) and LIS (black) systems relative to the APO (teal). d) The Leu362 residue in the EBP on TM7 exhibits a large relative displacement in the LSD (red) system compared to the LIS (black).
https://doi.org/10.1371/journal.pone.0243313.g005 (black) relative to the LSD system (red) that are consistent with (partial) activation (Fig 5a and  5b) [11,12]. As was the case for the LSD system (red), there is relatively little structural change in the D(E)/RY motif in the stable LIS structure (black) (Fig 5c). A comparison of these motifs between the LIS and LSD systems revealed RMSD PIF = 0.23 nm, RMSD D(E)/RY = 0.16 nm, and RMSD NPxxY = 0.20 nm, indicating small, but significant, changes in the structure of these regions. Experimental characterizations of the structural changes in the PIF, NPxxY, and D (E)/RY motifs in the lisuride-bound complex against which to compare these computational predictions remain to be conducted [10]. S1-S4 Videos show portions of the WTMD trajectory within the global free energy minima and illustrate good structural homogeneity in the PIF, NPxxY, D(E)/RY, and Leu362 motifs.
The overlap between the APO and LIS free energy landscapes allows us to read-off the free energy difference between the coordinates of the APO (RMSD LIS = 0.40 nm, RMSD LSD = 0.46 nm) and LIS (RMSD LIS = 0.31 nm, RMSD LSD = 0.41 nm) minima in the APO and LIS systems. In the APO system, this corresponds to the free energy cost associated with changing the conformation of the unliganded receptor from its global free energy minimum to a structure resembling the native structure of the lisuride-bound system in the absence of the ligand. We term this quantity DG APO APO!LIS . In the LIS system, this corresponds to the free energy cost associated with changing the conformation of the lisuride-bound system from its global free energy minimum to a structure resembling the native structure of the unliganded system with the lisuride ligand bound. We term this quantity DG LIS LIS!APO . Interrogation of our APO free energy landscape reveals DG APO APO!LIS ¼ ð38:4 � 6:0Þ kJ=mol, indicating that the lisuride-bound-like conformation is a somewhat low-lying member of the unliganded conformational ensemble. Uncertainties are estimated by 100 bootstrap resamples. Similarly, interrogation of our LIS free energy landscape reveals DG LIS LIS!APO ¼ ð23:6 � 5:0Þ kJ=mol, indicating that unligandedlike conformation is a relatively low-lying member of the lisuride-bound conformational ensemble.
The quantity DG LIS APO!LIS ¼ ðÀ DG LIS LIS!APO Þ ¼ À ð23:6 � 5:0Þ kJ=mol can be interpreted as the favorable thermodynamic driving force for conformational relaxation of the lisuride-5-HT 2B -TM complex by structural rearrangement of the lisuride-bound 5-HT 2B -TM receptor from an apo-like conformation to its thermodynamically stable lisuride-bound state. Physically, this process can be conceived of as the reversible work associated with the conformational relaxation of the 5-HT 2B -TM receptor after inserting the lisuride ligand into the binding pocket of the apo conformation (Fig 2b). The small value of this quantity is consistent with the relatively minor nature of the perturbation to the 5-HT 2B -TM free energy landscape induced by lisuride binding.

Ligand-induced conformational rearrangement free energies
In the previous section, we reported DG LIS APO!LIS by reading this off the LIS free energy landscape. Due to the disjoint nature of the APO and LSD free energy landscapes from WTMD (Fig 3d) we are unable to read off an analogous quantity DG LSD APO!LSD for LSD. Instead, we employed umbrella sampling to construct a path in the LSD system between the coordinates of the global minimum in the LSD system (RMSD LIS = 0.45 nm, RMSD LSD = 0.25 nm) and the global minimum of the APO system (RMSD LIS = 0.40 nm, RMSD LSD = 0.46 nm). We illustrate in Fig 6 the free energy pathway constructed between these states using umbrella sampling, which provides an estimate of DG LSD APO!LSD ¼ À ð110:3 � 0:4Þ kJ=mol. This is the favorable thermodynamic driving force for conformational relaxation of the LSD-5-HT 2B -TM complex by structural rearrangement of the LSD-bound 5-HT 2B -TM receptor from an apo-like conformation to its thermodynamically stable LSD-bound state. Physically, this quantity can be interpreted as the reversible work associated with the conformational relaxation of the 5-HT 2B -TM receptor after inserting the LSD ligand into the binding pocket of the apo conformation (Fig 2b).
The value of DG LSD APO!LSD ¼ À ð110:3 � 0:4Þ kJ=mol is approximately five times larger than DG LIS APO!LIS ¼ À ð23:6 � 5:0Þ kJ=mol. This reflects the much larger perturbation of the 5-HT 2B -TM conformational free energy landscape induced by the LSD agonist compared to the lisuride antagonist, and a much greater thermodynamic gradient for conformational rearrangement ultimately leading to receptor activation. The driving force for conformational relaxation of the lisuride-5-HT 2B -TM complex is the same order of magnitude as the values reported by Provasi et al. for the conformational relaxation of β 2 AR bound to a full agonist, weak partial agonist, two inverse agonists, and one neutral agonist [19]. Our calculated value for LSD binding is approximately ten times larger, indicating a far greater driving force for the LSDinduced activation of 5-HT 2B -TM than any of the β 2 AR agonists explored in that work.
It is also instructive to compare our calculated free energies for structural relaxation of the ligand-receptor complex (DG ligand APO!ligand , Fig 2b) with the experimentally determined ligand binding free energies (DG ligand a , Fig 2a). Millan et al. report an experimental dissociation constant for lisuride binding at T = 295 K of pK d = -log 10 K d = 9.87 [65], from which we compute a strongly favorable standard free energy of association-at a standard reference concentration McCorvy et al. [10] report an experimental dissociation constant for LSD binding at T = 310 K of pK d = 9.33, from which we compute a very similar standard free energy of association of DG LSD a ¼ RT ln K d ¼ À 55:4 kJ=mol. The calculated standard binding free energies are about two times larger than the conformational rearrangement free energy upon lisuride binding DG LIS APO!LIS ¼ À ð23:6 � 5:0Þ kJ=mol, and about two times smaller than the conformational rearrangement free energy upon LSD binding DG LSD APO!LSD ¼ À ð110:3 � 0:4Þ kJ=mol.

Conclusions
In this work, we have reported the use of enhanced sampling molecular dynamics simulations to determine the conformational free energy landscapes of the 5-HT 2B -TM thermostable mutant of the 5-HT 2B seratonin receptor in the apo state, LSD-bound state, and lisuridebound state. LSD is a prototypical agonist for 5-HT 2B , and its binding induces a large perturbation to the conformational free energy landscape of the receptor within the ligand-receptor complex to the degree that the structural ensembles of the apo and LSD-bound states become effectively disjoint. LSD binding shifts the global minimum of the ligand-receptor complex free energy landscape to the active state and induces a thermodynamic driving force for structural activation of ΔG � -110 kJ/mol. We also observe the presence of a metastable inactivelike LSD-bound structure that lies ΔG � 25 kJ/mol higher in free energy. Lisuride, on the other hand, is a prototypical antagonist and its binding induces a relatively smaller perturbation of the conformational free energy landscape. The structural ensembles for the apo and lisuride-bound states show a high degree of similarity, and the global minimum of the lisuridebound free energy landscape shows close structural similarity with the inactive apo form and lies only ΔG � 24 kJ/mol lower in free energy. Our work also reveals that the structural conformations adopted by long 1.5 μs unbiased molecular dynamics simulations do not correspond to the most stable structures identified in our enhanced sampling calculations, revealing the value of accelerated sampling and cautioning against over interpreting unbiased molecular simulations due to the long relaxation time scales associated with these systems. Our results quantify the driving forces for activation of 5-HT 2B -TM by LSD binding and demonstrate the absence of this driving force under lisuride binding, thereby shedding light on the molecular-level structural and thermodynamic basis for ligand-induced conformational specificity and functional selectivity. Future work may consider the binding of other drugs and ligands and the introduction of path sampling techniques such as forward flux sampling to quantify the kinetics of the conformational rearrangement events [66,67]. This work also establishes a framework for high-throughput virtual screening of putative 5-HT 2B ligands, and a principled platform for drug design through the rational engineering of ligand-bound free energy landscapes [19,68].