Received Date : 06-Apr-2014

Accepted Article

Revised Date : 29-May-2014 Accepted Date : 13-Jun-2014 Article type

: Research Article

A Conformational Analysis Study on the Melanocortin 4 receptor Using Multiple Molecular Dynamics Simulations

Mohsen Shahlaei1*, Atefeh Mousavi2 1

Novel Drug Delivery Research Center, School of Pharmacy, Kermanshah University of

Medical Sciences, Kermanshah, Iran 2

Student Research Committee, Kermanshah University of Medical Sciences, Kermanshah,

Iran

*Corresponding author: Mohsen Shahlaei, PhD Department of Medicinal Chemistry, Faculty of Pharmacy, Kermanshah University of Medical Sciences, 67346-67149, Kermanshah, Iran.Tel: +98-831-4276489; Fax: +98-8314276493; E-mail: [email protected] AND [email protected]

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process which may lead to differences between this version and the Version of Record. Please cite this article as an 'Accepted Article', doi: 10.1111/cbdd.12495 This article is protected by copyright. All rights reserved.

Abstract

Accepted Article

Taking into account the uncertainties involved in 3D model of biomolecule developed by homology modeling (HM), it is important to opportunely validate the initial structure before employing for different purposes such as drug design. Extended simulation times and the necessity of correct representation of interactions within the protein and the nearby molecules impose significant limitations on molecular dynamics (MD) based refinement of structures developed by HM. Consequently, there is a pressing requirement for more efficient methods for HM and subsequent validation of developed structure. Multiple MD simulation runs are well suited for producing ensembles of structures. In this context, a computational investigation was presented to study the structure of Melanocortin 4 receptor (MC4R) using molecular dynamics (MD) simulations in explicit phospholipids bilayer. Several MD runs with different initial velocities were employed to sample conformations in the neighborhood of the native structure of receptor, collecting trajectories spanning 0.21 ms.

The coherence between the results, different structural analysis and the convergence of parameters derived by principal component analysis (PCA) shows that an accurate description of the MC4R conformational space around the native state was achieved by multiple MD trajectories.

Keyword: Melanocortin 4 receptor, Molecular Modeling, Multiple Molecular Dynamics Simulation, Homology Modeling Uncertainties

This article is protected by copyright. All rights reserved.

Introduction

Accepted Article

G-protein coupled receptors (GPCRs), a well known superfamily of transmembrane proteins, have an important role in signal transduction. They are activated by a diverse set of signals, involved in signal transduction pathway, including small molecules and peptides (1). Around 50% of all available pharmaceutical agents on the market act on GPCRs (2). Despite their pharmaceutical application, only a few three-dimensional (3D) structures of GPCRs are available up to now (3-8). However, there are a number of other ways to overcome the lack in structural information on other GPCRs that one of the most important is HM technique (9, 10).

GPCRs superfamily consists of about 800 receptors which can be subdivided into distinct subcategories regarding the similarities in the sequence homology and evolutionary origin (11). GPCR subfamilies are defined by the substrates e.g. histamine receptors.

The human melanocortin 4 receptor (MC4R), one of the most important members of melanocortin receptor family, is a 332-amino acid GPCR (12). This receptor mediate signal transduction induced by melanocortin peptide via G-protein coupled signaling system. Human MC4R is found to express only in the brain and be involved in the Leptin pathway; it contributes to the controlling of human feeding behavior and energy homeostasis (13).

The knowledge of the 3D structure of a protein provides important insights into the molecular basis of protein function. In particular, it plays an important role in drug design and discovery, such as design of specific ligands. Although great progresses have been achieved in the field of structural biology by experimental techniques such as X-ray crystallography and nuclear magnetic resonance spectroscopy, it is still a time-consuming and costly task

This article is protected by copyright. All rights reserved.

without guaranteed success. Currently, over 100,000 experimental protein structures are

Accepted Article

deposited in the Protein Data Bank (PDB) (14). Nevertheless, the number of structurally determined proteins is low compared to the number of known protein sequences. Thus, no structural data are available for the large number of proteins. Therefore, it is an obvious demand to bridge this ‘structure knowledge gap’. This is the main reason for the increased interest in computational methods for protein structure prediction in the recent years. Up to now, it is not possible to reliably predict protein structures from sequence employing only ab initio approaches. Nevertheless, among all current computational methods and techniques, HM is one of the most common techniques to obtain an accurate 3D model of a protein from its amino acid sequence (15, 16) exploiting the similarity between the target protein and a (or more) template with known experimental structure.

A primary objective of HM is to provide atomistic details for the protein structures where there is no enough experimental information (17). HM is a knowledge based method relying upon the crystallographic structure of related receptor and experimental information.

However, homology models are imprecise by definition in particular in terms of predicting side-chain orientation or disordered protein regions as loops. Indeed, considering the uncertainties involved in HM, it is very important to opportunely validate the initial 3D model before employing its structure for drug-design purposes (18). Said another way, the extent to which homology models can be used with confidence is unclear. This is in part due to alignment errors, but primarily due to the lack of effective methods that can be used to refine the structural models obtained (19).

This article is protected by copyright. All rights reserved.

Present attempts to refine homology models are normally based either on energy

Accepted Article

minimization, limited conformational sampling using molecular dynamics (MD). Still, these approaches are generally considered to be ineffective (20-22). Failure of refinement of homology models based on simplified representations and/or limited sampling is not surprising as proteins are densely packed rendering searching their conformational spaces rather difficult (23). On the other hand, extended simulation times and the requirement of accurate representation of interactions within the protein and the surrounding environment impose significant limitations on MD-based refinement of homology models (23-26). Therefore, there is a pressing need for more efficient approaches for HM and subsequent structure validation. Multiple MD simulation runs are well suited for producing ensembles of structures at a fixed temperature, and, actually, they have been successfully employed to study the conformational landscape near a protein native state (27-31).

The shortcomings of HM and MD simulations combined with the general difficulty in identifying binding sites in GPCRs prompted us to propose a ligand-guided computational approach to evaluate the validity of homology models and to explore the suitability of their putative surface cavities as potential binding sites.

