Binding, Thermodynamics, and Selectivity of a Non-peptide Antagonist to the Melanocortin-4 Receptor.

The melanocortin-4 receptor (MC4R) is a potential drug target for treatment of obesity, anxiety, depression, and sexual dysfunction. Crystal structures for MC4R are not yet available, which has hindered successful structure-based drug design. Using microsecond-scale molecular-dynamics simulations, we have investigated selective binding of the non-peptide antagonist MCL0129 to a homology model of human MC4R (hMC4R). This approach revealed that, at the end of a multi-step binding process, MCL0129 spontaneously adopts a binding mode in which it blocks the agonistic-binding site. This binding mode was confirmed in subsequent metadynamics simulations, which gave an affinity for human hMC4R that matches the experimentally determined value. Extending our simulations of MCL0129 binding to hMC1R and hMC3R, we find that receptor subtype selectivity for hMC4R depends on few amino acids located in various structural elements of the receptor. These insights may support rational drug design targeting the melanocortin systems.

Introduction

The agonistic proopiomelanocortin-derived neuropeptides α-, β-, and γ-melanocyte-stimulating hormone (α-, β-, γ-MSH), the adrenocorticotropic hormone (ACTH) and the inverse agonistic agouti-related peptide (AgRP) are ligands for five closely related melanocortin receptors (MC1-5R), a subgroup of class A G-protein coupled receptors (GPCRs) (Gantz et al., 1993; Cone, 2005, 2006; Dores et al., 2014). The melanocortin system is involved in a plethora of physiological functions, including pigmentation, steroidogenesis, energy homeostasis, and sexual function (Cone, 2006; Rodrigues et al., 2015).

MC1R is expressed by the dermal melanocyte and regulates the synthesis of eumelanin (black–brown pigment) versus phaeomelanin (yellow–red pigment) (Herraiz et al., 2017). The ligand ACTH acts on MC2R (Liang et al., 2013) and is involved in controlling glucocorticoid levels (Fridmanis et al., 2017). The exact physiological role of the MC3R is not clear. Potential roles include control of energy homeostasis in addition to MC4R or immune response and circadian rhythms (Girardet et al., 2014; Mountjoy, 2015). Mice with targeted deletion of MC5R exhibit a defect in water and thermoregulation (Chen et al., 1997), and recent studies have suggested MC5R/α-MSH involvement in modulation of muscle glucose uptake and thermogenesis (Enriori et al., 2016).

The MC4R is of specific interest in the central melanocortin pathways, because it is known to regulate food intake, glucose homeostasis, and energy expenditure (Yang and Tao, 2017). The huge number of pathogenic heterozygous inactivating MC4R mutations are the most frequent genetic cause of obesity in humans (Farooqi et al., 2003). Consequently, MC4R agonists have been developed as potential therapeutics for obesity (Kievit et al., 2013; Chen et al., 2015; Kuhnen et al., 2016) and have been proven effective in treatment of male and female sexual dysfunctions, whereas MC4R antagonists are potential anxiolytics that may function as antidepressants or for treating cancer-induced anorexia (Chaki et al., 2003; Ercil et al., 2005; Pontillo et al., 2005; Chaki and Okubo, 2007; Lansdell et al., 2010).

Although several clinical trials with MC4R ligands have been initiated (reviewed in Ericson et al., 2017b), none of these ligands has yet been therapeutically approved, mainly because of adverse side effects related to low subtype or signaling selectivity (Yang and Harmon, 2017; Yang and Tao, 2017). For instance, the non-selective agonist bremelanotide has been shown to induce elevated blood pressure in clinical trial for male and female sexual dysfunction, while the MC4R selective peptide agonist setmelanotide (10-fold over MC3R) have shown no adverse effects on either blood pressure or heart rate in an anti-obesity study (Ju et al., 2018). However, setmelanotide treatment has impact on skin or hair coloring (Kuhnen et al., 2016). Additionally, MC4R agonist melanotan-II, a super-potent cyclic melanotropic peptide, increased blood pressure (Dorr et al., 1996). Animal models have revealed that MC4R and not MC3R, to be the receptor subtype responsible for stress-induced anorexia (Chaki and Okubo, 2007) and more recently an animal study provided a preclinical proof of concept of the HS014 selective antagonist in treating post-traumatic stress disorder (Sabban and Serova, 2018).

