Determinants of Orexin Receptor Binding and Activation—A Molecular Dynamics Study

We assess the stability of two previously suggested binding modes for the neuropeptide orexin-A in the OX2 receptor through extensive molecular dynamics simulations. As the activation determinants of the receptor remain unknown, we simulated an unliganded receptor and two small-molecular ligands, the antagonist suvorexant and the agonist Nag26 for comparison. Each system was simulated in pure POPC membrane as well as in the 25% cholesterol–POPC membrane. In total, we carried out 36 μs of simulations. Through this set of simulations, we report a stable binding mode for the C-terminus of orexin-A. In addition, we suggest interactions that would promote orexin receptor activation, as well as others that would stabilize the inactive state.


■ INTRODUCTION
Orexin System. Peptide Binding. The orexin receptors are key regulators in several neurological processes. The main physiological role of the orexinergic system is the regulation of sleep, and it also plays a role in energy homeostasis and the reward system. 1 Orexin receptors are prominent drug targets for novel hypnotics, 2,3 and hopefully, in the future, narcolepsy could be alleviated through orexin receptor activation. 4,5 Orexin receptors are also putative targets for cancer treatment, as orexinergic signaling directs certain cancer cell lines to apoptosis. 6 However, the treatment of both narcolepsy and cancer through orexin receptors would require agonistic ligands, whereas the pharmaceutical industry has concentrated on antagonists. A few potent small-molecular agonists have been reported and characterized, 7−9 but the poor understanding of the receptor activation determinants hampers the search for new agonists.
The orexinergic system comprises two peptide ligands, orexin-A and -B, and two G protein-coupled receptors (GPCR) OX 1 and OX 2 . 10 In aqueous solution, the orexin peptides adopt conformations comprising helical segments I and II (and III for orexin-A), with a structural hinge motif between the helices I and II (Figure 1). 11−13 Orexin-B, which is a linear 28-amino-acid amidated peptide, shows a 90°bend between the helices I and II, at N20−A22, i.e., 1.5 turns from the C-terminal. Orexin-A, a 33-amino-acid amidated peptide with two intramolecular disulfide bridges, is often observed with a similar bend at A23−N25, half a turn further from the C-terminus, resulting in bending in the opposite direction in comparison to orexin-B (Figure 1). In addition to the bent conformations, orexin-A appears in some NMR-models in a straight α-helical conformation, except for the disulfide-bridgestabilized hook at the N-terminus. For orexin-B, such conformation has not been reported.
Site-directed mutagenesis (SDM) on both receptor subtypes and both peptides has highlighted several key residues for peptide−receptor interaction. 14−20 The most notable observations is that the C-terminus of both orexin peptides is vital for bioactivity, whereas the N-terminus can be mutated or even deleted without the loss of activity, which has given rise to the assumption that the C-terminus enters the canonical class A GPCR binding cavity. With this information, molecular modelers have provided few suggestions on the binding mode for the orexin peptides, 21,22 with the assumption that the C-terminus binds in the α-helical conformation observed in solution, but validation of these suggestions has been challenging. Recently, we suggested that orexin peptides locked in an α-helical conformation can indeed activate the orexin receptors, 23 but as the method included an introduction of a bulky hydrocarbon linker between two helical turns and severely impaired the peptide's potency, no firm conclusions could be drawn.
GPCR Activation. Simulations of GPCRs. GPCR activation is a process where ligand binding stabilizes a receptor conformation that allows for G protein (or other effector protein) binding, activation and downstream signaling. 24 As GPCR ligands differ greatly, also the key interactions at the binding site must differ. Among the class A GPCRs, the activation cascade that links events at the binding site with the intracellular G protein binding site appears to converge into defined interaction patterns, 25 and large-scale movements such as the outswing of the TM6 (the sixth transmembrane helix), which opens up the G protein binding site. 26 Molecular dynamics (MD) simulations have shed light on different parts of the receptor function, such as ligand binding, receptor transitions both between and within active and inactive states, and allosteric modulation. 27 The opening of the G protein binding site has not been observed in classical MD simulations reaching up to tens of microseconds, but the reverse event has been reported. 28 In contrast, the binding site and the so-called core triad can adopt activelike conformations within a microsecond-scale simulation. 28 In the present study, we set out to simulate the OX 2 receptor with two previously reported alternative orexin-A binding modes, 21 a small molecular agonist Nag26 7 and a small molecular antagonist suvorexant. 2 In addition, we simulated an unliganded receptor. Our main point of interest was the ligand binding interactions and their stability. As our simulations did not include G proteins or mimetics thereof, we anticipated that the transition of the receptor into its active state was likely beyond our reach, but we hoped to observe such differences between the ligands that could indicate key interactions for agonism.