In this context, we recently applied HM to the study of different chemokine receptors (3234). These models include not only transmembrane α-helices but also all hydrophilic loops. Also, recently, we are reported a homology model of MC4R (35) . Here a computational study on the MC4R structure and the binding modes of MC4R agonists was proposed using developed 3D structure of MC4R reported previously and MD simulations. To the best of our

This article is protected by copyright. All rights reserved.

knowledge, there have been few studies on the interaction of drugs and MC4R, and still no

Accepted Article

theoretical research has been related to the dynamics of the MC4R system. Additionally, it seems to be possible to obtain more accurate model using molecular dynamics (MD) simulations. Our simulations have been carried out in a membrane environment, modeling the receptor embedded in a phospholipid bilayer solvated by water, since it is suggested to be the most adequate environment to simulate membrane receptors (36-39).

Although changes in diet and exercise underlie the current global increase in the prevalence of obesity, there is considerable evidence of a substantial genetic contribution to the regulation of body weight (13). Over the past few years, research efforts in drug discovery groups have focused on identifying novel MC4R agonists for the potential treatment of obesity and sexual dysfunction. Numerous reports on the identification of structurally related small-molecule MC4R agonists have emerged in the literature as well as structurally distinct MC4R agonists. Among these agonists, we focused on a potent and selective, small molecule agonist of the MC4R that are efficacious in animal models for both food intake and erectile dysfunction (40).

2. Materials and methods

2.1 MD simulation The starting structure for the simulations was taken from the HM (35) ( Fig. 1). MD simulations were carried out in two parts. In the first part, the minimum energy conformation of MC4R was used as starting structure for MD simulation. In this part components of simulation system included protein, lipid, and water molecules. In the second

This article is protected by copyright. All rights reserved.

phase, output conformation resulted from the first phase of MD as well as a joined ligand

Accepted Article

were used in MD process.

All of the MD simulations were carried out by the GROMACS 4.5.3 package (41) and in an explicit phospholipid bilayer. The Gromos G53A6 force field (42) was applied for the MD simulation of MC4R, while the parameters developed by Berger et al. were used to the 1palmitoyl-2-oleoyl-phosphatidylcholine (POPC) lipids (43). A pre-equilibrated POPC bilayer with 128 lipid molecules was downloaded from Tieleman's group electronic page: (http://moose.bio.ucalgary.ca/index.php?page=Structures_and_Topologies). The MC4R was embedded in the pre-equilibrated POPC lipid bilayer using inflategro perl script from Tieleman's group (44) The protein was inserted at the center of the POPC bilayer with its long axis perpendicular to the membrane-water interface. The α-helical domain of the receptor was placed at the same level as the lipid bilayer, and overlapping lipid molecules were discarded.

All systems were equilibrated under an isothermal–isobaric (NPT) ensemble for 1 ns. The LINCS method was applied to restrain all protein heavy atoms in any direction, and on the phosphorus atoms of the lipid head groups in the vertical (z) direction allowing a 2 fs integration step (45). In the next step, water was added to fill the gaps above and below the phospholipid leaflets and box edges and to hydrate the lipid head groups.

Water molecules were represented applying a simple point charge (SPC) model. Solvent (i.e., water), ions, lipid, protein, and ligand were coupled separately to a temperature bath. To balance the net charge associated with the presence of POPC, 7 water molecules were randomly replaced by Cl-. This number of cholorine ions must stay constant to preserve the

This article is protected by copyright. All rights reserved.

electroneutrality of the system of interest. Final protein-POPC membrane system was placed

Accepted Article

in a box with the dimensions of 103.6 × 99.1 × 117.8 (all in Å) with a total of 95880 atoms.

After the starting system has been constructed, the system was energy minimized before the MD simulation using 200 steps of the steepest descent method in order to relax any steric conflicts generated during the setup. After stabilization of temperature, an NPT ensemble was performed. In this phase a constant pressure of 1 bar was employed with a coupling constant of 5.0 ps and time duration. Water, lipids and protein were coupled separately to a temperature bath at 310°K with a relaxation time of 0.2 ps using a Berendsen thermostat (46).

The particle mesh Ewald (PME) method was used (47). This treatment of electrostatics has been reported to provide an accurate representation of lipid properties and is also commonly used in simulations of proteins (48). A 12 Å cut off for long-range electrostatic interactions and the Lincs algorithm for covalent bond constraints were applied (45) with the neighbor list updated every five simulation steps. Also, A 10 Å cutoff was used for the Van der Waals interactions. The molecular topology for ligand was generated by PRODRG2 server at http://davapc1.bioch.dundee.ac.uk/prodrg/ (49). All the molecular images were produced using VMD (50). The trajectories were analyzed using the standard tools included in the GROMACS distribution. The Cα root-mean-square-deviation (rmsd), which is a crucial parameter to evaluate trajectory stability, was computed using the starting structures of the MD simulations as a reference. Time interval for rmsd calculation was 5 ps. To improve the conformational sampling, twenty one independent 10 ns simulations (replicates) were performed, initializing the MD runs with different initial atomic velocities or using different initial structures. For each simulation (replicate), the equilibrated portions were joined

This article is protected by copyright. All rights reserved.

together to obtain concatenated trajectory of different lengths. Details of trajectories of

Accepted Article

replicates and concatenated trajectories are reported in Table 1.

For each protein, the equilibrated portions of each replicate were joined together to obtain concatenated trajectories, which are representative of different directions of sampling around the starting structure. The secondary structures were calculated using DSSP (51).

2.2 Principal component analysis Principal component analysis (PCA) is a method that can be used to compress a pool of data into principal components (PCs). PCA is a linear transformation that changes a set of (possibly) correlated variables into a set of orthogonal variables (PCs). In new array, PCs can be sorted on the basis of corresponding eigenvalues. The first PC is the PC with the largest eigenvalue and each succeeding PC was emerged in order of decreasing eigenvalues. The strongest correlations in the dataset tend to be captured by the first PCs. PCA can be applied for decreasing dimensionality by keeping only a few of the first PCs, hopefully without much loss of information. In a typical protein simulation, PCA of an ensemble of conformations can be calculated and diagonalized using the covariance matrix of the Cartesian coordinates of a representative subset of atoms in the molercule of interest, typically the Cα atoms. The covariance between two coordinates xi and xj is defined as:

This article is protected by copyright. All rights reserved.

