Investigation of the Binding Site of CCR2 using 4-Azetidinyl-1-aryl-cyclohexane Derivatives: A Membrane Modeling and Molecular Dynamics Study

Abstract

Chemokine receptor (CCR2) is a G protein-coupled receptor that contains seven transmembrane helices. Recent pharmaceutical research has focused on the antagonism of CCR2 and candidate drugs are currently undergoing clinical studies for the treatment of diseases like arthritis, multiple sclerosis, and type 2 diabetes. In this study, we analyzed the time dependent behavior of CCR2 docked with a potent 4-azetidinyl-1-aryl-cyclohexane (4AAC) derivative using molecular dynamics simulations (MDS) for 20 nanoseconds (ns). Homology modeling of CCR2 was performed and the 4AAC derivative was docked into this binding site. The docked model of selected conformations was then utilized to study the dynamic behavior of the 4AAC enzyme complexes inside lipid membrane. MDS of CCR2-16b of 4AAC complexes allowed us to refine the system since binding of an inhibitor to a receptor is a dynamic process and identify stable structures and better binding modes. Structure activity relationships (SAR) for 4AAC derivatives were investigated and reasons for the activities were determined. Probable binding pose for some CCR2 antagonists were determined from the perspectives of binding site. Initial modeling showed that Tyr49, Trp98, Ser101, Glu291, and additional residues are crucial for 4AAC binding, but MDS analysis showed that Ser101 may not be vital. 4AAC moved away from Ser101 and the hydrogen bonding between 4AAC and Ser101 vanished. The results of this study provide useful information regarding the structure-based drug design of CCR2 antagonists and additionally suggest key residues for further study by mutagenesis.

Introduction

The chemokine receptor family includes ~20 G-proteincoupled receptors, which play central roles in leukocyte migration and activation.1 Specific family members are also involved in viral entry and angiogenesis. Given this diverse range of functions, Chemokine receptors have been targeted as potential points for pharmaceutical intervention on the basis of chemokine receptor antagonism for the treatment of diseases as diverse as asthma, rheumatoid arthritis, multiple sclerosis, atherosclerosis, cancer, and HIV infection.2 Chemokines are small (8-10 kDa) water-soluble proteins consisting of 340-380 amino acid residues that play key roles in immuno-modulation and host defense. Chemokines selectively recruit monocytes, neutrophils, and lymphocytes to sites of vascular injury and inflammation,3-5 and different chemokines produce different leukocyte responses depending on the complementary natures of their chemokine receptors. 67 The basic feature of inflammation is the recruitment of leukocytes by tissue, and this process is mediated mainly by chemokines (chemotactic cytokines) via their receptors. The chemokine super family can be categorized into four groups (CC, CXC, CX3C, and C), according to the number and spacing of conserved cysteines in their amino acid sequences.28-10 Apart from their well-recognized roles in leukocyte recruitment, some chemokines and chemokine receptors play crucial roles during other cellular functions, such as, activation, proliferation, and differentiation.28-10 Specific family members are also involved in viral entry and angiogenesis.2

CCR2 is a G protein-coupled receptor that binds multiple ligands (macrophage chemoattractant proteins), such as, CCL2 (MCP-1), CCL8 (MCP-2), CCL7 (MCP-3), and CCL13 (MCP-4). The relative contributions made by these ligands to CCR2-mediated in vivo functions remain to be elucidated. 910 Of these ligands, MCP-1 has been studied most comprehensively and CCR2 is considered the elite receptor of MCP-1.910 The indispensable role of the CCR2/MCP-1 axis in the tissue recruitment of monocytes/macrophages has been established by CCR2 and MCP-1 genetic ablation11-13 and pharmacologic inhibition studies.1415

Both CCR2 and MCP-1 KO mice are apparently healthy but exhibit impaired abilities to recruit monocytes/macrophages to sites of inflammation11-13 and to produce cytokines such as TNF-α.1617 This impaired monocyte recruitment in CCR2 KO mice treated with a CCR2 antagonist can be explained by reduced release of inflammatory monocytes from bone marrow to blood17 and migration of blood monocytes from blood to inflamed tissues.11-13 Notably, the CCR2/ MCP-1 axis is not important for monocyte adhesion but it is important for migration into inflamed tissues, such as, atherosclerotic arteries.18 Consistent with its key role in monocyte trafficking, CCR2 has been shown to drive inflammation in a number of animal models of diseases such as: RA, CD, and cardiovascular diseases, including atherosclerosis and AIH, and transplant rejection. Recent pharmaceutical research has focused on the antagonism of CCR2. Several molecules including CCX915 (Chemo Centryx), INCB3284 (Incyte and Pfizer), MK0812 (Merck), MLN1202 (Millennium Pharmaceuticals), and MCP-1 antagonist (Telik) are currently under clinical study for roles in disease conditions, such as, arthritis, multiple sclerosis, and type 2 diabetes. Accordingly, because of its importance, CCR2 is considered an attractive target in the pharmaceutical field.19

Computational modeling has become an essential tool for guiding and enabling rational decision making during hypothesis- driven biological research. Knowledge of the 3D structures of receptors is important for understanding the molecular mechanisms underlying the diseases caused by mutations. Though homology models cannot provide protein structural data comparable with that determined experimentally, they often provide rational alternatives that aid the understanding of the binding modes of lead molecules. The crystal structure of CCR2 has not yet been reported, and in the absence of firm structural data, ligand-based approaches have proven to be particularly useful for studying G proteincoupled receptors (GPCR).20 However, few studies have been undertaken to examine CCR2 modeled structures.2122