■ METHODS
Overview. We carried out and analyzed a total of 36 μs of molecular dynamics simulations. We simulated OX 2 receptor without a ligand, with small molecules suvorexant and Nag26 and with two alternative binding modes of orexin-A. Each simulation was run both in a pure POPC (1-palmitoyl-2oleoylphosphatidylcholine) membrane and in a POPC−cholesterol (25%) membrane for 3 μs. For Nag26 and orexin-A simulations, we also produced 1 μs replicas. The simulation conditions are summarized in the Table 1, and the simulations are referred to by their numbers thereafter. For the system setup and analysis we employed Gromacs 5.1, 29 and for the MD simulations Gromacs 4.6.7. 30 Amino Acid Numbering. For clarity, we use one-letter amino acid codes for the orexin-A, and three-letter codes for the OX 2 receptor. For the receptor, we use the GPCRdb numbering scheme, 31 in which the most conserved amino acid in each transmembrane helix is denoted Nx50, where N is the number of the respective helix. Amino acids up-and downstream are numbered consecutively based on their overlap with an ideal α helix; if a residue is "missing", the corresponding number is skipped, and in the case of bulges, the residue furthest from the helical axis is denoted NxYY1, where NxYY is the number for the previous residue.
System Setup. The OX 2 Receptor. The OX 2 receptor crystal 32 (PDB id: 4S0V) is missing the N-terminus upstream of Pro50, residues 160−163 (the second intracellular loop, ICL2), 198−202 (top of the hairpin in the second extracellular loop, ECL2), 255−293 (ICL3), 336−337 (ECL3), and the Cterminus downstream of Cys381. Few side-chain atoms are also missing near the unresolved sections. To fill these gaps, we used Modeller 9.14 33 with default settings to build 30 homology models of the human OX 2 (residues 50−254 and 294−381) using the crystal structure 4S0V as a template, and selected the model with the lowest RMSD in comparison to the template. We left out the receptor termini and the ICL3 since long stretches such as these do not produce reliable models without suitable templates, which were not available at the time. We set all arginine, lysine, glutamic acid, and aspartic acid residues as charged, except for the Asp100 2x50 , which we set to be protonated. 34 We allowed Gromacs to decide the protonation states for histidine residues based on the optimal hydrogen bonding conformation with default settings.
Peptide Docking. Orexin-A 2−33 in the straight α-helical conformation 11 was docked to the receptor model with ZDOCK 35 and docking poses were refined with RDOCK 36 via Discovery Studio 4.5. 37 As pyroglutamate is not parametrized in these docking algorithms, we were forced to omit the first residue. It was added manually prior to simulation. From top-scoring docking poses, we selected manually two poses that were close matches to our previously published binding modes, 21 which were named the TM5-mode and the TM7-mode based on the orientation of H26.
Small Molecule Placement. Suvorexant was copied from the OX 2 crystal structure 4S0V. 32 Nag26 was docked to the crystal structure as described previously. 38 In short, Nag26 was docked with Glide induced-fit protocol, retaining the binding site water molecules. Shape and interaction similarities to  A, TM5  POPC  3  2  orexin-A, TM5  POPC  1  3  orexin-A, TM5  POPC−CHOL  3  4  orexin-A, TM5  POPC−CHOL  1  5  orexin-A, TM7  POPC  3  6 orexin-A, TM7 The Journal of Physical Chemistry B Article suvorexant directed the pose selection. The selected pose (but not the water molecules) was then copied into the binding cavity of the homology model with no further refinement. Membrane. We employed CHARMM-GUI 39−42 to generate two hexagonal membranes solvated in TIP3P 43 water to result in a 12 nm total height. We used a pure POPC membrane of 150 lipids per leaflet, and a POPC + CHOL (25%) membrane of 126 + 42 lipids per leaflet. The membrane size was selected to ensure at least 50 Å between the periodic images of the protein in the membrane plane, and 25 Å in the direction perpendicular to the membrane. Both membrane models employed here are simplified representations of reality. For computational feasibility, we selected only one phospholipid for the study. POPC membranes are frequently used in MD simulations as the lipid is commonly present in biomembranes. 44 The comparison of simulations, especially in the field of GPCRs, is facilitated by employing membranes similar to other studies. 28,45−50 While most MD simulations with GPCRs employ single-component membranes, we additionally replicated all simulations with cholesterol-containing membranes to probe for the effects of the membrane composition. Cholesterol is known to be present in the cellular membrane, and its effects on the membrane thickness and packing are clear. 51−54 Also, several interaction sites for cholesterol have been mapped onto the membrane-facing interface of GPCRs. 55 As the starting positions for the cholesterol in contact with the receptor can have a significant impact on their distribution through the simulation, we carried out coarse-grained MARTINI 56 simulations to define suitable starting locations. We used CHARMM-GUI [39][40][41][42]57 to embed the OX 2 receptor in POPC−cholesterol bilayers with 10%, 25% or 50% of cholesterol, and to derive necessary topologies for the simulation. All systems were energy-minimized for 5000 steps and subsequently equilibrated for 100 ns (NVE ensemble). Then, each system was simulated for 10 μs in triplicate under NPT conditions. A time step of 20 fs was used in the simulations. Temperature of 310 K and pressure of 1.0 bar were maintained with V-rescale thermostat and Berendsen barostat with coupling times of 1.0 ps. Long-range electrostatics were calculated with the reaction field method and the van der Waals interaction was cut off at 1.1 nm. All backbone atoms of the receptor were position restrained with a force constant of 1000 kJ/mol × nm 2 . The first 3 μs of each simulation were treated as equilibration, and the final 7 μs were subjected to analysis of cholesterol density around the protein. This revealed 8 sites for cholesterol molecules, and 15 sites for POPC molecules. These were back-mapped to an all-atom system.
Construction of Simulation Systems. We inserted the receptor complexes into the membranes with the Gromacs tool membed with default settings except for a 0.1 resize factor in the xy-dimension. To the POPC membrane, the protein complex was embedded as such; this resulted in a deletion of 3 + 3 POPC molecules, except for the peptide simulations, where the intracellular leaflet lost an additional lipid. For the POPC−cholesterol membrane, we first combined the receptor with the lipids identified by the coarse-grained simulation. Then, to ensure a symmetrical membrane with the intended 25% cholesterol content, we selected manually a suitable location for the insertion into the generated membrane and created a circular hole by deleting 18 POPC molecules per leaflet and 6 + 4 cholesterol molecules from extra-and intracellular leaflets respectively (as 5 out of 8 identified cholesterols were on the extracellular leaflet). The receptor and the close lipids were then inserted into the pregenerated hole.
In each system, we added NaCl to neutralize the charge for the receptor complex and to reach a salt concentration of 100 mmol/L. The resulting systems contained 293−294 POPC lipids or 246 POPC + 82 CHOL lipids, and ∼25 000 water molecules in addition to necessary ions and the receptor complex.
Force Field and Protocol. We used Amber99sb-ildn force field 58 together with Slipids 59 parameters for POPC and cholesterol. We added parameters for pyroglutamate manually, and used atom charges from TINKER 60 parameters. Suvorexant was parametrized with Antechamber, 61 with Gaussian 62 -derived RESP partial charges. Nag26 parameters were combined from Amber99sb-ildn (most bond and angle parameters), OPSL-AA 63 (bonds and angles whenever native Amber parameters were unavailable, and torsions for the sulfonamide) and Antechamber (torsions excluding the sulfonamide, AM1-BCC charges 64 ). The topologies for the small molecules are available as part of the Supporting Information.
All simulations were run under periodic boundary conditions in hexagonal prism-shaped boxes. For all simulations, the time step was 2 fs, center-of-mass motion was removed every 10 steps for membrane−protein−ligand complex and solvent separately, and Verlet pair-list was updated every 20 steps. Temperature was kept at 310 K with Nose−Hoover thermostat separately for membrane, solvent and the protein complex, and pressure of 1 bar was maintained independently for z-and xy-dimensions with semi-isotropic Parrinnello−Rahman barostat. Short-range electrostatics and van der Waals forces were cut off at 1 nm with long-range dispersion corrections for energy and pressure. PME was used for long-range electrostatics with cubic interpolation and Fourier spacing of 0.12 nm. LINCS constraints were used with all bonds.
Equilibration and Production. We used an equilibration scheme inspired by a previously presented work. 65 Each simulation was equilibrated in seven 10 ns phases; in the first phase we imposed a 1000 kJ/mol restraint on protein and ligand heavy atoms, as well as on cholesterol heavy atoms for the simulations in the cholesterol-containing membrane. Through phases 2−5, we tapered off the heavy atom restraints with 200 kJ/mol intervals, except for the cholesterols, where the restraint was removed after phase 1. At phase 6, the 200 kJ/mol restraint was applied to all protein and peptide Cα, as well as to all heavy atoms in the small molecular ligands. At phase 7, the restraint was lifted from those α carbons that reside in the receptor loops. In addition to visual examination of the receptor, ligand and membrane, we monitored the equilibration through the stabilization of potential energy, simulation box vectors and receptor Cα-RMSD.
The equilibration was followed by 3 μs production runs. For the Nag26 and orexin-A simulations in both membranes, we started additional 1 μs replicas from the equilibrated systems.
Analysis. For analysis, each trajectory was aligned based on the receptor Cα with the Gromacs flag progressive; i.e., each frame was aligned to the previous frame to create continuous motions. Analysis was carried out on these full trajectories, except for the analysis involving water and calculation of the preferred locations for cholesterol around the receptor. For these, trimmed trajectories with frames every 3 ns were used.