Generally, MC4R-targeting molecules can be linear or cyclic peptides, or small molecules (Goncalves et al., 2018). A comprehensive understanding of molecules that act at MC4R may help to design safe and effective anti-obesity drugs (Ju et al., 2018). Recent medicinal chemistry efforts have focused on using cyclized peptides, cyclotides, and constrained tetrapeptides to improve bioavailability, selectivity, and blood-brain-barrier uptake (reviewed in Zhou and Cai, 2017), whereas more recent small-molecule ligands development has focused on screening GPCR-privileged pharmacophores (reviewed in Ericson et al., 2017a).

Such drug-development problems have been solved for other GPCRs with the aid of robust structural information and structure-based drug design (SBDD) (Sabban and Serova, 2018). Whilst the last decade has witnessed a revolution in the field of GPCR structural biology (Granier and Kobilka, 2012), structural information regarding the melanocortin receptors and ligands has been limited to NMR data of AgRP or standalone transmembrane helix 2 (TM2) of the receptor (Sabban and Serova, 2018). Biomolecular [molecular dynamics (MDs)] simulations have been successful in providing missing structural information for GPCRs (Dror et al., 2011, 2012, 2013; Grossfield, 2011; Kruse et al., 2012; Vanni and Rothlisberger, 2012; Rose et al., 2014; Clark, 2017). Moreover, homology models based on high-quality templates in combination with microsecond-scale unbiased MD refinement (Clark, 2017) have provided robust experimentally testable structural and thermodynamic information on GPCRs (Milanos et al., 2016; Saleh et al., 2016, 2017a,b,c).

Using an established metadynamics protocol, we have been able to predict GPCR-ligand binding -affinities, -modes and free-energies accurately (Saleh et al., 2017b). Metadynamics simulations use a history-dependent bias to encourage the system to explore other areas of the free-energy landscape that are not accessible using plain MD simulation, even using special purpose supercomputers (Laio and Gervasio, 2008). Here, we have used microsecond-scale molecular-dynamics and metadynamics enhanced sampling to investigate binding mechanisms and thermodynamics of the selective non-peptide antagonist MCL0129 to human hMC4R.

These insights were used to unravel the molecular determinants of MCR subtype selectivity and the antagonistic properties of MCL0129.

Materials and Methods

Homology Modeling of the Initial hMC4R Conformation in an Inactive State

Human MC4R (hMC4R) is characterized by several specific properties with respect to the amino acid constitution and related structural features. In particular:

an extremely short second extracellular loop 2 (ECL2, constituted by approximately four amino acids),

a missing cysteine disulfide bridge between the ECL2 and TM3 which is highly conserved among GPCRs,

a regular α-helical conformation of TM5 because of a methionine instead of a proline (Pro5.50) that is highly conserved in class A GPCRs (Supplementary Figure S1) and induces a helical kink and bulge,

the highly conserved class A GPCR motif 7.49xxY in TM7 is a 7.49xxY motif,

a disulfide bridge in the third extracellular loop (ECL3) (Tarnow et al., 2003).

Because these features are of functional importance, we searched for a structural template that best represents the criteria described. Using the GPCR-SSFE 2.0 modeling server (Worth et al., 2017), we found the inactive-state structure of the lyso-phospholipid sphingosine 1-phosphate receptor [S1PR1, PDB entry 3V2Y (Hanson et al., 2012)] with a sequence similarity in the transmembrane region of 56% and a sequence identity to hMC4R of 30% (Blossum62 matrix, Supplementary Figure S1). Such high sequence similarity between the suggested template and MC4R is important for the accuracy of the structural model (Costanzi et al., 2016). The template structure is characterized by a leucine at position 5.50 in TM5 (Ballesteros and Weinstein, 1995) and in consequence shows a regular α-helical conformation, as also expected for hMC4R. Moreover, the disulfide bridge between ECL2 and TM3 is missing but a disulfide bridge is located in the position proposed for MC4R (Sabban and Serova, 2018) (Supplementary Figure S1). The N-terminus of S1PR1 is 40 amino acids long, in comparison to 39 amino acids in hMC4R. The template structure starts with residue Ser17, the final MC4R model contains amino acids Asn17-Pro321.