Where the bar denotes an ensemble average. As MD simulations tend to create huge

Accepted Article

quantities of data, PCA is an influential mathematical tool employed to detect correlations in MD trajectories (52, 53). This technique permits one to analyze the important characteristics of trajectories and filter out random, insignificant fluctuations. Said another way, the PCA of the protein positional fluctuations along the MD trajectories is performed to evidence the presence of correlated residue motions. The output of PCA is two matrices: A matrix contains the eigenvectors, and a matrix of the corresponding eigenvalues. These eigenvectors are sorted as the information content (variance) that they explain decreases from the first PC to the last one. In other words, the eigenvector with the highest eigenvalue is considered the first PC; the eigenvector with the second highest eigenvalue is considered the second PC and so on. The eigenvectors stand for the direction of motion, and the eigenvalues correspond to the amount of motion along the eigenvectors. It has been proposed that the majority of the motion for proteins can be accounted for in the first some PCs. When the ensemble of protein conformations obtains from a MDS procedure, each of the eigenvectors resulting from the diagonalization, or PCs, explains a collective mode of motion (i.e., involving all degrees of freedom) that is not linearly correlated with any other in the system. The extent of this motion, or fluctuation, is given by the corresponding eigenvalue (54). Therefore, the dynamics of a protein can be studied by projecting its atomic motion during a MD simulation onto its first some PCs. PCA is an efficient tool for monitoring protein dynamics because the observed motion is unconstrained and represents the atomic fluctuations of the protein.

To define the dimensionality of the essential subspace, the fraction of total motion described by the reduced subspace was considered and computed as the sum of the eigenvalues relative to the included eigenvectors, describing the amount of variance retained by the reduced representation of the total space. This article is protected by copyright. All rights reserved.

The root-mean-square inner product (RMSIP) (55, 56) was computed to measure the degree

Accepted Article

of overlap between the conformational subspaces of the protein explored by the simulation at replicates. The RMSIP is defined as follows:

Where ηia and ηjb are the ith and jth eigenvectors from two different replicates to be compared, a and b, respectively.

In order to estimate how much of the major essential modes of dynamics obtained from the covariance matrix resemble random diffusion, its cosine content (57, 58). ci of the eigenvector i was calculated as follows:

Where pi(t) is the amplitude of the motion along eigenvector i at time t. The cosine content

can take values between 0 (no cosine) and 1 (a perfect cosine), and provides an indicator of the extent of sampling.

3. Result and Discussion Clarification of ligand binding mechanisms is the essential step to introduce more selective and potent lead compounds for a given target protein. To do this, we need to construct a 3D model of protein and then let it feel its natural environment. MD simulation is one of the best methods for such refinement.

This article is protected by copyright. All rights reserved.

It is usually accepted that the most realistic and helpful approach for performing MD

Accepted Article

simulations of GPCRs is the use of a phospholipid bilayer solvated by water molecules. Nowadays, a large number of phospholipid bilayer models are introduced applying MD simulations. However, most of them include no more than 128 monomer of phospholipid. This seems to be enough for describing physicochemical characteristics of the phospholipid bilayer; however, it is insufficient for simulating a GPCR in the phospholipid bilayer environment. For this reason, a model of the POPC bilayer was used consist of 414 monomer of phospholipid. Such an amount of phospholipid monomers seems to be enough for inserting the MC4R receptor model into the bilayer.

Placing of the alpha helices of MC4R into the lipid bilayer core was carried out in such a way that alpha helices were perpendicular to the membrane plane and protein-overlapping lipids were removed.

In order to impair the incorrect interactions imposed to the 7 TMs during HM, a relaxation procedure was applied to the POPC–MC4R system. As a result, constraints between TM helices were reduced in sequential steps to avoid unwanted structural drifts during production phase of the MD. The POPC–MC4R system remained firm after the relaxation as very slight drift in energy, temperature, or lipids density was monitored during the MD (data not shown).

As discussed above, several independent simulations (replicates) for MC4R were performed in order to increase the conformational sampling, using an iterative procedure as explained in Methods section.

This article is protected by copyright. All rights reserved.

3.1 Essential dynamics and sampling

Accepted Article

One of the main problems pointed out by many researchers (59-63) which appears when trajectories are analyzed, is the convergence of conformational sampling. Since MD simulations are necessarily limited in time, conformations visited are only a subset of all the possible conformations that the protein can assume. In order to correlate the MD data with protein properties, one should ensure a sufficiently high sampling efficiency. When the sampling of MD trajectories is insufficient, protein motions along the PCs appear indistinguishable from the dynamics of random diffusion, not allowing an accurate description of the system of interest (64). Said another way, for short MD simulation procedures insufficient for convergence of sampling, the first several PCs can have the shape of a cosine function caused by random motion of the protein chain without potential barriers, as characteristic of Brownian motion (60, 65). The following analyses of the resulted trajectories have been performed discarding the portions of the replicates required by the system to reach stable values of Cα rmsd and of gyration radius, in order to guarantee that calculated parameters reflect the intrinsic characteristics of the system of interest.

Most of the variance in conformational fluctuations of protein structure is generally related to the first PCs, which, in the case of MC4R simulations, account for a major fraction of the total fluctuations in the concatenated trajectories, but explain a minimal part of the total fluctuations if a single replicate is considered. A reduction in the number of the significant PCs required to describe a significant portion of the total variance is indicative of a wider sampling of the conformational space in the concatenated trajectories.

This article is protected by copyright. All rights reserved.

Generally, in order to sample conformational space reliably, a wide region of the

Accepted Article

conformational space should be sampled and, at the same time, a partial overlap between different trajectories should be achieved.

In order to gain insights into the conformations visited by the system, a cluster analysis on the rmsd matrices and also an essential dynamics (ED) analysis were performed. Conformations from each independent replicate were clustered separately (Figure 1), suggesting that different initial velocities caused the trajectories to sample different conformational basins.

Cα rmsd values for average structures of twenty one replicates fall in the ~0.34–0.06 nm range (Figure 2), when cluster average structures are compared, indicating that the overall 3D structure is well conserved in the different replicates.