Article
The analysis of membrane thickness and area per lipid was carried out on unaligned coordinates of the final frame, as the receptor-based alignment tends to twist the membrane slightly out of the xy-plane.
The analysis of receptor and its interactions was conducted with Gromacs utilities in combination with VMD 66 and Matlab. 67 RMSD for protein and peptide was measured from the α carbons, and for the small molecules from the heavy atoms (i.e., non-hydrogen). We refer to RMSD as "external", when the alignment was based on the receptor Cα, and "internal" when the alignment considered only the ligand Cα (peptides) or heavy atoms (small molecules). Thus, external RMSD contains ligand translation, rotation and conformational changes, whereas internal RMSD describes only the conformational stability. Interactions were mapped residuewise (or by atom groups in small molecules) with gmx mindist with a cutoff of 4 Å, disregarding hydrogens, and hydrogen bonds with the tool gmx hbond with default settings. Pairwise interaction energies were extracted with gmx enemat. Waterbridged hydrogen bonds were scanned with an in-house script by framewise inspection of multiple hydrogen-bonding partners to a single water molecule.
Stable water molecule locations, as well as preferred locations for cholesterol were mapped with VMD Volmap as density with default settings. Cholesterol density across the membrane was mapped using the gmx densmap with grid size of 1 Å. For the density analysis, we aligned the trajectories on the protein, but allowed rotation only around the z-axis to prevent the membrane from tilting, as this would have hindered the separated analysis of the membrane leaflets.
Order parameters for the POPC lipids were calculated with gmx order, and the double bond was treated as described by Pluhackova et al. 68 The area per lipid and bilayer thickness were calculated with GridMAT-MD. 69 For the analysis, the simulation box was surrounded by periodic images in the xyplane, while calculations included only the lipids from the original box. The average membrane thickness was estimated as the area-weighted mean of the headgroup P−P distance at individual lipids. Membrane. First, we analyzed the membranes to verify that our simulations replicate the expected membrane properties. At the end of the simulations, the pure POPC membranes showed an average area per lipid (APL) of 60.9 ± 0.8 Å 2 , and membrane thickness of 36.9 ± 0.4 Å (see Supporting Information, Table S1). These values are in line with experimental data. 51−54 The unsymmetrical insertion of some peptide−receptor complexes did not result in marked differences in APL. The cholesterol containing membranes showed clear effects of cholesterol condensing with POPC APL of 51.1 ± 0.6 Å 2 and membrane thickness of 41.4 ± 0.5 Å. Cholesterol APL was 37.5 ± 0.9 Å 2 . The ordering parameters of the POPC tails closely replicate experimental values 70 for both pure and cholesterol-containing membrane; see Supporting Information, Figure S1. Despite the clear effect of cholesterol on the membrane properties, we observe no consistent differences in ligand binding interactions or receptor mobility between pure POPC membrane and the cholesterolcontaining membrane. This might reflect the orexin receptors' indifference to the membrane composition, or more likely catching such effects would require more replicas or extended simulation periods.
The class A GPCRs have been postulated to feature specific cholesterol-binding sites facing the membrane milieu. 55 On the extracellular leaflet, crystal structures of GPCRs have shown cocrystallized cholesterol molecules next to the TM3, at the TM4−5 interface and along the ECL3-side of the receptor, where the position at the TM6−7 interface has been most frequently populated. 55 On the intracellular leaflet, the most notable site is the suggested "Cholesterol Consensus Motif" (CCM) 71 located at the cleft formed by the TM1−4, which provides a binding site for one or two cholesterol molecules at least in the β 2 -adrenoceptor. As these sites offer a possible route for cholesterol-mediated effects, we sought to find and populate them through the coarse-grained simulations. The simulations suggested five sites on the extracellular leaflet (one facing the TMs 1 and 7, two in the crevice formed by the TM2−4, one close to the TM4−5 and one facing the TM6−7) and three sites on the intracellular leaflet (two next to the TMs 1, 2, and 4, and one facing the TM3−5) ( Figure 2). Many of the sites discovered by the coarse-grained simulations match the cholesterol binding observed in the crystal structures; both sides of the extracellular end of the TM7, at the extracellular TM4−5 interface, and at the intracellular TM3−5 interface. Two cholesterol molecules were located close to the CCM site enclosed by the TMs 1−4, but not completely within the cavity. Of the starting sites, only the extracellular leaflet positions on both sides of the TM7 and next to the TM4−5 were stable through equilibration, and only the two positions flanking the TM7 were consistently stable in the production simulations ( Figure 2, Supporting Information, Figure S2). Individual simulations show also other binding sites, most often in the extracellular leaflet, but these are not consistent across the pool of simulations (Supporting Information, Figure  S2). Interestingly, the TM7-flanking sites, which were the most stable, are among those observed through X-ray crystallog- The Journal of Physical Chemistry B Article raphy. 55 Aromatic residues in the region (Phe333 ECL3 , Trp345 7x33 , Phe348 7x36 , Trp351 7x39 , and Tyr53 1x31 appear to interact with cholesterol in these "pockets". The significance of these sites for orexin receptor function remains to be elucidated. Receptor Conformation. GPCR activation is heralded by the outswing of the intracellular end of the TM6, which is accompanied by the reorganization of interhelical interactions at the intracellular side of the receptor. 25 The interaction patterns have been analyzed for all GPCRs that have been crystallized both in the active and the inactive state. The analysis showed that there are five conserved "inactive" interactions, which in OX 2 would be Ile148 3x46 −Leu306 6x37 , Val75 1 x 5 3 −Tyr364 7 x 5 3 , Tyr364 7 x 5 3 −Phe371 8 x 5 0 and Asn365 7x54 −Arg372 8x51 . 25 Similarly, there are two conserved "active" interactions, which would be Leu310 6x41 −Val240 5x55 and Tyr364 7x53 −Ile148 3x46 in the OX 2 . In addition to these interactions, which locate at the intracellular side of the receptor, a triad of amino acids at the center of the receptor has been suggested to reorganize upon receptor activation. 72 This "core triad" comprises Phe313 6x44 , Pro235 5x50 and Val142 3x40 in the OX 2. During the transition from the inactive to the active state, the side chain of Phe 6x44 is suggested to pass the side chain of the amino acid 3x40 (valine in orexin receptor but isoleucine in adrenoceptors where the core triad reorganization has been observed) and move closer to Pro 5x50 . This movement is also tied to the outswing of the TM6.
Our 16 simulations showed a quite static receptor (receptor Cα RMSD of 2−4 Å in comparison to the equilibrated structure, see Supporting Information, Figure S3). No largescale conformational changes were visible, including the TM6, which retained its closed, inactive conformation. Residuewise RMSF showed that most fluctuation took place at or near the loop regions or chain termini (data not shown). Concerning the reorganization of contacts upon activation, none of the simulations showed either of the two active interactions described above. The five interactions typical to inactive receptors were mostly maintained, and even if they were periodically broken, it was not caused by activation-like conformational changes. The core triad of Phe313 6x44 , Pro235 5x50 and Val142 3x40 remained in the starting (inactive) conformation, except for the simulations 7 (TM7-mode) and 10 (suvorexant). The observed reorganization was not replicated in other simulations with these ligands, and not reflected by the large-scale conformation of the TM6, which remained in the inactive conformation.
The antagonist-bound orexin receptor crystal structures display three to four salt bridges lining the binding site. All structures show the salt bridges Asp 45x51 −Arg 6x59 , Glu 45x52 − His 5x40 , and Asp 2x64 −His 7x38 . In addition, the OX 2 crystal structures show a salt bridge between Glu 2x67 −Arg 7x27 . It has been suggested that these bridges might stabilize the inactive state, and that a disruption of the bridges could lead to receptor activation. 21,32 While the simulations conducted here showed nonuniform salt bridge behavior (see Supporting Information, Figure S4), the differences did not appear to depend on the ligand or the membrane composition.
Peptide Binding Modes. One of the aims of our study was to assess the stability of the two alternative orexin-A binding modes we had previously suggested, 21 called the TM5mode and the TM7-mode based on the orientation of the (arbitrarily selected) H26. In the previous study, neither binding mode was promoted by the assessment of docking scores and interactions toward residues that had been shown important by SDM.
While the 19-amino-acid fragment from the orexin-A Cterminus has been shown to suffice for the biological activity, we elected to use a full-length peptide, which we docked in the straight α-helical conformation to replicate our previous docking poses. The simulations, however, showed that this straight conformation was not stable unless hydrogen bonding between the N-terminal hook and the receptor ECL2 took place. These hydrogen bonds often involved the ECL2 main chain, and were not consistent across the simulations. In absence of hydrogen bonding between the orexin-A Nterminus and the receptor ECL2, the peptide hinge was allowed to bend, moving the peptide N-terminus to rest horizontally atop the receptor. This bending took place in the simulations 2 and 5, which feature the TM5-and TM7-modes, respectively, while other simulations with orexin-A showed the N-terminus in contact with the ECL2, following its fluctuations (Figure 3). Our simulations did not include the receptor Nterminus, as its structure and location had not been resolved at the time. The possible impact of the receptor N-terminal helix 73,74 on the peptide conformation and binding is discussed below, after the description of the peptide C-terminus binding. The Journal of Physical Chemistry B

Article
As the C-terminus of the peptide is the key for biological activity, we focused in its stability and interactions. To assess the stability of the predicted binding modes, we analyzed the peptide C-terminus (N25−L33) separately, since the bending of the hinge region and the conformational flexibility of the Nterminal hook dominate the RMSD of the full-length peptide. Overall, the receptor-bound C-terminus of orexin-A maintained its conformation and location better in the TM5-mode (i.e., H26 side chain points toward TM5) than in the TM7mode, and for both binding modes, the conformation was more stable in the pure POPC membrane than in the cholesterol-containing membrane ( Figure 4).
As the binding interaction of the peptide N-terminus are heavily influenced by the (possible) bending of the peptide, we assessed the binding interactions only for the peptide Cterminus (N25−L33). These interactions, averaged across the simulations, are shown in Figure 5, panels A and B for the TM5-and TM7-mode, respectively. Simulationwise heatmaps of the interaction frequencies, along with mapped interaction energies are available in the Supporting Information as Figures S5−S12. The stable location and conformation of the peptide C-terminus in the TM5-mode is mirrored in well-defined, stable interactions across simulations 1−4, whereas the higher mobility and instability of the TM7-mode is reflected by more variable interactions in the simulations 5−8.
As described, the orexin-A C-terminus remained quite stable in the TM5-mode (simulations 1−4), especially in the POPC membrane (1−2). In the cholesterol-containing membrane (3−4), the peptide rose some 1.4 Å further from the bottom of the binding cavity, which was reflected in the external RMSD. In the simulation 4, the peptide C-terminus also tilted 10°t oward the TM1. Despite these differences, the interactions remained similar in all four simulations ( Figure 6). N25 interacted mainly with residues at the top of the TM7, most often Tyr343 7x31 . In simulations 3 and 4, this interaction was often a hydrogen bond; in the simulation 3, between the side chains, either NH 2 ···OH or OH···O, and in the simulation 4, due to the tilting of the peptide, Tyr343 7x31 OH to N24 mainchain carbonyl. The simulation 1 showed an alternative interaction, namely a hydrogen-π bond from N25 NH 2 to Phe346 7x34 . In all four simulations, the N25 also interacted with His350 7x38 , either directly as seen in the simulation 1 and 2 (POPC membrane), or via a water-mediated hydrogen bond due to the increased distance in the simulation 3 and 4 (cholesterol membrane). H26 interacted with the side chains of the acidic residues in the ECL2 (Asp211 45x51 and Glu212 45x52 ) in the simulations 1 and 2 (a prominent interaction also in terms of interaction energy, as seen in Figures S5−S6 in the Supporting Information), while the peptide repositioning seen in the simulation 3 and 4 was accompanied by a local tightening of the helical turn, which brought the H26 side chain toward the ECL3. Further toward the C-terminus, the interactions were similar in all simulations; A28 side chain faced a pocket lined by His350 7x38 , Thr111 2x60 , Val114 2x63 , and the aliphatic chain of Asp115 2x64 . G29 allowed the close packing of His350 7x38 against the side of the peptide and the A28−G29 amide bond. I30 faced Met191 4x65 and Pro131 3x29 . To the TM2-side of the Pro131 3x29 lay the L31 side chain in a pocket lined also by Ile130 3x28 , Ala110 2x59 , Thr111 2x60 , and the aliphatic chain of Gln134 3x32 . The T32 side chain faced Tyr354 7x42 , Val353 7x41 , His350 7x38 , Ile320 6x51 , and Thr111 2x60 . Most often it hydrogen bonded intramolecularly to A28 (OH···O), unless the His350 7x38 adopted a rotamer where the unprotonated Nδ became available for hydrogen bonding. This took place in the simulations 3 and 4, where the His350 7x38 side chain flipped upward to follow the N25, in concert with the peptide movement. The L33 side chain lay in a pocket lined by Ile320 6x51 , Asn324 6x55 , and The Journal of Physical Chemistry B Article Phe227 5x43 . The NH 2 of the amidated C-terminus mainly bonded to the G29 or the I30 carbonyl. An interesting network of interactions took place at the bottom of the binding cavity, between the exposed carbonyls of L31 and T32, along with Gln134 3x32 and Tyr317 6x48 . In the beginning of the simulations, Gln134 3x32 side-chain nitrogen usually formed a hydrogen bond to the carbonyl of T32, while Tyr317 6x48 hydroxyl pointed down. This was most likely imposed by the starting configuration. During the simulations, the Gln134 3x32 hydrogen bond shifted to the carbonyl of L31, Tyr317 6x48 hydroxyl flipped upward, and a water molecule moved in to bridge interactions between the Gln134 3x32 side-chain carbonyl,  The Journal of Physical Chemistry B Article Tyr317 6x48 hydroxyl and the T32 carbonyl. The Gln 3x32 −L31 interaction was also prominent in terms of interaction energy ( Figures S5−S8 in Supporting Information). In the simulation 2, this reorganization did not happen, as Tyr317 6x48 hydroxyl immediately bound to the Gln134 3x32 side-chain carbonyl, leaving no space for the bridging water molecule. While the hydroxyl of the Tyr354 7x42 was consistently close to the carbonyl of L31, it participated in the network only transiently, preferring to bind with Thr111 2x60 ; this interaction was almost exclusively seen only in the TM5-mode (simulations 1−4). The carbonyl of L33 also took part only rarely.
The instability of the orexin-A helix in the TM7-mode rendered it difficult to sum up the interactions across the simulations 5−8. The simulation 7 showed an almost complete "melting" of the peptide C-terminus, and in simulation 6 and 8 the hinge region unfolded and fell into the cavity between the TMs 2 and 7. The common observations between the simulations were that N25 was close to Asp211 45x51 and Cys210 45x50 , I30 inhabited a pocket lined by Ile320 6x51 , Asn324 6x55 , Lys327 6x58 , and Phe346 7x34 , and L33 was close to His350 7x38 , often hydrogen bonded to the C-terminal amide. In terms of interaction energy, there was no apparent consistency between the simulations (Supporting Information, Figures S9− S12). The simulation 7 with the most prominent peptide unfolding aside, the three other simulations (5, 6, and 8) showed the side chain of A27 facing a salt bridge formed by Glu212 45x52 and Arg328 6x59 , A28 residing under the ECL2 and the side chain of L31 facing the pocket of Pro131 3x29 , Thr135 3x33 , Gln187 4x61 and Met191 4x65 . The C-terminal carbonyls did not form uniform interactions such as those described above for the TM5-mode.
Given the difference in stability, in terms of both peptide conformation and binding interactions, between the TM5-and the TM7-mode of binding, the TM7-mode seems less likely to be the correct binding mode. Therefore, if the C-terminus of the orexin-A indeed binds in an α-helical conformation, we suggest the TM5-mode as the binding mode for orexin-A at the OX 2 . The binding mode should also be valid for the OX 1 receptor and for orexin-B, since the orexin peptide C-terminus is conserved and the receptor structures are similar. However, with current data, we cannot rule out the possibility of nonhelical bioactive conformation. Endothelin, for example, adopts a conformation that is mostly helical but the Cterminus is extended into the receptor cavity. 75 On the other hand, endothelin displays a similar extended C-terminus also in solution, 76 whereas orexin peptides show consistently helices also at the peptide C-terminus.
Concerning the receptor amino acids that SDM studies have shown to be important for orexin-A binding and/or receptor activation, 14,15 the peptide C-terminus in the TM5-mode interacts strongly with Thr111 2x60 , Gln134 3x32 , and His350 7x38 , comes into direct close contact with Phe227 5x43 , Asn324 6x55 and Tyr354 7x42 and forms water-bridged interactions to Tyr317 6x48 . Thr111 2x60 and Tyr354 7x42 are also consistently hydrogen bonded with each other. Asp211 45x51 and Glu212 45x52 in the ECL2 interact with the middle section of the orexin-A, but defined binding partners cannot be summarized from the simulations as the peptide flexibility results in quite different midsection conformations and locations. Trp214 in the ECL2 and Tyr223 5x39 line the binding pocket at the ECL2−TM5 junction, but the peptide does not interact with them in our simulations. However, both pack between the TMs 4 and 5, which suggest that they might contribute to the receptor overall conformation or the binding site conformation. Thr231 5x461 and Tyr232 5x47 are mostly out of reach for the peptide, especially as Tyr232 5x47 mostly resides between the TMs 5 and 6. There are also no frequent water-mediated interactions with the peptide, but the residues could be a part of the hydrogen-bond network at the bottom of the binding site or otherwise important for the activation cascade.
The conformations of the bound peptide in the TM5-mode simulations reflect the NMR-derived models of the conformation in water; the peptide C-terminus and the middle section retain their helical conformation, while the hinge region allows the peptide to bend. However, in the simulation 2, and to some extent also in the simulation 3, the peptide bends toward the TM5, which is not among the bending directions observed in the NMR models of orexin-A. 11,12 The direction is similar to the orexin-B solution structure, 13 except that the bending takes place one helical turn further from the C-terminus. Recent studies 73,74 have been able to resolve the structure of the OX 1 and OX 2 N-terminus, which appears to form a short amphipathic helix either atop the receptor or pointing away from the receptor, residing atop the membrane. While it remains unclear whether these different locations represent different functional states or if either one is a crystallization artifact, the N-terminal helix was shown, through mutagenesis, to be important in the peptide-induced receptor activation. 73 This suggests direct interaction with the peptide, which in turn would promote the location next to the ECL2. If the N-terminal helix would indeed lie atop the receptor in interaction with ECL2, as seen in the OX 1 structure, 73 it would prevent the binding of orexin-A in the straight α-helical conformation such as the one we used for docking and as the starting conformation for the simulations. However, the peptide bent toward the TM5 would fit snugly underneath the N-terminal helix, presenting R15 and Y17 among others toward the downward-facing polar residues of the N-terminal helix (Figure 7). Interestingly though, the N-terminal helix in the crystal structure features lipophilic residues (Phe31, Leu32, Tyr34, Leu35, and Trp36) facing toward the solvent, which seems unfavorable. A roll of 180°would bring them in contact with L16, Y17 and L20 in the orexin-A. It is also possible that the N-terminal helix adopts a distinct peptide-bound Small Molecule Binding. In addition to orexin-A, we carried out simulations with bound small molecules, the antagonist suvorexant (simulations 9−10) and the agonist Nag26 (simulations 11−14). Suvorexant binding mode seen in the crystal structures was stable through the simulations in both membranes. The external RMSD fluctuated around 1 Å, and the internal RMSD between 0.5−1 Å (Figure 8). The chlorinated benzoxazole remained in the pocket lined by Cys107 2x56 , Ala110 2x59 , Thr111 2x60 , Val114 2x63 , Trp120 23x50 , Ile130 3x28 , Pro131 3x29 , and Gln134 3x32 , with the oxazole ring Figure 8. RMSD of the small molecular ligands. In green and orange, the "external" RMSD, and in blue and red, the "internal" RMSD. Figure 9. Heatmap of binding interactions of the small molecular ligands. For the mapping, the small molecules were divided into groups as displayed below the heatmaps. The heatmaps display the fraction of simulation time that any interatomic distance was below 4 Å. H:xx denotes the presence of a hydrogen bond and W:xx the presence of a water-mediated hydrogen bond. These are shown only when the frequency was over 20%.

The Journal of Physical Chemistry B
Article packing next to Pro131 3x29 ( Figure 9A). A water molecule was often seen binding with the oxazole nitrogen, but it did not form a stable water-mediated interaction. The sevenmembered ring occupies a central location with the methyl facing Phe227 5x43 . The linker carbonyl hydrogen bonds to Asn324 6x55 and via water to His350 7x38 . The methylated benzene ring interacts with His350 7x38 , Val353 7x41 , and Tyr354 7x42 in TM7, and is flanked by Gln134 3x32 and the benzoxazole ring on the other side. Underneath, the Gln134 3x32 −Tyr354 7x42 hydrogen bond is maintained throughout the simulations, which is not seen consistently in any other simulations. The triazole ring is next to Ile320 6x51 and Phe227 5x43 , and traps 2−3 water molecules into a pocket between the TMs 5 and 6, where the water molecules form a hydrogen-bond network between Ser321 6x52 , Thr231 5x461 , and His224 5x40 .
Concerning the pool of crystallized orexin receptor antagonists, the binding mode of suvorexant is near identical between the OX 1 and OX 2 subtypes, 32,73 and closely replicated by the OX 1 -selective SB-674042. 73 The OX 2 -selective EMPA also occupies a similar location and displays a horseshoe-like conformation. 74 All antagonists display the intact Gln 3x32 − Tyr 7x42 hydrogen bond in the crystal structures.
The small molecular agonist Nag26 was more mobile than suvorexant both in terms of conformation (internal RMSD 1.5−3 Å) and location (external RMSD 3−6 Å) (Figure 8), which is not surprising given the number of rotatable bonds in the ligand scaffold. While the molecule mostly retained the initial location, especially the extremities (groups Na1−2 and Na8−9 in Figure 9), showed significant movement. This was reflected in the interaction patterns, which were not always consistent between simulations. The average binding interactions are shown in Figure 9B, and the interaction heatmaps for individual simulations, along with interaction energies, are available in Supporting Information as Figures S13−S18. In some simulations, the ligand retained its initial horseshoepacked conformation, whereas in others, the ligand extended within the binding site ( Figure 10). The dimethylamide group was initially located deep in the binding cavity, next to Thr231 5x461 and flanked by Tyr317 6x48 , Val138 3x36 , and Ser321 6x52 . At this site, the carbonyl oxygen often formed a water-mediated hydrogen bond with the side chain of Thr231 5x461 . In the simulations 13 and 14 (cholesterol membrane), the group retained this position, whereas in the simulations 11 and 12, the group rose a bit and then pushed in-between the TMs 4 and 5. This rise brought the amide oxygen in water-mediated contact with Asn324 6x55 or Arg328 6x59 . The push also moved Thr231 5x461 and Phe227 5x43 toward the TM4, resulting in a marked counterclockwise rotation of the upper section of the TM5 and tightening of the helical turn at the bulge near Tyr232 5x47 . The sulfonamide moiety, which is crucial for activity, 7 lay next to Pro131 3x29 in all simulations but found different binding partners. In the simulation 13, where the Nag26 conformation was the most stable, the sulfonamide formed both direct and water-bridged hydrogen bonds from the oxygens to Gln134 3x32 . Despite close proximity, Gln187 4x61 and Thr135 3x33 did not hydrogen bond with the sulfonamide. The amide nitrogen formed a stable water-bridge to Glu212 45x52 , which in turn remained bonded to Arg328 6x59 through the simulation. The simulation 11 showed a more mobile sulfonamide moiety, which hydrogen-bonded first with Arg328 6x59 and then moved higher and toward the center of the binding site and switched to Lys327 6x58 . This movement also allowed the amide nitrogen to bond with Asp211 45x51 . Arg328 6x59 and Glu212 45x52 remained in interaction through the simulation, but the ligand movement moved both side chains toward the TM5 and appeared to constrict the entrance of the binding site. Simulations 12 and 14 showed Arg328 6x59 hydrogen bonding with the sulfonamide and Glu212 45x52 /Asp211 45x51 , while the amide nitrogen interacted with Cys210 45x50 main-chain oxygen via water. Water also bridged transient interactions between the sulfonamide oxygens and the side chains of Asp211 45x51 and Glu212 45x52 . The hydrogen bonding between the sulfonamide nitrogen and Lys327 6x58 /Arg328 6x59 seen in simulations 11, 12 and 14 was also a strong contributor to the interaction energy (Supporting Information, Figures S13, S14, and S16). The rest of the molecule was quite mobile; if the ligand retained or readopted the packed conformation, the methylated ring folded in-between the TM7 and the rest of the ligand, favoring interactions with Phe346 7x34 . In the case of ligand extension, the ring moved close to the ECL1, which brought the polar atoms of the chain linker in hydrogen-bonding distance of the exposed carbonyl of Cys210 45x50 , and also Thr111 2x60 .
It is noteworthy that we employed only one starting conformation for the Nag26 agonist. As also our peptide simulations highlighted, the starting conformation can heavily impact the simulation trajectory, and it is possible that another starting conformation for the Nag26 would have yielded different results. Then again, our simulations showed the Nag26 sampling multiple conformations within the binding site, unable to find a markedly stable binding mode. This could mirror the weaker binding interactions in comparison to suvorexant; the Nag26 K i is 140 nM at the OX 2 , while the corresponding value for suvorexant is 0.35 nM. 2,7 Alternatively, the initial conformation might be far from correct, but given the resemblance between our docking pose and the structurally related small molecular antagonist EMPA (not shown), that seems unlikely.
Water Molecules. The simulations showed distinct preferences for water molecule locations. Within the receptor core, Asp100 2x50 was often accompanied by few water molecules ( Figure 11). In some simulations, a semistable chain of water molecules connected these waters to the binding cavity. Two to three water molecules usually resided between the TMs 5 and 6, bound to Thr231 5x461 and Ser321 6x52 , sometimes also to His224 5x40 and Asn324 6x55 . Often, the The Journal of Physical Chemistry B Article simulations showed also stable water molecules that took part in the ligand−receptor interactions; these were discussed in detail above for each complex.
Many of the observed stable water molecules are also present in the high-resolution crystal structures, 32,73,74 such as between the TMs 5 and 6 (Thr231 5x461 and Ser321 6x52 ), around Asp100 2x50 , at the TM6 kink and atop the TM6, and in the intracellular cavity between TMs 3, 4, and 5. Water molecules were also seen around the ECL2 in matching positions to the crystal structures, and close to the TMs 5 and 6. Interestingly, the water molecule that takes part in suvorexant binding, bridging the carbonyl oxygen to His350 7x38 , was also stable in simulation 1 (orexin-A in the TM5-mode), where the water molecule bridged a three-way interaction between Asn324 6x55 , His350 7x38 and the main-chain carbonyl of N25. In simulation 13, the carbonyl oxygen in the linker chain of the Nag26 occupied a similar position. The water molecules between the TMs 5 and 6 provide the most interesting case; there were two-to-three water molecules in all simulations, unless the TM5 made the counterclockwise turn, moving Thr231 5x461 further from Ser321 6x52 and thus breaking the water-binding pocket. With the Nag26 simulations, this was linked with the dimethylamide pushing between the TMs 4 and 5, but in the simulation 2 (orexin-A in the TM5-mode), this turning appeared to take place without clearly defined ligand-dependency.
Activation. The activation of a GPCR is a chain of events propagating through the receptor from the ligand-binding site to the G-protein-binding site at the intracellular face. In our simulations, the intracellular conformation of the receptor remained in the inactive state, as anticipated, which was mirrored by the absence of the common contacts observed in the structures of activated GPCRs. 25 Concerning the transitional steps along the activation cascade, the core triad, 72 which lies underneath the binding site and is thought to be relevant in relaying the activation signal, also remains in the initial, inactive conformation. 72 It is noteworthy that orexin receptors have a valine residue at 3x40, whereas a bulkier residue such as isoleucine is often found in the receptors that display a clear "activation" of the core triad. The M 2 receptor also features Val 3x40 , but there both active 77 and inactive 78 receptor structure resemble the activated core triad of e.g. β 2 receptor. At the binding site, the activation determinants differ between receptors, rendering comparison of orexin receptors to the available structures of active GPCRs dubious. However, a common theme in activating interactions is a change in interhelical interactions across the binding site. With this in mind, we noted that suvorexant simulations showed a consistent hydrogen bond between Gln134 3x32 and Tyr354 7x42 , which was not seen in other simulations, but is replicated by other crystallized antagonists. 73,74 In contrast, the simulations with orexin-A in the TM5-mode showed this hydrogen bond only 0−8% of the simulation time. On the other hand, the TM5-mode propagated Tyr354 7x42 − Thr111 2x60 hydrogen bonding. We speculate that enforcing the Gln134 3x32 −Tyr354 7x42 hydrogen bond could stabilize the inactive state of the receptor, and that the alternative Tyr354 7x42 −Thr111 2x60 hydrogen bond could enable receptor activation, even though Nag26 did not show as strong effect as the orexin peptide did. The comparison between suvorexant and Nag26 highlights the agonist interacting with Thr231 5x461 , a location which is central to the activation of adrenoceptors, for example. The Nag26 also interacted with the ECL2, and the salt bridge Glu212 45x52 −Arg328 6x59 , which crosses the The Journal of Physical Chemistry B Article binding site. One could speculate that this interaction with the salt bridge could induce an inward tilt of the extracellular end of the TM6, which, in turn, could pivot the TM6 and open the intracellular binding site. Nag26 was also seen in two simulations pushing between the TMs 4 and 5, which induced a tightening in the TM5 helix and changed the interactions between the TMs 4, 5, and 6. This could also serve in destabilizing the inactive receptor conformation, thus facilitating receptor activation.

■ CONCLUSIONS
We performed extensive molecular dynamics simulations on the OX 2 receptor. On the basis of the simulations, we suggest that orexin-A binds the receptor in a configuration where N25 faces the TM7 (Tyr343 7x31 ), H26 points toward the TM5 but interacts either with the ECL2 or the ECL3, A27 faces the ECL1, A28 and G29 pack next to His350 7x38 , I30 and L31 flank the TM3 and Pro131 3x29 with the L31 side chain occupying a pocked beneath the ICL1, T32 faces the TM7 and L33 the TMs 5 and 6. The exposed carbonyls in the last helical turn participate in a hydrogen-bonding network formed by a stable water molecule, Gln134 3x32 and Tyr317 6x48 . The peptide N-terminus interacts with the ECL2 hairpin, if extended, or bends toward the TM5 to lie atop the binding site. The presented binding mode is, however, contingent on the assumption that the bioactive conformation of the orexin peptide C-terminus is helical and not unfolded.
The small molecule simulations highlight the stability of suvorexant binding, and the contrasting flexibility of the Nag26 agonist. The antagonist and agonist have such different interactions patterns that it is difficult to identify definite key interactions for agonism or antagonism, but certain differences arise, such as the suvorexant-induced stabilization of Gln134 3x32 −Tyr354 7x42 hydrogen bond and the binding interactions of Nag26 toward Thr231 5x461 , Tyr317 6x48 , and the ECL2 amino acids Cys210 45x50 , Asp211 45x51 , and Glu212 45x52 . Nag26 is also observed pushing the extracellular end of the TM5 into a counterclockwise rotation, which influences the interhelical interactions between the TMs 4, 5, and 6, for example, breaking up a water-binding site. Alex Bunker: 0000-0002-1236-9513 Henri Xhaard: 0000-0002-3000-7858

Author Contributions
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
This work was supported by the Academy of Finland (Grant #287963), the Finnish Cultural Foundation, the Finnish Pharmaceutical Society, and the Orion Research Foundation. We thank the Drug Discovery and Chemical Biology network for organizing computational infrastructure in the HX group and CSC−IT Center for Science, Finland, for generous computational resources. The COST action is thanked for organizing European networks on GPCRs (CM1207/GLIS-TEN and CA18133/ERNEST) and multi-target binding (Mu.Ta.Lig/CA15135).