Structural modifications to generate a human MC4R homology model were performed with the software Sybyl X2.0 (Certara, Princeton, NJ, United States). For template preparation, the T4-lysozyme and the co-crystallized ligand were removed, loop lengths were adjusted by deletion of template amino acids (e.g., for ECL2) and gaps of missing residues in the loops of the template structure were closed manually (e.g., in ICL2 and ICL3). Moreover, missing residues between the N-terminal helical part and TM1 (between Ala39-Leu47) were added manually. Amino acids of the receptor-template were substituted with residues of human MC4R according to a sequence alignment between S1PR1 and hMC4R (Supplementary Figure S1), followed by conjugate gradient minimization of side chains until convergence to a gradient of 0.1 kcal/mol.Å with constraint backbone atoms of the transmembrane helices. This preliminary model was refined by MD simulations (300 K, 2 ns) of side chains and loops, followed by energy minimization of the entire model until convergence to a gradient of kcal/mol.Å. The quality and stability of the model (Figure ) were validated by checking the geometry with the program PROCHECK (Laskowski et al., 1993).

The final conformation of the hMC4R model with docked MCL0129, sodium binding and an unusual organization for the distal end of helices 7–8.

General Setup of the MD Simulations

Melanocortin receptor models were inserted into a 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) membrane bilayer according to the orientation in the OPM database (Lomize et al., 2006). Parameters for the simulations were generated using the CHARMM-GUI web interface (Jo et al., 2008; Wu et al., 2014; Lee et al., 2016) using the CHARMM36 force field (Huang and MacKerell, 2013) with CgenFF (Vanommeslaeghe et al., 2010) for ligands with partial charges generated from AM1-BCC calculations (Jakalian et al., 2002). Ligand molecules were added 7 Å above the receptor’s extracellular surface. The appropriate number of sodium and chloride ions was added to the systems to simulate a physiological salt concentration of 150 mM. Particle-mesh Ewald (PME) (Darden et al., 1993) was used to treat electrostatic interactions, using a cut-off distance of 10 Å. The resulting system was geometry-optimized and then equilibrated for 10 ns followed by a production run. All simulations used the TIP3P water model. All simulations were performed using GROMACS 5.0.4 (Pronk et al., 2013) and the Plumed plug-in 2.1 (Tribello et al., 2014) for the metadynamics simulations. The simulation systems compromised a box of 75 Å × 75 Å × 110 Å.

Metadynamics simulations were performed following a recently published protocol (Saleh et al., 2017b). We used a combination of the well-tempered variant (WT) (Laio and Parrinello, 2002; Barducci et al., 2008) of metadynamics and funnel metadynamics (FM) (Limongelli et al., 2013). A metadynamics history-dependent bias was applied along the component of the z-distance between the relatively immobile Cα of Trp6.48 deep in the ligand-binding region and the center of mass of positively charged nitrogen atoms of MCL0129. This distance was used as the single collective variable. The funnel restraint was then applied to the relative position on the xy-plane (restrained to 8 Å radius when fully unbound and allowing 13 Å radius if around the vestibule or deeper) to ensure better sampling for the relevant region of the free energy, because the ligand can otherwise move extensively in the extracellular solvent without affecting the free energy. Gaussian hills with an initial height of 0.48 kcal mol-1 were applied every 1 ps. The hill width was chosen to be 1 Å. The Gaussian functions were rescaled in the WT scheme using a bias factor of 20.

Representative structures were extracted along the simulation of initial docking for each 2 Å window, and used as starting coordinates for the multiple walker technique (Raiteri et al., 2006). This ensured faster convergence of the free-energy surface and enhanced the parallelization. Free energies were calculated using the PLUMED plug-in (Tribello et al., 2014), as described in our protocol (Saleh et al., 2017b).

Homology Modeling of the Inactive State hMC1R and hMC3R Structures