As discussed above, when the cosine content of the first few PCs is close to 1, the largestscale motions in the protein dynamics cannot be distinguished from the random diffusion in protein. Said another way, the analysis of cosine content can exclude the interpretation of random diffusion of atoms as a correlated motion. Therefore, the quality of ED analysis can be examined by checking the cosine content of the first principal components. Thus, the average cosine content of first 20 PCs for the concatenated trajectories was further analyzed as a function of the trajectory duration, allowing the evaluation of convergence of the conformational sampling. (Fig. 3). Cosine content for the first 20 PCs of concatenated trajectories with high time of simulation is very low which signifies that the simulations do not represent the intrinsic dynamical behavior of the receptor. It is to be mentioned that if the

This article is protected by copyright. All rights reserved.

sampling of MD trajectories is insufficient, protein motions along the PCs appear identical

Accepted Article

from the dynamics of random diffusion and then the cosine content of the numerical value of cosine content for first few PCs is close to 1, a perfect cosine.

As a general rule in this figure, with increases in time of simulation, the value of average cosine content for first 20 PCs is decreased, confirming that using a large number of MD simulations with different initial velocities is an effective approach to obtain a reliable conformational sampling.

In the figure 4 a graph of required eigenvector numbers to explain 70% of the total motion is represented for each simulation. As can be seen, PCA on single replicate of MC4R shows that, in general, about 5 eigenvectors are required to describe more than 70% of the total variance (Fig4. A), whereas, in the concatenated trajectory, this value is increased as duration time of simulation is increased (Fig. 4 B).

Generally, the first three PCs describe a consistent part of the protein dynamics, and they could be used as a reference subspace to analyze protein motion (30, 66). The projection of frames of concatenated simulation into the subspace defined by three first PCs was reported in Figure 1, confirming and clarifying results obtained from cluster analysis. In particular, a wide sampling of the conformational space with re-sampling of similar conformations was achieved in our simulations, showing that the essential subspace is well explored when concatenated trajectories are considered.

This article is protected by copyright. All rights reserved.

Passing all the quality assurance tests by the developed model, proposes that this model is a

Accepted Article

good estimation of MC4R structure, and can be used to describe various protein–ligand interactions and to study the relationship between the structure and function.

3.2 Two dimensional projections of the conformation ensemble Because of the complicated interaction of many factors, complex statistical analysis of conformational ensembles methods, such as PCA, appears to be the most suitable computational means to describe the flexibility.

PCA was carried out on the final concatenated ensemble (REP1_21), with the Cα atoms of the residues fit to starting structure. Resulting eigenvectors were ordered by descending eigenvalues, which represent the variance of the motion along the PCs. In Fig. 5, the first twelve largest eigenvalues are shown for each concatenated run reported in Table 1. As can be seen, more or less, for all concatenated trajectories the three first PCs can be recognized with prominent eigenvalues of comparable magnitude, which is s higher than the magnitude of other eigenvalues. These PCs with the largest eigenvalues represent correlated motions of the MC4R with the most significant standard deviations of the motion along the corresponding orthogonal directions.

A more comprehensive representation of the conformational motions of REP1_21 in receptor is provided by projections of the ensemble of conformations onto the planes defined by the most important PCs (Fig. 6). In figure 6, the spatial distributions of occupancies of the various conformational states are reported over the planes spanned by the first and second and by the first and third principal components. As can be seen in figure 6A, both distributions are smooth and without evidence of considerable clustering. Similar behavior

This article is protected by copyright. All rights reserved.

can be seen using Fig.6B, which show no significant clustering of conformational occupancy.

Accepted Article

The widths of the distributions shown in Fig. 6 are not directly dependent on the number of residues. The size of the histogram has been selected in proportion to the width of the distribution; so the maxima observed for MC4R not artifacts of the sampling but represent regions of increased occupancy. Since the occupancy of a conformational state is inversely proportional to the free energy of that state, the lack of clustering (67-69) that the global minima are broad or that there are generally no significant local minima.

the average structure of REP1_21 was extracted from concatenated trajectory and was used for further analysis.

we decided to perform a second MD simulation on the ligand–protein complex for further investigation the effects of ligand binding on the conformational behavior of protein.

Hence, a previously reported agonist of MC4R was docked in active site of protein (scheme 1)(40). Ligand–MC4R complex was selected as a representative for MD simulation. This complex was used as the starting conformation for an additional MD run. This complex was placed into the phospholipid bilayer as described in the computational methods section.

The aim of the second MD simulation was to get more precise ligand–receptor models in a

state close to the natural conditions and to explore the binding modes of the ligands further.

The stability of MC4R was investigated by a 15 ns MD simulation and the following Cα RMSD calculation. Results show that MC4R becomes equilibrated at 5 ns and afterwards

This article is protected by copyright. All rights reserved.

(figure 6). Then, from the curve of Figure 7, it can be seen that MC4R was stabilised by

Accepted Article

POPC and the RMSD value of protein is about 0.4 nm.

From Figure 7, it can be seen that the ligand seems stabilised judging by its RMSD value. The binding mode of the ligand established after MD simulation (Figure 8).

The analysis of the binding mode of the ligand obtained after MD simulation (Figure 8) suggests that the phenyl ring of the His283 is involved in π-π interactions with fluorobenzene ring of ligand.

The analysis of the results of molecular dynamics simulations of protein ligand system confirmed the existence of a suitable binding site located inside the transmembrane domain. Also, amino acid residues located in the TM2 and TM3 seem to be essential for binding with ligand.

It was found that some lipophilic residues, in particular Thr178, Tyr162 located in TM3, and Leu265 and Ser282 located in ECL3, could be essential for ligand recognition through lipophilic interactions.

To determine the protein residues that are more affect by MC4R interaction, the rms fluctuations (rmsf) of the Cα atoms of the MC4R were calculated and plotted as a function of the residue number for free and bound MC4R. As is shown in Figure 9, for bound MC4R the residues located in the hydrophilic loops have the greatest rmsf values, while the rmsf of the transmembrane region is less than 0.15 nm. The greatest rmsf values were obtained for

This article is protected by copyright. All rights reserved.

residues located in the longest hydrophobic loop (I3); in particular, a maximum value of

Accepted Article

0.65nm was obtained for Pro230 of the MC4R-ligand complex.

