Free-Energy Simulations Support a Lipophilic Binding Route for Melatonin Receptors

The effects of the neurohormone melatonin are mediated by the activation of the GPCRs MT1 and MT2 in a variety of tissues. Crystal structures suggest ligand access to the orthosteric binding site of MT1 and MT2 receptors through a lateral channel between transmembrane (TM) helices IV and V. We investigated the feasibility of this lipophilic entry route for 2-iodomelatonin, a nonselective agonist with a slower dissociation rate from the MT2 receptor, applying enhanced sampling simulations and free-energy calculations. 2-Iodomelatonin unbinding was investigated with steered molecular dynamics simulations which revealed different trajectories passing through the gap between TM helices IV and V for both receptors. For one of these unbinding trajectories from the MT1 receptor, an umbrella-sampling protocol with path-collective variables provided a calculated energy barrier consistent with the experimental dissociation rate. The side-chain flexibility of Tyr5.38 was significantly different in the two receptor subtypes, as assessed by metadynamics simulations, and during ligand unbinding it frequently assumes an open conformation in the MT1 but not in the MT2 receptor, favoring 2-iodomelatonin egress. Taken together, our simulations are consistent with the possibility that the gap between TM IV and V is a way of connecting the orthosteric binding site and the membrane core for lipophilic melatonin receptor ligands. Our simulations also suggest that the open state of Tyr5.38 generates a small pocket on the surface of MT1 receptor, which could participate in the recognition of MT1-selective ligands and may be exploited in the design of new selective compounds.


S1. Docking of compound 1 in the MT1 receptor
The crystal structure of the MT1 receptor in complex with 2-phenylmelatonin (PDB id 6ME3) 1 was used for docking calculations after deletion of the co-crystallized ligand. The docking grid was centered on 2-phenylmelatonin, defining a bounding box of 20x10x10 Å and an enclosing box of 40x30x30 Å. Ligand docking was performed with Glide 7.9 2,3 in standard precision mode and with default scaling factors, 2 setting MAXKEEP and MAXREF parameters, which control the number of poses to retain after the rough scoring stage and the number of poses to refine, to 50,000 and 4000, respectively. The scoring window cutoff was increased to 1000 to widen the selection of the initial poses for the rough scoring stage. The best ligand pose according to the GScore scoring function was merged into the protein structure and the complex ( Figure S1) was energy-minimized to a gradient of 0.01 kJ·mol −1 ·Å −1 with OPLS3e force field 4 implemented in MacroModel 12.0 5 in an implicit water solvation model, 6 using the Polak-Ribière conjugate gradient method. 7 Figure S1. Docking pose of compound 1 in the MT1 receptor. The naphthalene nucleus resides in the orthosteric binding site, the amide side chain and the methoxy oxygen interact with Gln181 ECL2 and S3 Asn162 4.60 , respectively, as observed for co-crystallized ligands. 1 The biphenyl-carboxylate moiety is accommodated at the protein-membrane interface and protrudes through the channel between TM helices IV and V via the oxybutyloxy spacer.

S2. Minimization and equilibration protocol
The systems were minimized using the sander module implemented in Amber16 8 with three minimization cycles of 5,000 steps, each one passing from the steepest descent method to the conjugate gradient, and with protein backbone restrained with a constant of 50 kcal·mol −1 ·Å −2 . The coordinates and the topology files of the minimized structures were then converted in the Gromacs coordinates and topology input files with the Acpype parser. 9 Systems equilibration of the minimized complexes comprised different stages:

S4. Unbinding CV for steered molecular dynamics: definition of plane η
Plane η is defined through the coordinates of three centers of mass, H1, H2 and H3, defined in counterclockwise order, observing from the outside of the receptor, as in Figure 2 of the main text.
The centers of mass are related to the following alpha carbons:

S5. Path CVs (PCVs) optimization: technical details
Path optimization 15 was performed through four consecutive steered molecular dynamics (SMD) simulations 16

Qualitative agreement of the US windows a priori probabilities
The global PMF, given by the contributions of all the US windows employed in the WHAM calculation, is a free-energy profile over the coordinate , divided in intervals j: where the constant F is fixed to have PMF( =1) = 0 kcal·mol −1 .
The reliability of PMF was evaluated through the qualitative agreement between a priori probabilities of neighboring US windows ( Figure S5). The a priori probability 0 is related to the biased probability that the i th simulation falls into the j th interval according to the following relationship: where represents the probability distortion due to the bias potential and is a normalization coefficient. Figure S5. Dissection of the global potential of mean force (PMF) into the free-energy profiles of the PCV-US windows, scaled as free-energy differences with respect to =1 (bound state). Tyr187/200 5.38   Table S1. Parameters adopted for monodimensional and bidimensional well-tempered metadynamics simulations on Tyr187/200 5.38 dihedral angles.