It is reported that the ligand MCL0129 is selective among the hMCR subtypes (Sabban and Serova, 2018) and we identified residues in the suggested hMC4R ligand binding site that are similar, identical, or specific compared to hMC1R, hMC2R, hMC3R or hMC5R (Supplementary Figure S1). The diverse amino acids may help explain ligand selectivity between the MCR subtypes; to test this hypothesis we also performed binding dynamics with MCL0129 at hMC1R and hMC3R. The hMC1R and hMC1R models were built based on the final hMC4R/MCL129 docking complex that resulted from the 5 μs unbiased simulation. For this purpose, the ligand was deleted from the hMC4R/MCL0129 complex. Missing residues [such as in hMC3R position Asp110 (ECL1) or Pro226-Ala227 (ICL3)] were inserted manually. Amino acids of the hMC4R model were substituted with corresponding residues of human hMC1R and hMC3R, respectively, according to the sequence alignment (Supplementary Figure S1), followed by conjugate-gradient minimization of side chains until converging at a termination gradient of 0.1 kcal/mol.Å with constraint backbone atoms of the transmembrane helices. These initial hMC1R and hMC3R models were refined by MD simulations of side chains and loops (300 K, 1 ns), followed by energy minimization of the entire model until convergence to a gradient of 0.05 kcal/mol.Å.

Results

In this study, spontaneous binding of the selective MC4R non-peptide antagonist MCL0129 (1-[(1S)-1-(4-fluorophenyl)-2-[4-[4-(2-methoxynaphthalen-1-yl)butyl] piperazin-1-yl]ethyl]-4-propan-2-ylpiperazine) was observed during 5 μs unbiased MD simulation (Chaki et al., 2003) (see section “Materials and Methods”). These simulations were analyzed to

elucidate the pathway and mechanism of MCL0129 binding to hMC4R,

detect binding intermediates and calculate binding affinities of MCL0129 using metadynamics simulations and to

explore the basic structural features of hMC4R in a stabilized antagonist-bound inactive state.

Structural Features of the hMC4R Homology Model

The initial hMC4R model is based on the crystal structure of a lipid GPCR, the sphingosine 1-phosphate receptor 1 (S1PR1, PDB entry 3V2Y) (Hanson et al., 2012), which shares high sequence similarity (55%) and identity (30%) with hMC4R in the transmembrane region (Supplementary Figure S1). Strikingly, both receptors lack two prominent structural features: (a) Pro5.50 in TM5 and (b) a disulfide bridge between TM3 and ECL2. The absence of Pro5.50 [superscripts are according to the Ballesteros and Weinstein numbering scheme (Ballesteros and Weinstein, 1995)] (Supplementary Figure S1), which causes a kink and a bulge in TM5 of most class A GPCRs (Worth et al., 2017) [as in known structures of adrenergic-receptors (Sabban and Serova, 2018) and discussed in detail by Sansuk et al. (2011)], leads to a straight helix conformation in S1PR1 (Troupiotis-Tsailaki et al., 2017).

The straight TM5 conformation is maintained in 5 μs unbiased MD simulations of the hMC4R model. On the cytoplasmic side, sodium spontaneously binds to a site defined by Ser581.46, Asp902.50, Ser2957.46, and Asn2947.45 (Figures , ). Binding of sodium to the highly conserved Asp2.50 has also been observed in the NMR structure of the MC4R TM2 (Yun et al., 2015) and in the crystal structures of the protease-activated receptors 1 and 2 (Zhang et al., 2012; Cheng et al., 2017), as well as in the high-resolution structure of the A2A adenosine receptor (Figure ) (Liu et al., 2012). Finally, a 5–6 Å outward movement of the distal end of TM7 relative to the template structure was observed, accompanied by a displacement of helix 8 (H8). H8 swivels inwards between TM1 and TM7 to be fixed by polar interactions, including a salt bridge between K3148.55 of H8 and Asp2987.49 of TM7 (Figure ).