Importantly, simulation of MC4R in the presence and the absence of the ligand led to the different final receptor conformation and orientation. There are some differences in the backbone dihedral angle values, and these differences considerably affect the position of some residues in the Ramachandran Plot (data not shown).

3.5 Interaction energy analysis To obtain a detailed description of the effects of individual residues on binding affinity, a perresidue decomposition of the total energy was carried out to evaluate the energetic influences of critical residues on the binding on two system, protein and ligand before MD and after MD. The residues were separately studied, and the interaction, Van der Waals, and electrostatic energies of these residues and investigated ligand were determined. These include Thr178, Ala176, Asn 285, His283, Asn274, Tyr276, Cys277, Asn74 (Table 2). Energy of interaction between individual residues and ligands was calculated with the Calculate Interaction Energy protocol encoded in Discovery Studio. The complexes were typed with CHARMM force field and the dielectric model was assigned to Implicit DistanceDependent Dielectrics. All other parameters used were kept at their default settings. The residue-based decomposition of interaction energies in two systems identified several critical residues of MC4R. Table 2 lists the calculated energies contributions of these key residues of interest for two systems. As can be seen, major favorable energy contributions originate predominantly from Thr178, His264, Val179, Ile182, Leu288, and Tyr276 before MD and Thr162, Ser282, Asn274 and Ala176 after MD.

This article is protected by copyright. All rights reserved.

Thus, these critical residues may provide guidance for the rational design to discover more

Accepted Article

potent MC4R agonists.

Conclusion In the absence of experimental structures, homology models play an important role in the structural biology and structure-based drug discovery process. Despite clear and rapid achievements in crystallization and NMR methods, it is reasonable to assume that experimental structures for most of the therapeutically relevant targets will not be available in the near future. Thus, in silico homology modeling provides a viable cost-effective alternative to generate reasonably accurate models for drug discovery. To date, homology modeling has been successfully used to identify hits, to suggest accurate binding modes and ligand-receptor interactions, aid in mutagenesis experiments, rationalize SAR data and guide the optimization toward more potent ligands. However problems associated with template identification, accurate sequence alignment and refinement methods hinder a wider use of in silico generated models in the drug discovery process. Thus, there is a need to develop reliable ways to reduce modeling errors, specifically at pharmacologically relevant sites, and when target and template proteins share lower sequence similarities.

The MD simulations results presented in this study provide the first detailed work of the Melanocortin 4 receptor structure inserted into the phospholipid bilayer. A molecular model of the Melanocortin 4 receptor containing seven transmembrane helices connected by three extracellular and three intracellular hydrophilic loops was developed and described. This structure provides insights into the binding mode of MC4R agonists.

This article is protected by copyright. All rights reserved.

Considering the uncertainties involved in structure of MC4R developed by homology

Accepted Article

modeling (HM), the initial 3D model of biomolecule before employing its structure for different purposes and was validated using multiple MD simulation. Multiple MD simulation runs are well suited for producing ensembles of structures. In this context, several MD runs with different initial velocities were employed to sample conformations in the neighborhood of the native structure of receptor, collecting trajectories spanning 0.21 ms.

The coherence between the results different structural analysis and the convergence of parameters derived by principal component analysis (PCA) shows that an accurate description of the MC4R conformational space around the native state was achieved by multiple MD trajectories.

A focus on the dynamic and energetic characterization of the MC4R was carried out specialized in the recognition of the ligand. A MD simulation was performed on the MC4R and the ligand complex to explore effects of the presence of ligand on the protein binding pocket and types of interactions. The conformational changes of the MC4R occurring during MD simulations were studied, and the stable binding modes of the ligand were determined. According to the models presented, the involvement of the Thr 178, Ala176 located in TM IV, Asn 285 located in TM VIII, His 283, Asn 274, Tyr 276, Cys 277 and Asn74 and some other residues in ligand recognition was determined. The obtained results allowed us to explore the important features of binding modes of this ligand and to determine the amino acid residues involved in the recognition of the MC4R.

This article is protected by copyright. All rights reserved.

In conclusion, we here provide the first detailed report of the Melanocortin 4 receptor

Accepted Article

structure inserted into the phospholipid bilayer and its complex with one of the most important agonists know so far.

The results presented in this study provide guidelines for experimental site-directed mutagenesis and if confirmed, they may be useful in the designing of new selective MC4R ligands with agonistic properties.

Since the above work is an in silico study, the developed model can be useful to develop new drug candidate against obesity. The above work aims to serve all those researcher and patient who are currently experiencing this disorder.

Acknowledgment The authors gratefully acknowledge the constructive suggestions by the Dr. Elena Papleo. The authors also gratefully acknowledge Vice Chancellor for Research and Technology, Kermanshah University of Medical Sciences for financial support. This article resulted from the Pharm. D thesis of Atefeh Mousavi, major of Pharmacy, Kermanshah University of Medical Sciences, and Kermanshah, Iran

References 1.

Ballesteros J, Kitanovic S, Guarnieri F, Davies P, Fromme BJ, Konvicka K, et al.

(1998) Functional microdomains in G-protein-coupled receptors the conserved arginine-cage motif in the gonadotropin-releasing hormone receptor. Journal of Biological Chemistry;273: 10445-53.

This article is protected by copyright. All rights reserved.

Drews J (2000) Drug discovery: a historical perspective. Science;287: 1960-4.

3.

Okada T, Sugihara M, Bondar A-N, Elstner M, Entel P, Buss V (2004) The retinal

Accepted Article

2.

conformation and its environment in rhodopsin in light of a new 2.2 Å crystal structure. Journal of molecular biology;342: 571-83. 4.

Warne T, Serrano-Vega MJ, Baker JG, Moukhametzianov R, Edwards PC, Henderson

R, et al. (2008) Structure of a &bgr; 1-adrenergic G-protein-coupled receptor. Nature;454: 486-91. 5.

Cherezov V, Rosenbaum DM, Hanson MA, Rasmussen SG, Thian FS, Kobilka TS, et

al. (2007) High-resolution crystal structure of an engineered human β2-adrenergic G protein– coupled receptor. Science;318: 1258-65. 6.

Jaakola V-P, Griffith MT, Hanson MA, Cherezov V, Chien EY, Lane JR, et al. (2008)