Most structure/activity relationship (SAR) studies have been focused on potent CCR2 antagonists.2324 We have developed 3D quantitative SAR (3D-QSAR) models and used them to perform QSAR studies on CCR2 antagonists.25 However, no study has yet been performed to investigate the time dependent behavior of CCR2 and potent antagonists using molecular dynamic (MD) simulations, which prompted us to initiate this analysis.

In this study, the SARs of 4-azetidinyl-1-aryl-cyclohexanes (4AACs; potent CCR2 antagonists) were investigated using MDS (CCR2-potent 4AAC derivative) studies. A homology model of CCR2 was constructed by using CXCR4 (PDB code: 3ODU) as a template26 and the developed model was then refined. Furthermore, a potent 4AAC derivative was docked inside the proposed CCR2 binding site and the final CCR2-ligand complex was subjected to a 20 ns MD simulation in a DPPC/TIP3P membrane/water environment. The results from this analysis are useful for evaluating the stability of the binding site and for identifying perturbations in interaction profiles undetectable in docking studies.

 

Materials and Methods

Sequence Analysis of CCR2. The entire amino acid sequence of human CCR2, which is composed of 374 amino acids, was retrieved from the UniProt database (accession number P41597). In order to identify an adequate template for modeling CCR2 chemokine receptor, a BLAST algorithm2728 (Basic local alignment search tool for proteins) search was carried out against the protein data bank (PDB).29 The top template obtained was then selected for sequence analysis using ClustalW 2.030 with default parameters.