Structural–functional specificities observed in the hMC4R model. In the resulting hMC4R homology model (A) formation of a stable sodium binding pocket constituted by hydrophilic amino acid residues Ser1.46 (Ser58, TM1), Asp2.50 (Asp90, TM2), Ser7.46 (Ser295, TM7), and Asn7.45 (Asn294, TM7) was observed (dotted yellow lines indicate potential hydrogen bonds), surrounded by further hydrophilic residues in TM1 (Asn62, Glu61) and in TM2 (Ser94). Sodium binding under participation of the highly conserved amino acid Asp2.50 was already reported for the determined structures of protease-activated receptor (PAR) 1 [PDB entry 3VW7 (Zhang et al., 2012)] and PAR2 [PDB entry 5NDD (Cheng et al., 2017)], or as shown in (B) for the high-resolution structure of the A2AR [PDB entry 4EIY (Liu et al., 2012)]. Potential water molecules additionally surrounding the Na+ ion in a cluster like arrangement. We further note that during our MD simulations (C) the C-terminus of TM7 moved (∼5.6 Å) toward the membrane (red arrow, MC4R white backbone) in comparison to the template structure of TM7 (blue backbone-ribbon of S1PR1). This process was accompanied by an inward displacement of H8 (red arrow). This motion resulted in potential polar interactions including salt bridges between positively charged residues in H8 and a negatively charged side chains in TM7 (Asp298) below to the sodium binding pocket (shown in A). In the template structure (C, blue backbone) H8 localization and the position of residue 7.49 are similar to the commonly observed features in determined class A GPCR structures.

Of note, the crystallization construct used for structure determination of the S1PR1 template is characterized by three stabilizing mutations located directly at the interface of the TM7-H8 junction, at residues Lys2506.29, Ser2516.30, and Leu2526.31. To test whether these substitutions cause the observed differences between the S1PR1-template structure and the final hMC4R model, 2 μs MD simulations were performed on the S1P1R crystal structure with the original wild-type amino acids at these three positions. We find that the S1PR1 wild-type model indeed features high structural flexibility within this region, with similar RMSD to that observed in the hMC4R model (>5 Å RMSD, Supplementary Figure S2).

Multi-Step Binding of the Antagonist MCL0129 to the hMC4R Model

In the 5 μs unbiased MD simulations described above, one MCL0129 molecule was placed in the extracellular solvent phase approximately 7 Å from the hMC4R surface. A funnel-like restraint centered on the receptor [as described in reference (Saleh et al., 2017b)] was set to focus sampling of the ligand close to the binding pocket. It is worth noting that this work reports the first use of a funnel-like restraint, which has so far been used exclusively for metadynamics sampling, in successful binding using unbiased MD simulation (Limongelli et al., 2013; Troussicot et al., 2015).

Spontaneous binding of MCL0129 was observed to the ligand-binding pocket after ≈3 μs (Figure ). The antagonist undergoes a four-step binding mechanism with three intermediate steps before settling at the hMC4R in its final position. Both the N-terminus (Nterm) and ECLs participate in recognition and spatial orientation of ligand binding. The highly flexible Nterm acts as a hook that traps the MCL0129 to guide it into the final binding mode, where the helical part of the Nterm caps the antagonist as a lid that prevents it leaving the pocket.

(A) Root-mean-square deviation of MCL0129 binding to hMC4R with a scheme for the structure of the non-peptidic antagonist MCL0129. (B) Free-energy profile for MCL0129 binding to hMC1R, hMC3R, and hMC4R with (C) representatives structures for the observed binding steps at MC4R (states 1–4).

In order to verify the final binding mode of MCL0129 obtained from 5 μs unbiased MD simulation, we used metadynamics simulation of 2 μs collective sampling. We followed the recently published funnel-metadynamics-based protocol (Saleh et al., 2017b) to characterize the global minimum for the hMC4R-MCL0129 complex. The metadynamics simulation was run independently from the unbiased MD simulation, using only the final pose of the MCL0129-hMC4R as starting geometry. We find that the final binding mode of MCL0129 obtained from unbiased MD simulation corresponds to the global minimum with ΔGbinding = -11.5 ± 0.3 kcal mol-1 (Figure ) and Ki ≈ 6.5 nM [experimental Ki = 7.99 nM (Pontillo et al., 2005)]. This binding mode (state 1) is defined by contacts of MCL0129 to amino acids from the Nterm, ECL1-3 and TM1-7 except for TM5 (Figures and Supplementary Figure S3). Specifically, hydrogen bonds are formed with Glu29Nterm, Thr1012.61, and Asp1263.25, in addition to aromatic interactions with Phe511.39, Phe2847.35, and Tyr2686.58 (Figure and Supplementary Table S1).