The 2.6 angstrom crystal structure of a human A2A adenosine receptor bound to an antagonist. Science;322: 1211-7. 7.

Wu B, Chien EY, Mol CD, Fenalti G, Liu W, Katritch V, et al. (2010) Structures of

the CXCR4 chemokine GPCR with small-molecule and cyclic peptide antagonists. Science;330: 1066-71. 8.

Chien EY, Liu W, Zhao Q, Katritch V, Han GW, Hanson MA, et al. (2010) Structure

of the human dopamine D3 receptor in complex with a D2/D3 selective antagonist. Science;330: 1091-5. 9.

Wolf S, Böckmann M, Höweler U, Schlitter J, Gerwert K (2008) Simulations of a G

protein-coupled receptor homology model predict dynamic features and a ligand binding site. FEBS letters;582: 3335-42. 10.

Michino M, Abola E, Participants GA, Brooks III CL, Dixon JS, Moult J, et al. (2009)

Community-wide assessment of GPCR structure modeling and docking understanding. Nature reviews Drug discovery;8: 455.

This article is protected by copyright. All rights reserved.

11.

Marshall FH, Foord SM (2010) Heterodimerization of the GABA< sub> B

Accepted Article

Receptor—Implications for GPCR Signaling and Drug Discovery. Advance pharma;58: 6391. 12.

Fan W, Boston BA, Kesterson RA, Hruby VJ, Cone RD (1997) Role of

melanocortinergic neurons in feeding and the agouti obesity syndrome. 13.

Farooqi IS, O'Rahilly S (2008) Mutations in ligands and receptors of the leptin–

melanocortin pathway that lead to obesity. Nature Clinical Practice Endocrinology & Metabolism;4: 569-77. 14.

Westbrook J, Feng Z, Chen L, Yang H, Berman HM (2003) The protein data bank

and structural genomics. Nucleic Acids Res;31: 489-91. 15.

Tramontano A, Leplae R, Morea V (2001) Analysis and assessment of comparative

modeling predictions in CASP4. Proteins: Structure, Function, and Bioinformatics;45: 2238. 16.

Martí-Renom MA, Stuart AC, Fiser A, Sánchez R, Melo F, Šali A (2000)

Comparative protein structure modeling of genes and genomes. Annual review of biophysics and biomolecular structure;29: 291-325. 17.

Hénin J, Maigret B, Tarek M, Escrieut C, Fourmy D, Chipot C (2006) Probing a

model of a GPCR/ligand complex in an explicit membrane environment: the human cholecystokinin-1 receptor. Biophys J;90: 1232-40. 18.

Zhu J, Fan H, Periole X, Honig B, Mark AE (2008) Refining homology models by

combining replica‐exchange molecular dynamics and statistical potentials. Proteins: Structure, Function, and Bioinformatics;72: 1171-88. 19.

Marti-Renom M, Fiser A, Madhusudhan M, John B, Stuart A, Eswar N, et al. (2003)

Current Protocols in Bioinformatics. Current Protocols in Bioinformatics. John Wiley & Sons, Inc.

This article is protected by copyright. All rights reserved.

20.

Schonbrun J, Wedemeyer WJ, Baker D (2002) Protein structure prediction in 2002.

Accepted Article

Current opinion in structural biology;12: 348-54. 21.

Venclovas Č, Zemla A, Fidelis K, Moult J (2001) Comparison of performance in

successive CASP experiments. Proteins: Structure, Function, and Bioinformatics;45: 163-70. 22.

Raval A, Piana S, Eastwood MP, Dror RO, Shaw DE (2012) Refinement of protein

structure homology models via long, all‐atom molecular dynamics simulations. Proteins;80: 2071-9. 23.

Fan H, Mark AE (2004) Refinement of homology‐based protein structures by

molecular dynamics simulation techniques. Protein Science;13: 211-20. 24.

Baker D, Sali A (2001) Protein structure prediction and structural genomics.

Science;294: 93-6. 25.

Moro S, Deflorian F, Bacilieri M, Spalluto G (2006) Ligand-based homology

modeling as attractive tool to inspect GPCR structural plasticity. Current pharmaceutical design;12: 2175-85. 26.

Evers A, Klabunde T (2005) Structure-based drug discovery using GPCR homology

modeling: successful virtual screening for antagonists of the alpha1A adrenergic receptor. Journal of medicinal chemistry;48: 1088-97. 27.

Moraitakis G, Purkiss AG, Goodfellow JM (2003) Simulated dynamics and biological

macromolecules. Rep Prog Phys;66: 383. 28.

Straub JE, Rashkin AB, Thirumalai D (1994) Dynamics in rugged energy landscapes

with applications to the S-peptide and ribonuclease A. J Am Chem Soc;116: 2049-63. 29.

Caves LS, Evanseck JD, Karplus M (1998) Locally accessible conformations of

proteins: multiple molecular dynamics simulations of crambin. Protein Sci;7: 649-66.

This article is protected by copyright. All rights reserved.

30.

Papaleo E, Pasi M, Riccardi L, Sambi I, Fantucci P, Gioia LD (2008) Protein

Accepted Article

flexibility in psychrophilic and mesophilic trypsins. Evidence of evolutionary conservation of protein dynamics in trypsin-like serine-proteases. FEBS Lett;582: 1008-18. 31.

Zuckerman DM, Lyman E (2006) A second look at canonical sampling of

biomolecules using replica exchange simulation. Journal of chemical theory and computation;2: 1200-2. 32.

Shahlaei M, Madadkar-Sobhani A, Mahnam K, Fassihi A, Saghaie L, Mansourian M

(2011) Homology modeling of human CCR5 and analysis of its binding properties through molecular docking and molecular dynamics simulation. Biochimica et Biophysica Acta (BBA)-Biomembranes;1808: 802-17. 33.

Shahlaei M, Madadkar-Sobhani A, Fassihi A, Saghaie L (2011) Exploring a model of

a chemokine receptor/ligand complex in an explicit membrane environment by molecular dynamics simulation: the human CCR1 receptor. Journal of chemical information and modeling;51: 2717-30. 34.

Shahlaei M, Fassihi A, Papaleo E, Pourfarzam M (2013) Molecular Dynamics