Homology Modeling of CCR2. A number of homologous structures were identified as templates in the PDB. CXCR4 was reported as the top template. CXCR4 structure was used as the template structure and a sequence alignment was performed for modeling. Using the sequence alignment obtained against the top template structure, a homology model of CCR2 was generated using the Modeller9v4 program, 31-33 which computed a model composed of non-hydrogen atoms by comparing the sequence to be modeled against known related structures. A 3D model was obtained by optimization of a molecular probability density function (PDF) using a variable target function procedure in Cartesian space by employing conjugate gradients and MD with simulated annealing. Fifty such 3D models were generated for further computational study. Out of the generated models, the model with both the lowest Molpdf (molecular probability density function) score and the lowest root mean square deviation (RMSD) was selected. The best-modeled structure of CCR2 was prepared using a structure preparation tool, the Biopolymer module within SYBYL8.1.34 This model was refined by performing energy minimization in a vacuum assumption to relax from the strain in the model. A steepest-descent algorithm and Gromos53a6 force field were used in this refinement procedure. The selected model was further validated using PROCHECK,35 ERRAT,36 and ProSA (https://prosa.services.came.sbg.ac.at/prosa.php).

Dataset Used in This Study. A series of potent CCR2 antagonists and their biological activities (IC50 values) have been previously reported by Zhang et al.37 The geometries of ligand molecules were optimized using the Tripos force field methodology 38 using the distance-dependent dielectric and Powell's conjugate gradient method. 39 Partial atomic charges were applied using Gasteiger-Hückel charges.40 One of the most active compounds in the series reported by Zhang et al., compound 16b of 4AAC (IC50 = 5 nM, pIC50 = 8.3) was selected as the template molecule (Table 1 and Table 2). The IC50 (nM) values of molecules used in this study were converted to corresponding pIC50 = (-log IC50) values.

Table 1.Structure and biological values of 4-Azetidinyl-1-arylcyclohexane derivatives

Table 2.Structure and biological values of 4-Azetidinyl-1-arylcyclohexane derivatives

Binding Site Construction and Docking Analysis. The Autodock 4.0 program was utilized for making the docking calculations. Autodock uses the Lamarckian genetic algorithm (LGA) and is considered to be the best method in terms of the accuracy of its structural predictions and its ability to identify lowest energy structures.41 Hydrogen atoms and the active torsions of ligand were assigned using Autodock tools (ADT). The binding site for the receptor structure (CCR2) was assigned using previously published results. As it is well from previous reports that Glu291 is the crucial residue for ligand binding, we kept this residue as the starting point and the binding site was extended (by 5 Å residues) up to 5 Å and modeled to guide the ligand molecule en route for possible orientation in the binding site. Autogrid was employed to generate grid maps around the active site using 60 × 60 × 60 points and a grid spacing of 0.375 Å. The docking parameters modified from defaults were: the number of individuals in the population (set at 150), the maximum number of energy evaluations (set at 2,500,000), the maximum number of generations (set at 27,000), and the number of GA runs (set at 50). The final structures were clustered and ranked according to the Autodock scoring function. Finally two conformers were selected according to the scoring function, top cluster and using crucial residues determined by previous mutational studies and experimental studies.

Setup of the System [CCR2-ligand] in Bilayer Environment. To study the time dependent behavior of CCR2-16b of 4AAC complexes, MDS was done using the GROMACS simulation package42 in an explicit phospholipids bilayer. Protein, ligand, lipid and water molecules were used as components for the simulation. The GROMOS96 force field42 was used. The lipid bilayer was developed using specific topology files as described by Tieleman43 (http:// moose.bio.ucalgary.ca), was also used in the present study.

To obtain a better starting structure, the InflateGRO script from Prof. Tieleman’s website (http://moose.bio.ucalgary.ca)44 was used to pack lipids around an embedded protein. The starting point with this method is a pre-equilibrated bilayer into which the protein ligand complex has already been “inserted” in the lipid bilayer but fully overlaps with lipids. The system cannot be energy minimized due to extreme overlap. The bilayer was expanded using a scaling factor of 4 and a distance cut-off of 14 Å. Overlapping lipid molecules was deleted. A grid size of 5 Å was used for estimating the area per lipid. Energy minimization was done using strong position restraints to ensure that the structures (Proteinligand) don’t change at all. Energy minimization was followed by compression using a scaling factor of 0.95. The compression and energy minimization steps were repeated until the area per lipid converges to the reference value (0.62-0.64) for the lipid species used.

The system was then solvated with TIP3P water box45 by increasing the van der Waals radius of carbon atoms from 0.15 to 0.5, to prevent water molecules from being placed in the hydrophobic section of the bilayer. Once the system got solvated, the van der Waals radius of C was again set to 0.15 Å. The system was further neutralized by adding 27 Na+ and 37 Cl- counter ions with a concentration of 0.15 M. Periodic boundary conditions were applied and the ligand topologies as well as the parameters were obtained from PRODRG server.46 The parameterization methodology implemented in the PRODRG program is widely used in MD simulations of drugs. However, the applicability of PRODRG charges for molecular simulations has been questioned in the literature, because PRODRG doesn’t reproduce topologies in the force field due to inconsistent charges and charge groups.47 The partial charges are crucial in particular for a reasonably accurate description of non-covalent interactions due to the long-range nature of electrostatic interactions and PRODRG underestimates charges for some atom types. To overcome this issue, we cross-checked the topologies and we manually assigned charges.

After energy minimization, the system was subjected to a short NVT equilibration phase and was followed by a longer NPT phase. Generally a short NVT equilibration phase is followed by a longer NPT phase. We were dealing with a heterogeneous system [water and DPPC] as solvents, and such heterogeneity requires a longer equilibration process. Water has to reorient around the lipid head groups and any exposed parts of the protein and the lipids have to orient themselves around the macromolecule. During equilibration, the macromolecule (protein-ligand) was position restrained so that it allows the solvents to equilibrate around our protein structure without any structural changes in the protein. NVT ensemble was done at a constant temperature of 310 K for a period of 100 picoseconds (ps) and followed by NPT ensemble for a period of 1 ns. A modified Berendsen coupling scheme was employed in both the ensembles and the particle mesh Ewald (PME)48 method was used to calculate long-range electrostatics. The Parinello and Rahman coupling scheme was used as barostat for pressure coupling and semi-isotropic pressure coupling was applied with reference pressure (1.0 bar), which is intended for membrane simulations. All bond lengths were constrained using the LINCS algorithm49 and SETTLE algorithm was used to constrain the geometry of water molecules.50 The Final production run was performed for a period of 20 ns.

 

Results

Sequence Analysis of CCR2. A BLAST search revealed CXCR4 (PDB code: 3ODU) to be the top template. The sequence identity between the template (3ODU) and the query sequence (CCR2) was 35％ and the identity between the active sites of the template and query was 60％, which was encouraging. Whereas, the identity between the traditional bovine rhodopsin and CCR2 was found to be 21％ and the sequence identity with β2-adrenergic receptor was found to be 24％. Analysis showed high levels of homology between the target (CCR2) and template (CXCR4) sequences and better than that of traditional bovine rhodopsin and the more recent β2-adrenergic receptor templates. The reason for the high homology, overall and with the active site in particular was that the template sequence was a close homologue (CXCR4). A crucial step in the modeling procedure was to obtain acceptable alignment between target and template sequences. The alignment obtained using ClustalW 2.0 is shown in Figure 1. The E-value represents the number of different alignments with scores equivalent to or better than the scores expected to occur for a random database search. Generally, a lower E-value indicates that alignment is real and not due to chance. The expectation value (E-value) for CCR2 was found to be 2e-33.

Figure 1.Alignment obtained between the query (CCR2) and template (PDB code: 3ODU) sequences for modeling. Identical residues are marked (*), similar regions are marked (:). Secondary structure (TM domains) of the CCR2 receptor is indicated below the sequence in a cylindrical shape.

Homology Modeling of CCR2. CXCR4 (A chain) was used to develop the 3D models and a modeler program was used to derive 3D-models of CCR2. The unique feature of GPCR’s is the presence of seven transmembrane helices. As CCR2 is also a GPCR, we examined whether the seven transmembrane (TM) helixes were properly transformed in the models according to that of the template (PDB code: 3ODU) structure. Fifty models were developed for CCR2 and eventually the model with the lowest MolPdf as well as the lowest RMSD compared with the template structure was selected for further computational analysis. A key structural factor in CC chemokines is the presence of disulphide bridges between conserved cysteines residues and more than 90％ of the members of the GPCR super family have conserved disulfide bridges. Since the template structure also has the disulfide bridges in it, we checked whether the selected model has the disulfide bridges in it. As in the template, we found that the disulfide bridges were produced between Cys32-Cys277 and Cys113-Cys190 of the selected CCR2 model. The model (CCR2) selected was further refined by simple energy minimization in vacuum and the energy minimized CCR2 model using MDS is shown in Figure 2.

Figure 2.Homology model (green) of CCR2 obtained and its refinement by MDS superimposed on the template structure CXCR4 (red).

The model was further validated using a Ramachandran plot of PROCHECK35 to visualize the backbone dihedral angles ψ against φ of amino acid residues in the selected protein structure. The Ramachandran plot for the CCR2 model showed that most of the residues were in the favored (92.9％) and in the additionally allowed (7.1％) region. It is evident that, almost all models developed using modeling software pass the PROCHECK validation regardless of model quality. This is because PROCHECK uses the criteria of visualizing backbone dihedral angles ψ against the φ angles of amino acid residues in the protein structure, and thus, is more relevant for X-ray crystallographic studies. When the template and target structures are dissimilar, more caution should be taken and additional validation criteria must be used to determine the quality of the developed model.

Accordingly, additional parameters, such as, ERRAT36 and ProSA https://prosa.services.came.sbg.ac.at/prosa.php) energy plots were used to check the quality of the model. ERRAT plot is a program for verifying protein structures. Error values are plotted as a function of the position of a sliding 9- residue window. The error function is based on the statistics of non-bonded atomic interactions between atom types. Models with higher ERRAT scores are of higher quality and the ERRAT score of CCR2 was found to be 87.04. Similarly, we further validated our models using ProSA, which evaluates the energy of the structure using distance pair potentials. Residues with a negative ProSA score confirm model reliability. The ProSA energy scores for the template was −2.34 and for the model (CCR2) was −2.54. High ERRAT scores and low ProSA energy scores indicate that models are of high quality. Overall, our results indicate the selected models were satisfactory. Ramachandran and ProSA energy plots for refined CCR2 model are shown in Figure 3(a) and Figure 3(b). This refined and validated model of CCR2 was used for docking analyses.

Figure 3.Validation of generated model. (a) Ramachandran plot of the developed CCR2 model refined by MDS. Different color codes indicate most favored (red), generously allowed (dark yellow), additionally allowed (light yellow), and disallowed (white) regions. (b) ProsA energy plot for the developed CCR2 model refined by MD simulation.

Binding Site of CCR2. A putative binding site of CCR2 was determined as described previously.62122 Previous receptor homology modeling suggests that CCR2 antagonists bind to an extended pocket bounded by TM2, TM3, TM5, TM6, and TM7.51 Furthermore, it has been proposed in mutagenesis studies that Glu291 of TM7 is an important residue in the binding pocket.621 With knowledge of these previously published results, the binding pocket was determined. The binding pocket was found to be composed mainly of residues Leu45, Tyr49, Trp98, Ser101, Ala102, Tyr120, His121, Tyr124, Phe125, Met205, Arg206, Asn207, Trp256, Tyr259, Gln288, Glu291, Thr292, and Met295, which compares well with previous studies.621 The residues in the binding pocket that guide docking are shown in Figure 4(a) (front view) and Figure 4(b) (top view).

Figure 4.Proposed binding pocket of CCR2 for docking simulation. (a) Front view of the proposed binding pocket (CCR2). Binding site residues are colored based on atom types. (b) Top view of the CCR2 binding pocket. Transmembrane (TM) helices are colored green, whereas the constructed binding pocket residues were color-coded with different colors depending on atom types. All TM regions are labeled in red on the tops of helices. Figure generated using the Pymol program (http://www.pymol.org).

Docking Studies of Highly Active CCR2 Antagonist (compound 16b of 4AAC). A series of highly active CCR2 antagonists have been previously reported by Zhang et al.37As the binding pocket projected in this study is similar to that of the previously published reports by Mirzadegan et al. and Marshal et al.,622 one of the potent CCR2 antagonist reported by Zhang (IC50 = 5 nM, pIC50 = 8.3) was docked into the binding site. Fifty conformations were generated and the top ranking conformational clusters from our docking study were evaluated. The conformations within the top ranked cluster were then selected for analysis. It was observed that the selected cluster consists of 14 conformations within the cluster. It is well known that the top scoring binding mode is not always the correct binding mode in docking runs, even when ligands were docked into their own x-ray crystal structure. In cross docking experiments, the top scoring binding mode is rarely the native binding conformation. Docking against homology models makes it even more unlikely that the highest scoring pose is actually the correct pose. So, we carefully analyzed the conformations within the top cluster and two conformations were selected based on scoring function, interaction with crucial amino acid residues and a salt bridge contact with Glu291.

The compound 16b of 4AAC bound with a binding energy of −8.95 kcal/mol and an intermolecular energy of −10.34 kcal/mol. In addition, the ligand established crucial interactions with important residues in the CCR2 binding site. The basic nitrogen in the azetidinyl ring formed an electrostatic interaction (i.e., salt bridge contact) with the crucial and conserved Glu291, and the distance between the glutamic acid residue and the basic nitrogen was 4 Å. Furthermore, it was observed that the ligand formed hydrogen bond interactions. The hydrogen atom of amine nitrogen next to the azetidinyl ring hydrogen bonded with Tyr49 and amine nitrogen close to the indazole ring of 4AAC hydrogen bonded with Trp98 at a distance of 2.4 Å. Moreover, the nitrogen of the indazole ring hydrogen bonded with Ser101 at a distance of 2.9 Å. In addition to the hydrogen bonds and salt bridge interaction, we looked for residues that are in the vicinity of ligand. We found that residues, such as, Lys34, Leu45, Trp98, Ala99, Ala102, Trp106, Phe116, Val189, Cys190, Asn199, Thr203, Arg206, Asn266, Gln269, Asp284, Thr287, and Gln288 are in the vicinity of ligand (within 4 Å from the ligand molecule). The docked mode of the ligand molecule and its interaction with active site residues is shown in Figure 5.

Figure 5.Docked mode of compound 16b inside the proposed binding pocket of CCR2. Front view of the docked conformation of 16b. Binding site residues (carbon cyan, nitrogen blue and oxygen red) and the ligand molecule (carbon yellow, nitrogen blue and oxygen red) are colored based on atom types. Figure generated using the Pymol program (http://www.pymol.org).

In addition we also selected another binding mode to validate the stability of the binding mode using MDS. This would also be critical while analyzing the SAR’s. It was observed that the ligand bound with a binding energy of −7.88 kcal/mol and an intermolecular energy of −8.4 kcal/ mol. As the salt bridge distance between the ligand and Glu291 is crucial, we selected the pose with the larger distance between them. Moreover, this binding mode shared a similar ligand space, which is more or less similar to the previous mode with only a slight variation in its orientation of substitutions. The docked mode of the ligand molecule and its interaction with active site residues is shown in Supplementary Figure S1.

Despite the fact that the binding modes selected in the present study and the interaction between 4AAC and CCR2 are similar to previous results, the present study may be limited because our selection of binding mode was based on the top cluster, previous mutagenesis results, and binding energy. This study can also be done manually by assigning the salt bridge contact as the starting point of interaction between the ligand and Glu291 residue as it is crucial for activity.621 However, our selected docked modes obtained the essential salt bridge contact with Glu291 with additional hydrogen bonds and also identified residues in the TM1, TM2, TM3 and TM7 with additional residues from TM5 and TM6 as similarly proposed by Carter and Tebben.51 These findings encouraged the reliability of our results.

MDS of Receptor-Ligand-Membrane Complexes. Since CCR2 is a membrane protein, an accurate computational representation of CCR2 receptor should be embedded in a membrane aqueous environment. In vivo binding of inhibitor to a receptor is a dynamic process, whereas the docked orientation of ligand considers the receptor to be rigid. To overcome this bias, MDS of CCR2-ligand (selected poses from docking) complex embedded in a lipid layer was done to mimic a real biological environment. Accordingly, we analyzed a membrane MD simulation of the docked model of one of the 4AAC derivatives complexed with CCR2. This strategy employed here helped us to identify the most stable and low energy conformation of CCR2-4AAC complex. We assume the selected conformation of ligand inside the receptor could be the bioactive one. As the rest of the molecules in the dataset have similar core structure and variations in only substitutions, by assuming the stable and average conformation of ligand-receptor complex as the bioactive one, it will allow us to understand the SAR’s of 4-azetidinyl-1-aryl-cyclohexane derivatives. Furthermore, SAR’s of diverse CCR2 ligands could be analyzed in the same manner to obtain additional insights of the CCR2 binding site. In addition, the model obtained and relaxed in a membrane environment could be more appropriate for the structure-based design of novel CCR2 antagonists.

Setting up of CCR2-16b of 4AAC-lipid for the Production Run. The selected docked conformations of 16b of 4AAC complexed with CCR2 were implanted in a rectangular box composed of a DPPC bilayer solvated with water molecules. Overlapping lipids were deleted. It has been observed that, three lipid molecules from the upper leaflet and three lipid molecules from the lower leaflet were deleted using InflateGRO. After the deletion of overlapping lipids, the protein-ligand-lipid complex was then shrunk until it reached an area per lipid of 0.65 Å2, which is closest to the reference value for DPPC lipids. The protein-ligand-lipid complex was then solvated by adding water molecules and appropriate counter ions [27 Na+ and 37 Cl−] were added to neutralize the system. To overcome structural artifacts, the system was then energy minimized. The system composed of protein, ligand, lipids, water molecules and ions was then equilibrated and a final production run of 20 ns was performed. The tightly packed, prepared system (CCR2-16b of 4AAC-lipids-water molecules) included in the production run is shown in Figure 6.

Figure 6.Side view of CCR2-ligand-membrane-aqueous environment forproduction runs of 20 ns MDS. Seven transmembrane helices are represented as cartoons, the ligand molecule as a surface, lipid molecules as lines, and water molecules as sticks. Figure generated using the Pymol program (http://www.pymol.org).

Structural Analysis of CCR2-4AAC Complexes Throughout Production Runs. A number of parameters needed to be examined to substantiate the physical stability of the CCR2 model. The stability of the system was examined structurally and energetically during MD simulation as a function of time. The total energy and RMSDs of the protein structure versus the initial structure were used as parameters to monitor the stability of the model. The total energy plot of the system indicated that the total energy decreased gradually until 3 ns and then stabilized at this level. Trajectory based analysis revealed the total energy of the system decreased from −1.81e+05 KJ/mol to −2.16e+05 kJ/mol. It was observed that most of the structures are around −2.14e +05 kJ/mol, indicating the system were energetically stable. The CCR2 model was also evaluated on the basis of structural stability using RMSD calculated by structural variations with respect to time; the initial period was considered as period of equilibration. A gradual rise in RMSD until 0.38 nm (1.17 ns) was seen and then it decreased gradually to 0.34 nm till 5 ns. This was then followed by a plateau, which for most of the structures occurred at around 0.34(± 0.02) nm till 20 ns. These results indicated the structural stability of the model. The RMSD of ligand molecule and the movement of ligand atoms throughout the simulation were also analyzed. It was observed that the initial conformation gradually moved to 0.2 (+0.2) nm till 2 ns, and then gradually decreased to 0.16 ns around 3.9 ns. This was then followed by a slight increase in its RMSD, which reached 0.2 nm around 6 ns and maintained around 0.2 (± 0.02) nm till 20 ns. These results suggest that this ligand is stable throughout the simulation. The total energy plot and the RMSD plot of the receptor and ligand molecules are shown in Figures 7(a), 7(b), and 7(c), respectively.

Figure 7.(a) Root mean square deviation (RMSD) plot for CCR2 Cα from initial structures throughout the 20 ns simulation as a function of time. (b) Total energy plot of the MD simulation and variations in system total energy over 20 ns. (c) Root mean square deviation (RMSD) plot for DRG as ligand throughout the 20 ns simulation.

It was discussed in our initial docking study that the selected binding mode was based on previous results and crucial interaction with active site residues. Moreover, further MD simulation study found that the selected binding mode was found to be stable structurally as well as energetically. However, it is well known fact that the top scoring binding mode is not always the correct binding mode. Though our selected binding mode was based on scoring function, salt bridge contact and crucial interactions, in order to overcome bias and to validate the selection of binding mode, another binding mode was implemented for 20 ns MDS. The RMSD plot of ligand molecule indicates that the ligand is moving throughout the production run and found to be unstable (Supplementary Figure S2(b)). Moreover the movement of ligand molecule makes the protein structure to be unstable and the RMSD continues to rise throughout the simulation (Supplementary Figure S2(a)). These results validate the selection of binding mode.

Binding Mode (CCR2-16b of 4AAC) and Trajectory Analysis. It was observed from the RMSD plot of ligand molecule that the ligand was moved from initial position inside the proposed binding site and it maintained around 0.2 ± 0.02 nm after the equilibration period. Our results indicate the practical relevance of MDS after initial rigid docking analysis inside the proposed binding site. The selected docked mode after initial rigid docking observed that the ligand established hydrogen bonding with Tyr49, Trp98 and Ser101 and moreover the salt bridge distance is 4 Å (Figure 5). However, the binding mode of 4AAC in the binding pocket of CCR2 (Figure 8), which were derived from the average low energy structure from MD simulation leads to the introduction of some new residues and some new interactions into the binding pocket compared with docking. There are some differences in interaction of 4AAC and these differences somewhat affect the position of ligand molecule in the binding pocket of CCR2. We did monitor the detailed protein-ligand interaction over time through minimum distance analyses between important residues, and the salt bridge distance between the crucial Glu291 and the basic nitrogen in the azetidinyl ring of ligand molecule.

Figure 8.Front view of the docked mode of compound 16b inside the proposed binding pocket of CCR2 after MDS. Binding site residues (carbon cyan, nitrogen blue and oxygen red) and the ligand molecule (carbon yellow, nitrogen blue and oxygen red) are colored according to atom type.

In the case of the binding mode of 16b inside the CCR2 binding pocket (average low energy structure) derived from MD simulation, some new residues were noted in vicinity of the ligand. Our analysis of the 20 ns MDS of the interaction between CCR2 and the 16b of 4AAC revealed additional residues, such as, Ala103, Asn104, Tyr120, His121, Tyr188, Val189, Cys190, Pro192, Tyr193, Phe194, Asn207, Ty259 and Ile263 had moved into the vicinity of ligand (within 4 Å). More importantly, the hydrogen bonding between the ligand and Ser101 obtained from initial docking vanished after MDS, suggesting that this bonding might not be crucial and not in vicinity of the ligand. However, Tyr49, which hydrogen bonds with 16b of 4AAC, was found to stay in the vicinity of the ligand, which confirms the importance of this residue as suggested by previous mutational studies. The hydrogen bond between Trp98 and 16b after initial docking got eliminated because of the ligand movement; however this residue stayed in vicinity throughout the simulation, which confirms the importance of the residue. The distance between residues (Tyr49, Trp98 and Ser101) and 16b was monitored throughout the simulation and shown in Figure 9.

Figure 9.Minimum distances between the 16b and CCR2 active site residues surrounding 4 Å were computed using MD simulation trajectories (a) Tyr49 and 16b, (b) Trp98 and 16b and (c) Ser101 and 16b. Minimum distances were calculated and plotted as a function of time.

Similarly, the initial distance between the glutamic acid residue and the basic nitrogen in the azetidinyl ring was found to be 4 Å (i.e., salt bridge contact). The distance between them in the average structure after MDS was found to be 3.5 Å. This revealed the presence of a strong electrostatic interaction. The results of the present study show the importance of this interaction for activity and complement previous results. The distance between Glu291 and the basic nitrogen of azetidinyl ring throughout the simulation is also plotted and is shown in Figure 10. It was also observed that Tyr120 and His121 were not in the vicinity of ligand in the initial docked mode. However, these residues are back in vicinity and suggested the importance of these residues and complements previous mutational studies (Figure 11).

Figure 10.Minimum distances between Glu291 (Oxygen atom) and 16b (basic nitrogen in azetidinyl ring) computed using MD simulation trajectories. Minimum distances were calculated and plotted as a function of time.

Figure 11.Minimum distances between (a) Tyr120 and 16b and (b) His121 and 16b. Minimum distances were calculated and plotted as a function of time.

The indazole ring of 16b moved inside the pocket lined by mostly hydrophobic residues such as Trp98, Trp106, Phe116 and Tyr188. This interaction stabilized the position of ligand inside the pocket. The hydroxyl group of the cyclohexyl ring was found to be hydrogen bonded with Cys190. Moreover, the N-methyl of the R1 substituted phenyl ring oriented towards Thr203, Arg206 and Asn207, implying the importance of these residues and suggesting the need for mutational studies for further analysis. The minimum distance between these residues and 16b was also monitored and shown in Figure 12. Residues Tyr259 and Ile263 were not in vicinity after initial docking moved towards the ligand and mutational studies on these residues could also be informative. The distance between these residues and 16b had also been monitored throughout the simulation and is shown in Figure 13. On the other hand, residues such as Lys34, Ala99, and Ser101 moved out of ligand vicinity during the course of simulation.

Figure 12.Minimum distances between (a) Thr203 and 16b, (b) Arg206 and 16b and (c) Asn207 and 16b. Minimum distances were calculated and plotted as a function of time.

Figure 13.Minimum distances between (a) Tyr259 and 16b and (b) Ile263 and 16b. Minimum distances were calculated and plotted as a function of time.

SAR Studies of 4-Azetidinyl-1-aryl-cyclohexane Derivatives. The average structure (CCR2-ligand) resulting from the MD run was further used to exploit the SAR’s of 4AAC derivatives. We considered the binding mode of 16b from 4AAC (IC50 = 5 nM, pIC50 = 8.3) that resulted from MDS as a bioactive conformer and the representative molecule for other 4-azetidinyl-1-aryl-cyclohexane derivatives reported by Zhang et al. We then compared the binding mode of this compound with those of the other molecules, since other molecules in the dataset are similar with varied substitution and SARs were explored.

Our trajectory analysis and the binding mode of 16b from the average structure indicated that the basic nitrogen is necessary for activity because of its electrostatic interaction with Glu291 in the binding site. It was also evident that compounds with hydrophobic substitution around R1 enhance activity and that an optimal bulky substitution at the R1 substituted phenyl ring is needed for activity and that the 1- hydroxyl groups of the cyclohexyl ring enhances activity. These SAR’s explain the greater activities of compounds 13b (IC50 = 6 nM, pIC50 = 8.2), 16e (IC50 = 6 nM, pIC50 = 8.2), 16h (IC50 = 7 nM, pIC50 = 8.1), and 16b (IC50 = 5 nM, pIC50 = 8.3). Compounds from series 14 (14a, 4b, 14c, 14e, 14f, 14j, 14k, 14m, 14n) were not active because they lacked desirable hydrophobic substitution at R1. Similarly compounds, such as 16q (IC50 = 480 nM, pIC50 = 6.3), 16r (IC50 = 1400 nM, pIC50 = 5.8) and 16t (IC50 = 10000 nM, pIC50 = 5) that lack hydrophobic substitution at R1 had poor activities. On the other side of the 16b from 4AAC, an indazole ring is needed and this hydrophobically interacts with Trp98. The indazole ring is desirable to access the hydrophobic residues Ala102, Ala103, Trp106, Phe116, Tyr188, and Val189 in the CCR2 binding pocket. Though compounds such as 13g, 13h, 14q, 14r, 14s, 14t, 14u, 14v, 16l, 16m substituted with hydrophobic groups around R1 they lacked activity. The reason is because of increased substitution around the indazole ring, which is sterically restricted to the receptor. Based on our analysis we concluded that a molecule with hydrophobic substitution around R1 and optimal substitution around indazole is desirable for potent activity.

Comparison of Other Potent Antagonists with the MD Simulated Model. The selected low energy model obtained after a prolonged run for 20 ns MDS was considered to be stable enough to allow meaningful analyses of CCR2 binding site and CCR2 antagonists. We considered this more relaxed structure as a valid one and compared the binding site of CCR2 with some of the potent CCR2 antagonists. Moreover, the orientation of crucial active site residues after MDS will assist us in identifying the possible orientation of some of the potent antagonists of CCR2. We analyzed the binding mode of selected CCR2 antagonists based on pharmacophore properties against some of the key residues in the active site. The importance of electrostatic interaction between Glu291 and the basic nitrogen in the ligand molecule was considered as the starting point. Although there could be some limitations with this approach, it is potentially useful in identifying the binding mode of some new antagonists.

Accordingly, we chose a representative set of potent CCR2 antagonists (Table 3). Initially, the pharmacophoric features of highly active compound 71(IC50 = 3.2 nM, pIC50 = 8.4) of the Teijin lead23 was compared with the binding site. The basic nitrogen in the pyrrolidine ring is probably forming an electrostatic interaction with Glu291. The 2,4-dimethyl- phenyl ring, which is highly hydrophobic, is expected to hydrophobically interact with hydrophobic pocked lined by residues Trp98, Ala102, Ala103, Trp106, Phe116, Tyr188, and Val189. The linear chain mediating the pyrrolidine ring and 2-NH2, 5-trifluoro methyl group might interact with Tyr120 and Tyr259. The 2-NH2 of the phenyl ring may access the pocket lined by residues Thr203, Arg206 and Asn207. Similarly, we examined the potent CCR2 antagonist, INCB3344 (IC50 = 10 nM, pIC50 = 8) originally identified by Brodmerkel et al.15 This ligand is somewhat similar to 16b of 4AAC and expected to occupy the same volume and probably with similar orientation. It was clear that Glu2916 is a crucial residue and it supposed to form an electrostatic interaction with the basic nitrogen of ligand molecule. Using these results as a starting point and from the perspectives of binding site residues it is possible to identify the probable binding mode of some other potent CCR2 antagonists.

Table 3.Representative set of CCR2 antagonists used to derive the binding mode from the perspectives of binding site

 

Discussion

CCR2 is a CC chemokine receptor and an important player in the trafficking of monocytes/macrophages. It also functions in other cell types and is relevant during the pathogeneses of various diseases such as arthritis, multiple sclerosis, and type-2-diabetes,52 and thus CCR2 is considered an important receptor. Several research groups have studies the SAR’s of CCR2 antagonists experimentally232437 or computationally.19202553 However, no previous study was conducted on the time dependent interactions between CCR2 and antagonists using MDS. This paucity of data prompted us to undertake the present study to determine the effects of active site residues in the vicinity of ligand. Specifically, the objective of this study was to analyze the behavior of CCR2- 4AAC complexes in a membrane environment and to characterize the reasons for antagonistic activity for some other antagonists (Teijin lead and INCB3344) from the perspective of binding site after MDS.

As a starting point of this study, the homology model of CCR2 was developed. Previously homology models of CCR2 were developed using the traditionally used bovine rhodopsin and the recently reported β2-adrenergic receptors as templates.2154-56 In this study we used the recently reported CXCR4 (PDB code: 3ODU) as a template and modeling was done because of the non-availability of crystal structure. A molecular docking study was then performed with 16b of 4AAC against the proposed binding site.62122 One conformation was selected and the results are in line with previously published results. Glu291 has been previously established by mutational studies to be a crucial residue in terms of the activities of CCR2 antagonists,621 and this is confirmed by our results, which show that Glu291 forms a salt bridge with antagonists at the active site. We also observed that hydrogen bonds were observed between the ligand and Tyr49, Trp98 and Ser101 of CCR2. These results are in line with previous mutagenesis studies.2157

However, the above-mentioned results were obtained from the docking study and are not conclusive. Binding of an inhibitor to a receptor is a dynamic process, and the rigid docking results obtained can only represent instantaneous states, and thus, only a glimpse of the possible molecular structures involved in docking. Nevertheless, dynamic ligand binding to receptor is more realistic and MDS studies allow stable structures and better binding modes to be identified. Thus, in the present study, we performed membrane MDS to better understand the SAR relations between CCR2 and 4- azetidinyl-1-aryl-cyclohexane derivatives. Therefore, we consider the relaxed model derived in a membrane environment is more appropriate for the structure based drug design of novel CCR2 antagonists.

A sampling protein dynamics in a single, 20 ns MD run is rather short by today standards. We point out the limitations arising from the small amount conformational sampling and are as follows. It could happen that the short simulation time does not allow interacting protein-ligand at their full performance. In some cases ligand binding could result in protein folding or structural changes that would be seen with long runtime simulation, which could not be visible with small simulation time. Short runtime simulation shows variable RMSD throughout simulation. In contrast long runtime simulation allows protein to stabilize its coordinates and proper folding. Different conformation of ligand can be observed in short as well as long time simulations. Inadequate conformational sampling can be done in short time simulations. Typically, the ligand may bind and stabilize only a subset of the many conformations sampled by its dynamic receptor, causing the population of all possible conformations to shift towards those that can best accommodate binding. Beyond these limitations, this study is helpful in identifying the crucial residues, binding site movement, position of ligand and interaction of ligand with crucial residues, which is quite encouraging.

The MDS analysis performed in the present study provides crucial insights into the character of the CCR2 binding site. The hydrogen bond between Ser101 and the 16b of 4AAC vanished and the ligand moved towards the hydrophobic residues Trp98, Ala102, Ala103, Trp106, Phe116, Tyr188, and Val189. However, the hydrogen between Tyr49 was retained and it was maintained during MDS. More importantly, the distance between Glu291 and the basic nitrogen in azetidinyl ring decreased from 4 Å (before MDS) to 3.5 Å (after MDS). Furthermore, the N-methyl of the R1 substituted phenyl ring oriented towards Thr203, Arg206 and Asn207, implying the importance of these residues and suggesting the need for further analysis by mutational studies. Overall, we found that Tyr49, Trp98, Tyr120, and Glu291 are crucial for activity, which concurs with previous mutagenesis results.62157 In addition, we found that Ala102, Thr203, Arg206, Asn207 and Tyr259 are probably crucially required for activity and we believe that this warrants further mutagenesis studies.

Our SAR studies showed that electrostatic interaction between 4AAC derivatives and Glu291 is highly desirable and that hydrophobic substitution of the cyclohexyl ring (R1 substitution) is required. Furthermore, optimal substitution around the indazole ring is required for better activity and bulky substitution hampers activity due to steric effects. In addition, we examined Teijin and INCB3344 (both potent CCR2 antagonists) with respect to binding by MDS. Our results could also be helpful in determining the binding mode of new antagonists.

We have been working on chemokine receptors using in silico methodologies.2558-61 In particular, we have developed 3D-QSAR models for CCR2 using ligand based and receptor guided methodologies.25 In this study, we used Teijin derivatives23 to derive CoMFA and CoMSIA models, and we found that Ala102 and Asn207 are probably crucial for activity. Recently, we compared the binding sites of CCR2 and CCR5 to facilitate the development of dual antagonists targeting both CCR2 and CCR5.58 In the present study, we found that Ser101 is important for CCR2 antagonism from our initial rigid docking study. However, the dynamic process showed that 4AAC derivative moved away from Ser101 and the hydrogen bond vanished. N-di-methyl substitution moved toward the region lined by residues Thr203, Arg206, and Asn207 whereas other residues, such as, Ala102, Thr203, and Arg206 stayed in the vicinity of the ligand (within 4 Å). These findings suggest the results obtained are realistic and that mutational studies on these residues are warranted.

 

Conclusion

In the present study, the time-dependent behavior of CCR2/4-azetidinyl-1-aryl-cyclohexane derivative complexes was studied by MDS. CCR2 is a key player in monocyte trafficking and also has cell-dependent functions relevant to the pathogeneses various diseases. In silico methodologies, such as, comparative modeling, docking, and MDS studies (20 ns) were performed in an explicit lipid layer environment. Homology modeling was performed using CXCR4 as template. Docking of 4-azetidinyl-1-aryl-cyclohexane derivative at the proposed CCR2 binding site was performed. MDS was then performed with the docked structure for 20 ns. MDS showed the importance of the real time dynamic analysis. Relations between the properties of potent antagonists and binding site residues were sought and the possible binding orientations of antagonists were deduced.

MD simulation results found that Tyr49, Trp98, Tyr120, and Glu291 are crucial for activity, which agrees with previous mutagenesis reports. Our results also identify the importance of Ala102, Thr203, Arg206 and Asn207 for CCR2 antagonism, and suggest their importance in further experimental mutational studies. The present simulation time is very short but we still obtained some conclusive results. However, to elucidate for in-depth structural folding as well as the mechanism of inhibition by non-peptide ligands, MD simulation times of one microsecond to one millisecond are required for membrane proteins. In future, computational mutational studies on newly identified residues could also be performed. In reality, the cell based in vitro studies as well as thermodynamic calculation on the newly identified residues are the only source to validate these computational results and confirm their importance in CCR2 antagonism. The results of the present study are especially valuable for further studies using site-directed mutagenesis and structure-based drug design.