Binding mode of MCL0129. (A) Side view showing MCL0129 in magenta and potentially contacting residues as blue and cyan sticks. Residues identified by mutagenesis to interact specifically with agonists are indicated in green. Those side chains reported to interact with agonist and overlap with here suggested MCL0129 interactions are indicated cyan. (B) Top view for the binding mode of MCL0129-hMC4R.

Antagonist Versus Agonist Interactions at hMC4R

The binding mode of MCL0129 matches experimentally identified interactions of peptidic ligands, including the endogenous inverse agonist AgRP with MC4R (Yang et al., 2000; Nickolls et al., 2003; Tao and Segaloff, 2003; Chen et al., 2007; Tao, 2010). The acidic residues Asp1223.25 and Asp1263.29, which have been predicted to be directly involved in AgRP binding (Yang et al., 2000; Nickolls et al., 2003; Tao and Segaloff, 2003; Chen et al., 2007), form an ionic interaction with the positively-charged piperazine moiety of MCL0129 in a sandwich-like mode (one above and one below).

Comparison of the MCL0129 antagonist binding site with the proposed agonist-binding site (reviewed in Tao, 2010; Ericson et al., 2017a) provides an explanation for the different pharmacological effects. The spatial overlap of binding at the extracellular sides of TM1, TM2, TM3, and TM7 (Figure ) suggests that MCL0129 prevents access of the agonist to the binding pocket, and thus acts by direct competition. MCL0129, however, does not interact with the specific residues Ile1253.28, Leu1333.36, Trp2586.48, Phe2616.51, and Met2927.43 (Figure ) from the putative agonist-binding pocket (Yang et al., 2000; Nickolls et al., 2003; Tao and Segaloff, 2003; Chen et al., 2007; Tao, 2010). We assume that these additional contacts located deeper in the TM bundle define agonist- vs. antagonist-ligand binding properties.

MCL0129 Selectivity for MCR Subtypes

Amino acid variations of the MCL0129 binding pocket between all five melanocortin receptors were analyzed to evaluate subtype selectivity of MCL0129 (Pontillo et al., 2005) (Figure ). Briefly, three residues Ser25Nterm, Val461.34, and Tyr2686.58 are found exclusively in hMC4R, while residues Glu29Nterm and Pro481.36 are shared between hMC4R and MC2R or MC3R, respectively (Supplementary Figure S1). To test whether these differences contribute to subtype selectivity, we analyzed MCL0129 binding to hMC1R and hMC3R by an additional set of MD simulations. For this purpose, hMC1R and hMC3R subtype models were built based on the refined hMC4R structure obtained after 5 μs MD simulations. Next, these subtype models were simulated with the final docking pose of MCL0129 obtained for hMC4R for 2 μs of unbiased MD simulation each. The high RMSD values of 7 and 13 Å observed for the MCL0129 ligand with hMC1R and hMC3R, respectively, reveal high positional variability, indicating lower affinity to the MCL0129 (Supplementary Figure S2B). The subsequent metadynamics simulations highlight that MCL0129 indeed binds very weakly (with binding free energies less than -5 kcal mol-1 for hMC1R and hMC3R, respectively), in agreement with the experimental observation that neither subtype binds MCL0129 up to 10 μM (Pontillo et al., 2005).

Conservation of ligand binding site residues among the hMCR group. (A) Red: highly conserved in sequence; Orange: low sequence conservation (same residue in one further hMCR subtype or maintained biophysical properties like hydrophobic side chains); Blue: hMC4R specific amino acids (no sequence conservation); Residues shown as sticks form potential contacts between hMC4R and MCL0129; Residues shown as lines are side chains suggested to be involved in action of agonists, but not for MCL0129 binding. (B) Interaction diagram for the MCL0129-hMC4R complex color-coded based on the conservation of these interacting residues among the melanocortin receptors.

Discussion

This study was intended to gain insight into the structure and specific ligand-binding properties of MC4R, a high-priority drug target for treating obesity (Zhou and Cai, 2017; Ju et al., 2018). In the absence of an experimental structure, we used microsecond-scale MD simulations to study binding of the non-peptide antagonist MCL0129 to a homology model of human MC4R (hMC4R). To elucidated subtype selectivity, we evaluated the binding properties of MCL0129 to structural models of the homologous hMC1R and hMC3R receptors.