Simulation of Chemokine Receptors in Lipid Bilayer: A Case Study on C–C Chemokine Receptor Type 2. Chem Biol Drug Des;n/a-n/a. 35.

Shahlaei M, Mousavi A (2014) A 3D Model for Human Melanocortin 4 Receptor

Refined with Molecular Dynamics Simulation. Journal of Reports in Pharmaceutical Sciences (J Rep Pharm Sci);3: 42-53. 36.

Peter Tieleman D, Berendsen HJ, Sansom MS (1999) An alamethicin channel in a

lipid bilayer: molecular dynamics simulations. Biophysical journal;76: 1757-69. 37.

Bachar M, Becker OM (2000) Protein-induced membrane disorder: a molecular

dynamics study of melittin in a dipalmitoylphosphatidylcholine bilayer. Biophysical journal;78: 1359-75.

This article is protected by copyright. All rights reserved.

38.

Capener CE, Shrivastava IH, Ranatunga KM, Forrest LR, Smith GR, Sansom MS

Accepted Article

(2000) Homology modeling and molecular dynamics simulation studies of an inward rectifier potassium channel. Biophysical journal;78: 2929-42. 39.

Elmore DE, Dougherty DA (2001) Molecular Dynamics Simulations of Wild-Type

and Mutant Forms of the< i> Mycobacterium tuberculosis MscL Channel. Biophysical journal;81: 1345-59. 40.

Palucki BL, Park MK, Nargund RP, Ye Z, Sebhat IK, Pollard PG, et al. (2005)

Discovery of (2< i> S)-< i> N-[(1< i> R)-2-[4-cyclohexyl-4-[[(1, 1dimethylethyl) amino] carbonyl]-1-piperidinyl]-1-[(4-fluorophenyl) methyl]-2-oxoethyl]-4methyl-2-piperazinecarboxamide (MB243), a potent and selective melanocortin subtype-4 receptor agonist. Bioorg Med Chem Lett;15: 171-5. 41.

Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJ (2005)

GROMACS: fast, flexible, and free. J comput chem;26: 1701-18. 42.

Oostenbrink C, Villa A, Mark AE, Van Gunsteren WF (2004) A biomolecular force

field based on the free enthalpy of hydration and solvation: The GROMOS force‐field parameter sets 53A5 and 53A6. J comput chem;25: 1656-76. 43.

Berger O, Edholm O, Jähnig F (1997) Molecular dynamics simulations of a fluid

bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys J 72: 2002-13. 44.

Kandt C, Ash WL, Peter Tieleman D (2007) Setting up and running molecular

dynamics simulations of membrane proteins. Methods;41: 475-88. 45.

Hess B, Bekker H, Berendsen HJ, Fraaije JG (1997) LINCS: a linear constraint solver

for molecular simulations. J computat chem;18: 1463-72. 46.

Berendsen HJ, Postma JPM, van Gunsteren WF, DiNola A, Haak J (1984) Molecular

dynamics with coupling to an external bath. J chem phys 81: 3684.

This article is protected by copyright. All rights reserved.

47.

Darden T, York D, Pedersen L (1993) Particle mesh Ewald: An N⋅ log (N) method for

Accepted Article

Ewald sums in large systems. J Chem Phys 98: 10089. 48.

Patra M, Karttunen M, Hyvönen MT, Falck E, Lindqvist P, Vattulainen I (2003)

Molecular dynamics simulations of lipid bilayers: major artifacts due to truncating electrostatic interactions. Biophys J 84: 3636-45. 49.

van Aalten DM, Bywater R, Findlay JB, Hendlich M, Hooft RW, Vriend G (1996)

PRODRG, a program for generating molecular topologies and unique molecular descriptors from coordinates of small molecules. J Comput Aided Mol Des;10: 255-62. 50.

Humphrey W, Dalke A, Schulten K (1996) VMD: visual molecular dynamics. J mol

graph 14: 33-8. 51.

Kabsch W, Sander C (1983) Dictionary of protein secondary structure: pattern

recognition of hydrogen‐bonded and geometrical features. Biopolymers;22: 2577-637. 52.

Amadei A, Linssen A, Berendsen HJ (1993) Essential dynamics of proteins. Prot:

Struct Funct Bioinform;17: 412-25. 53.

Papaleo E, Mereghetti P, Fantucci P, Grandori R, De Gioia L (2009) Free-energy

landscape, principal component analysis, and structural clustering to identify representative conformations from molecular dynamics simulations: the myoglobin case. J Mol Graph Model 27: 889-99. 54.

Amadei A, Linssen A, Berendsen HJ (1993) Essential dynamics of proteins. Proteins:

Structure, Function, and Bioinformatics;17: 412-25. 55.

De Groot B, van Aalten D, Amadei A, Berendsen H (1996) The consistency of large

concerted motions in proteins in molecular dynamics simulations. Biophys J;71: 1707-13. 56.

Van Aalten DM, De Groot BL, Findlay JB, Berendsen HJ, Amadei A (1997) A

comparison of techniques for calculating protein essential dynamics. J Comput Chem;18: 169-81.

This article is protected by copyright. All rights reserved.

57.

Hess B (2002) Convergence of sampling in protein simulations. Physical Review

Accepted Article

E;65: 031910. 58.

Hess B (2000) Similarities between principal components of protein dynamics and

random diffusion. Physical Review E;62: 8438. 59.

Balsera MA, Wriggers W, Oono Y, Schulten K (1996) Principal component analysis

and long time protein dynamics. J Phys Chem;100: 2567-72. 60.

Hess B (2002) Convergence of sampling in protein simulations. Phys Rev E;65:

031910. 61.

Maisuradze GG, Leitner DM (2006) Principal component analysis of fast-folding λ-

repressor mutants. Chem Phys Lett;421: 5-10. 62.

Grossfield A, Feller SE, Pitman MC (2007) Convergence of molecular dynamics

simulations of membrane proteins. Proteins;67: 31-40. 63.

Maisuradze GG, Leitner DM (2007) Free energy landscape of a biomolecule in

dihedral principal component space: Sampling convergence and correspondence between structures and minima. Proteins;67: 569-78. 64.

Shu F, Ramakrishnan V, Schoenborn BP (2000) Enhanced visibility of hydrogen

atoms by neutron crystallography on fully deuterated myoglobin. Proc Natl Acad Sci USA 97: 3872-7. 65.

Hess B (2000) Similarities between principal components of protein dynamics and

random diffusion. Phys Rev E;62: 8438. 66.

Marianayagam NJ, Jackson SE (2005) Native-state dynamics of the ubiquitin family:

implications for function and evolution. J Roy Soc Interface;2: 47-54. 67.

Grubmüller H (1995) Predicting slow structural transitions in macromolecular

systems: conformational flooding. Physical Review E;52: 2893.

This article is protected by copyright. All rights reserved.

68.

Kosztin I, Barz B, Janosi L (2006) Calculating potentials of mean force and diffusion

Accepted Article

coefficients from nonequilibrium processes without Jarzynski’s equality. The Journal of chemical physics;124: 064106. 69.

Luchko T, Torin Huzil J, Stepanova M, Tuszynski J (2008) Conformational Analysis

of the Carboxy-Terminal Tails of Human< i> β-Tubulin Isotypes. Biophys J;94: 1971-82.

This article is protected by copyright. All rights reserved.

Accepted Article

Table 1 Summary information about the different replicates used in this study

Run name REP1 REP 2 REP 3 REP 4 REP 5 REP 6 REP 7 REP 8 REP 9 REP 10 REP 11 REP 12 REP 13 REP 14 REP 15 REP 16 REP 17 REP 18 REP 19 REP 20 REP 21 REP 1_2 REP 1_3 REP 1_4 REP 1_5 REP 1_6 REP 1_7 REP 1_8 REP 1_9 REP 1_10 REP 1_11 REP 1_12 REP 1_13 REP 1_14 REP 1_15 REP 1_16 REP 1_17 REP 1_18 REP 1_19 REP 1_20 REP 1_21

Duration of equilibrated time (ns) 3100 1700 4600 4200 3300 4200 6000 4500 8100 4100 6800 800 1900 2000 900 2000 6000 4600 2700 2300 7600 4800 9400 12700 16000 20200 26200 30700 38800 42900 49700 50500 52400 54400 55300 57300 63300 66000 68300 75900 83500

starting structure HM HM HM HM HM HM HM HM HM HM HM HM HM HM HM HM HM HM HM HM HM Concatenated trajectory including the first 2 replicas from HM Concatenated trajectory including the first 3 replicas from HM Concatenated trajectory including the first 4 replicas from HM Concatenated trajectory including the first 5 replicas from HM Concatenated trajectory including the first 6 replicas from HM Concatenated trajectory including the first 7 replicas from HM Concatenated trajectory including the first 8 replicas from HM Concatenated trajectory including the first 9 replicas from HM Concatenated trajectory including the first 10 replicas from HM Concatenated trajectory including the first 11 replicas from HM Concatenated trajectory including the first 12 replicas from HM Concatenated trajectory including the first 13 replicas from HM Concatenated trajectory including the first 14 replicas from HM Concatenated trajectory including the first 15 replicas from HM Concatenated trajectory including the first 16 replicas from HM Concatenated trajectory including the first 17 replicas from HM Concatenated trajectory including the first 18 replicas from HM Concatenated trajectory including the first 19 replicas from HM Concatenated trajectory including the first 20 replicas from HM Concatenated trajectory including the first 21 replicas from HM

This article is protected by copyright. All rights reserved.

Table 2 Interaction, Van Der Waals and electrostatic energies for selected residues of MC4R

Accepted Article

and studied ligand.

Residue

Before Final MD

Afetr Final MD

Interaction

Van Der Waals

Electrostatic

Interaction

Van Der Waals

Electrostatic

Energya

Energy

Energy

Energy

Energy

Energy

Tyr 268

-5.951

-3.255

-2.698

1.988

-0.020

2.009

Val 179

-7.972

-0.440

-7.531

1.807

-0.003

0.810

His 264

-9.122

-1.317

-7.805

0.843

-0.115

0.959

Ile 182

-5.152

-3.676

-1.475

7.059

-4.404

0.007

Leu 288

-5.754

-0.104

-5.649

-2.249

-1.81

-2.187

Ser 180

-1.235

-1.057

-0.178

0.000

0.000

0.000

Thr 178

-17.557

-4.946

-12.611

-1.252

-0.009

-1.513

Ala176

-4.407

-3.707

-0.700

-4.450

-2.759

-1.690

Asn 285

-4.636

-1.474

-3.161

-1.731

-1.328

-0.403

His 283

-4.471

-4.143

-0.329

-2.144

-1.459

-0.685

Asn 274

-4.772

-0.043

-4.772

-6.033

-0.008

-6.024

Tyr 276

-5.835

-1.448

-4.387

-2.332

-0.2.336

-2.095

Cys 277

-2.364

-0.034

-2.330

-2.275

-4.683

-2.270

Asn 74

-1.735

-1.412

-1.721

-1.223

-0.888

-0.334

Ala 175

10.041

12.026

-1.984

-2.414

-1.473

-0.940

Ser 282

29.038

26.799

2.238

-13.564

-2.267

-11.297

Leu 265

1.777

-0.240

2.017

-0.832

-0.0001

-0.832

Trp 174

12.171

-0.927

13.098

-2.020

-0.205

-1.815

Thr 162

44.345

45.568

-1.222

-17.706

-0.023

-17.682

a

All of the enrgies are Kcal/mol

This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article

Figure 7 RMSDs of compound 9 relative to MC4R derived by molecular dynamic calculation. The RMSD of MC4R relative to itself is also presented in the blue curves

Figure 8 The binding mode extracted from the results of the MD simulation as (A) 3D representation, and (B) 2D representation.

This article is protected by copyright. All rights reserved.

Accepted Article

Figure 9 RMSF of Cα atoms from their time-averaged positions for free (blue) and bound (red) MC4R.

H H

N HN

O

O

NH NH

N O

F

Scheme 1 2D structure of agonist applied in this study

This article is protected by copyright. All rights reserved.

A Conformational Analysis Study on the Melanocortin 4 Receptor Using Multiple Molecular Dynamics Simulations.

Taking into account the uncertainties involved in 3D model of biomolecule developed by homology modeling (HM), it is important to opportunely validate...
1MB Sizes 0 Downloads 7 Views