The initial hMC4R model is based on the crystal structure of the homologous sphingosine 1-phosphate receptor 1 (Hanson et al., 2012), which is characterized by the absence of a proline-induced kink in TM5 and a disulfide bridge between TM3 and ECL2. The absence of Pro5.50 causes a regular α-helical conformation different from other GPCRs that have a kink and a bulge in TM5 (Pro5.50 is approximately 80% conserved) (Sansuk et al., 2011). The 5 μs unbiased MD simulations of the hMC4R suggest that this straight TM5 conformation is maintained. The regular α-helical conformation of TM5 particularly defines the residues that form the putative ligand-binding site suggested here (Figure ). A different set and spatial localization of residues would become accessible for interaction with the ligand in a kinked and slightly rotated TM5.

Moreover, we identified a sodium-binding site defined by amino acid residues in TM 1, 2 and 7, which is known to impact ligand binding (allosteric effects) and protein stability (Liu et al., 2012). Key players in this motif are amino acids Ser581.46, Asp902.50, Ser2957.46, and Asn2947.45 (Figure ). This observation is in general agreement with NMR data on MC4R and crystal structures of other GPCRs (Zhang et al., 2012; Yun et al., 2015; Cheng et al., 2017), with the difference that the putative hMC4R sodium-binding site is located between TM2 and TM7, close to TM1, while in the other GPCRs it is more located towards TM3.

Finally, we find that in our simulations, a characteristic outward movement of TM7 by a 5–6 Å in the distal end relative to the template or other GPCR crystal structures. This novel structural feature is accompanied by a displacement of H8 in-between TM1 and TM7 (Figure ), which indicates high flexibility of TM7 and H8. In agreement, in the homologous template structure this part was thermo-stabilized by site-directed mutations (Supplementary Figure S2).

In accordance with previous computational and experimental work carried out for other receptors (Dror et al., 2011; Bock et al., 2012; Kruse et al., 2012; DeVree et al., 2016; Milanos et al., 2016; Saleh et al., 2016, 2017a,c), our study suggests a stepwise binding mechanism for the selective non-peptide antagonist MCL0129 to hMC4R. The Nterm plays a key role in ligand binding and recognition, trapping MCL0129 into its final binding mode.

This binding side of MCL0129 is defined by contacts with ECL1-3 and TM1-4 and TM6-7, for which we assume a straight helical conformation in the absence of Pro5.50. MCL0129 is stabilized by hydrogen bonds with Glu29Nterm, Thr1012.61 and Asp1263.25, in addition to aromatic interactions with Phe511.39, Phe2847.35, and Tyr2686.58 (Figures , and Supplementary Table S1).

Partial overlap of the antagonist binding site with the proposed agonist binding site suggests that MCL0129 directly blocks the agonist binding (Figure ). In addition, several of the contacts proposed to be essential for agonist binding are missing in the hMC4R/antagonist complex, thus explaining the antagonistic properties at hMC4R. We identified several hMC4R residues as potential selectivity determinants between hMC4R and other hMCRs (three are hMC4R-specific and two are highly variant), located in both the 7-TM bundle and the Nterm. Such differences between MCR subtypes may affect future efforts to design selective ligands for the MCR’s.

Finally, the ability of MD simulations combined with free-energy calculations to track and identify intermediate steps of ligand binding has increasingly contributed to our understanding of GPCR-ligand interactions and helped to identify potential binding sites and strategies for SBDD, which remains challenging for many GPCRs that are still not therapeutically accessible. Our MC4R-related insights obtained from computational simulations may pave the way for rational ligand design in the field of medicinal chemistry targeting the melanocortin systems.

Author Contributions

NS: project idea, modeling and simulation studies, data analysis, wrote the manuscript, and final approval. GK: project idea, modeling studies, data analysis, wrote the manuscript, and final approval. NH and TC: data analysis, wrote the manuscript, and final approval. PH and PS project idea and coordination, data analysis, wrote the manuscript, and final approval.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.