Computational modeling of human coreceptor CCR5 antagonist as a HIV-1 entry inhibitor: using an integrated homology modeling, docking, and membrane molecular dynamics simulation analysis approach

Chemokine receptor 5 (CCR5) is an integral membrane protein that is utilized during human immunodeficiency virus type-1 entry into host cells. CCR5 is a G-protein coupled receptor that contains seven transmembrane (TM) helices. However, the crystal structure of CCR5 has not been reported. A homology model of CCR5 was developed based on the recently reported CXCR4 structure as template. Automated docking of the most potent (14), medium potent (37), and least potent (25) CCR5 antagonists was performed using the CCR5 model. To characterize the mechanism responsible for the interactions between ligands (14, 25, and 37) and CCR5, membrane molecular dynamic (MD) simulations were performed. The position and orientation of ligands (14, 25, and 37) were found to be changed after MD simulations, which demonstrated the ability of this technique to identify binding modes. Furthermore, at the end of simulation, it was found that residues identified by docking were changed and some new residues were introduced in the proximity of ligands. Our results are in line with the majority of previous mutational reports. These results show that hydrophobicity is the determining factor of CCR5 antagonism. In addition, salt bridging and hydrogen bond contacts between ligands (14, 25, and 37) and CCR5 are also crucial for inhibitory activity. The residues newly identified by MD simulation are Ser160, Phe166, Ser180, His181, and Trp190, and so far no site-directed mutagenesis studies have been reported. To determine the contributions made by these residues, additional mutational studies are suggested. We propose a general binding mode for these derivatives based on the MD simulation results of higher (14), medium (37), and lower (25) potent inhibitors. Interestingly, we found some trend for these inhibitors such as, salt bridge interaction between basic nitrogen of ligand and acidic Glu283 seemed necessary for inhibitory activity. Also, two aromatic pockets (pocket I – TM1-3 and pocket II – TM3-6) were linked by the central polar region in TM7, and the simulated inhibitors show important interactions with the Trp86, Tyr89, Tyr108, Phe112, Ile198, Tyr251, Leu255, and Gln280 and Glu283 residues. These results shed light on the usage of MD simulation to identify more stable, optimal binding modes of the inhibitors.


Introduction
Human immunodeficiency virus type-1 (HIV-1) infection, which eventually leads to acquired immunodeficiency syndrome (AIDS) was first discovered by Barre-Sinoussi et al. (1983). AIDS remains a lethal disease, especially in sub-Saharan Africa and Southeast Asia. The chemokine receptor CCR5 has proven to be an exciting pharmaceutical target in the contexts of HIV-1 and inflammatory diseases. CCR5 plays an integral role in the R5-tropic HIV-1 entry process by serving as a critical coreceptor for viral envelop protein gp120 (Kwong et al., 1998;Wu et al., 1996). In fact, individuals with a homozygous 32-base pair deletion in the gene encoding CCR5 do not express the functional receptor and are resistant to R5-tropic HIV-1 infection (Samson et al., 1996). This fact has inspired a great amount of research in the past decade to identify an antiHIV-1 therapeutic that targets the CCR5-mediated entry mechanism (Barber, 2004;Kazmierski et al., 2003Kazmierski et al., , 2005Meanwell & Kadow, 2003;Shaheen & Collman, 2004). This effort recently resulted in FDA approval for small molecular inhibitor, maraviroc (Selzentry) (Dorr et al., 2005) for the treatment of HIV-1 infection. HIV-1 primarily infects vital cells in the human immune system, such as helper T-cells (specifically CD4 + T-cells), macrophages, and dendritic cells (Coffin et al., 1986;Furtado et al., 1999).
It enters these cells by the adsorption of glycoprotein (Fackler & Peterlin, 2000) on their surfaces via receptors followed by fusion of the viral envelope with the cell membrane and intracellular HIV-1 capsid release (Piot, Bartos, Ghys, Walker, & Schwartlander, 2001). HIV-1 infection leads to low levels of CD4 + T-cells via three main mechanisms: First, by the direct viral killing of infected cells; second; second, by increasing rates of infected cell apoptosis; and third, by the killing of infected CD4 + T-cells by CD8 cytotoxic lymphocytes (Abdelwahab et al., 2003). CCR5 receptor is a member of the G-protein coupled receptor (GPCR) superfamily (Strader, Fong, Tota, Underwood, & Dixon, 1994), and has been identified as a primary coreceptor on CD4 + cells for the entry of macrophage-tropic (M-tropic or R5) HIV-1 isolates (Mackewicz, Barker, & Levy, 1996). Furthermore, the homozygous deletion of 32-base pairs in the gene of CCR5 receptor prevents the expression of functional CCR5 on the cell surface and confers resistance to HIV-1 infections, and in one heterozygous individual with one intact CCR5 allele advance to AIDS was comparatively slow (Liu et al., 1996;Marmor et al., 2001).
However, no structural information regarding the precise CCR5 binding site of any known ligand is available, and we considered that a detailed understanding of the binding modes of existing inhibitors would help design more potent antagonists. A number of studies have been undertaken on CCR5 modeling using computational methods, in which approaches, such as quantitative structure-activity relationship (QSAR) analysis, pharmacophore modeling, virtual screening, molecular docking, molecular dynamics simulation, and membrane simulation, were used to deduce the interaction mechanisms of novel inhibitors. Virtual screening structure-based methods are expected to provide better results than ligand-based methodologies. In a previous study, structure-based methods were used to model the physics of protein-ligand interactions, and a combination of 2D and 3D structure-based screening techniques were used to identify 10 CCR5 binders from a library of 1.6 million compounds (Kellenberger et al., 2007). Virtual screening and QSAR studies based on 1-(3,3-diphenylpropyl)piperidinylamide derivatives have shown that for high affinity binding, chemical and structural requirements can be identified using one dimensional (1D) physicochemical properties, 2D topological descriptors, and 3D properties, such as steric, electrostatic, hydrophobic, and hydrogen bond acceptor/donor fields around a family of aligned molecules (Afantitis et al., 2006;Aher, Agrawal, Bharatam, & Garg, 2007). Furthermore, combined molecular modeling techniques, such as comparative receptor modeling, 3D-QSAR, docking, and virtual screening have been used to identify novel entry inhibitors (Carrieri et al., 2009), and a docking-based 3D-QSAR approach has been utilized (Xu et al., 2004) with protein modeling, MD simulation, automated docking, and 3D-QSAR analysis to investigate the detailed interaction of CCR5 with its antagonists. Recently, 4D-QSAR study was performed for Toxoplasma gondi adenosine kinase inhibitors, which is a common cause of secondary central nervous system infection in immunocompromised people such as AIDS patients (da Cunha, Mancini, & Ramalho, 2012). Previous report by Kothandan, Gadhe, and Cho (2012) on dual antagonists docking into the CCR2 and CCR5 suggested the conserved interaction patterns among chemokine receptors. Recently, the design of CCR5 allosteric antagonist's mutagenesis and modeling approach was successfully used, which shows N-terminal region and ECL2 are important in the interactions (Metz et al., 2012). Computational structure-based drug design approach was carried out for antiHIV-1 using gp120 V3 loop (Andrianov, 2008;Andrianov & Anishchenko;. Some groups used homology modeling, docking, and molecular dynamic simulation approaches to understand the mechanism of inhibitors' interaction with proteins (Guimarães, Oliveira, da Cunha, Ramalho, & França, 2011;Oliveira et al., 2011).
In this paper, we describe CCR5 homology modeling, inhibitor (14, 25, and 37) docking, and a molecular dynamics simulation of ligands-proteins (14-CCR5, 25-CCR5, and 37-CCR5) complexes in a hydrated lipid environment. A CCR5 homology model was developed using the recently reported structure of CXCR4 (PDB ID: 3ODU) (Wu et al., 2010) as a template. Automated docking of the highly active (14), medium active (37), and less active (25) ligands was performed and these complexes were simulated in a lipid environment for 20,000 picoseconds (ps) to check its structural effects in the vicinity of CCR5-ligand. The inhibitors chosen for simulations are from a wide range of activity and cover entire data-set. The aim of present compounds selection is to study some general trends. The strategy used in this paper could be appropriate to understand structure-activity relationship of other derivatives and deepens the knowledge of interaction mechanism between ligand-CCR5.

Materials and methods
Template selection and homology modeling zNo crystal structure has been reported for human CCR5, and thus, to identify a suitable template for comparative modeling; the human CCR5 sequence was retrieved from the UniProt database (accession code: P51681) (http:// www.uniprot.org/uniprot/P51681). The basic local alignment search tool (BLAST) at NCBI (http://blast.ncbi.nlm. nih.gov/) was employed to obtain the known homologous protein structure with the higher identity as a template to model CCR5 (Altschul, Gish, Miller, Myers, & Lipman, 1990). A BLAST search retrieved the recently reported CXCR4 structure (Wu et al., 2010) (PDB ID: 3ODU, .25 nm resolution) as a template with higher sequence identity with the target sequence. Multiple sequence alignment of the target (CCR5) and template structure (3ODU) was done using CLUSTALW (Thompson, Higgins, & Gibson, 1994) (http://www2.ebi.ac.uk/ CLUSTALW) program, and the resulting alignment was checked for correct transmembrane (TM) helices and conserved key residues of GPCR, as suggested by Baldwin , Schertler, and Unger (1997). The actual homology model of CCR5 was established using modeller9v4 software (Eswar et al., 2003) using CXCR4 as a template. The modeller calculates a model composed of nonhydrogen atoms, by comparing the sequence alignment to be modeled with those of known related structures. A 3D model was obtained by optimizing a molecular probability density function (PDF) using a variable target function procedure in Cartesian space employing conjugate gradients and molecular dynamics with simulated annealing. A 100 3D models of CCR5 were generated, and the model with the lowest molecular PDF (Molpdf) score and lowest root mean square deviation (RMSD) value was selected for further computational study. The stereochemical quality of the selected CCR5 model was checked using PROCHECK (Laskowski, MacArthur, Moss, & Thornton, 1993), ERRAT (Colovos & Yeates, 1993), and ProSA (https://prosa.services.came.sbg.ac.at/prosa.php). Previously, numerous groups derived homology model to know the interactions of the inhibitors using docking and virtual screening approaches (Bhargavi, Kalyan Chaitanya, Ramasree, Vasavi, & Murthy, 2010;Sarita, Kiran, Bhargavi, Vasavi, & Uma, 2012).

Automated molecular docking
The 3D structures of the oxamino-piperidino-piperidine amide (OPPA) ligands (14, 25 and 37) (Palani et al., 2001(Palani et al., , 2002 were generated using the SYBYL 8.1 (2008) molecular modeling package. The OPPAs other derivatives are reported in Table 1. The resultant 3D structure of the ligand was saved in PDB format for docking analysis.
Autodock 4.0 (Morris et al., 1998) combined with the Lamarckian genetic algorithm, which is used for the search ofglobally optimized conformation, was used for the docking study. Autodock is well known for its ability to produce the lowest energy conformer of a ligand in protein. The grid spacing was .0375 nm in each dimension, and each grid map consisted of 60 Â 60 Â 60 grid points. With prior knowledge of Glu283 as a crucial residue, the center of the grid was set to the position of the C α atoms of Glu283 in a CCR5 cavity. Global optimization was started with a population of 100 randomly positioned individuals, with a maximum of 2.5 Â 10 6 energy evaluations and 2.7 Â 10 5 generations. During the docking experiment, 100 poses were generated. Other parameters were left at their default values. At the end of the docking experiment, cluster analysis was carried out to obtain the global minimum conformation of ligands (14, 25, and 37) in the putative CCR5 cavity. Docking solutions of ligand all-atoms RMSD within .1 nm of each other were clustered together and ranked by the lowest docking energy. Representative binding mode was selected based on the lowest energy cluster of the docking experiment and with the previous knowledge of important residues determined by mutational study.

CCR5 and lipid bilayers assembly
The simulated systems consisted of 14-CCR5, 25-CCR5, and 37-CCR5 complexes embedded in equilibrated 128 dimysristoylphosphatidylcholine (DPPC) lipid bilayers. TM helices were inserted in the preequilibrated DPPC bilayers and the optimal position of TM helices relative to the membrane was determined.
Initially, TM domains were placed in the membrane by overlapping the membrane and the protein, and then all lipids within a cutoff distance of .12 nm from protein atoms were deleted, as described by other authors (Kandt, Ash, & Peter Tieleman, 2007). The Inflategro method developed by Kandt et al. (http://moose.bio.ucalgary.ca) was used to pack lipid bilayers around CCR5 using alternate compression and minimization steps. It was observed that the lipid molecules which were overlapping with the CCR5 and could not be minimized. Then the bilayer was expanded using a scaling factor (4). During these steps, strong position restrains were applied to CCR5 and ligands (14,25,and 37) to maintain their position in the XYZ coordinate system. Expansion of the lipids system was done using a scaling factor of 4 and a distance cutoff of 1.4 nm to estimate areas using a lipid grid size of .5 nm. System compression was done using a scaling factor of .95, followed by energy minimization with strong position restraint. The resultant system was then hydrated in a preequilibrated box of water molecules. To avoid the insertion of water molecules into the protein and lipids, the script at the website of Kandt et al. (2007) was used (http://moose. bio.ucalgary.ca). The van der Waals radius of carbon was artificially changed from .015 to .0375 nm, which made the addition of water molecules to lipid less likely. The central cavity formed by the TM helices was filled with water molecules. The water molecules misplaced in the center of the bilayers (formed by the highly hydrophobic DPPC tails) were removed by visual inspection. After solvation, the van der Waals radius of carbon was reset to .015 nm. The final simulation box dimensions were 6.74 Â 6.76 Â 6.12 nm.
General setup for molecular dynamics simulation All MD simulations were performed using the GROMACS4 package (Hess, Kutzner, Van Der Spoel, & Lindahl, 2008;Van Der Spoel et al., 2005) and the 53A6 GROMOS96 force field (Oostenbrink, Villa, Mark, & Van Gunsteren, 2004;Scott et al., 1999). The DPPC lipid molecule parameters were taken from the 53A6 GROMOS96 force field. The single point charge model (Hermans, Berendsen, Van Gunsteren, & Postma, 1984) was used for water molecules. Simulations were performed in the isobaric and isothermal ensembles. The temperature of the system was kept around 323 K, by coupling to a Berendsen heat bath (Berendsen, Postma, Van Gunsteren, DiNola, & Haak, 1984) with a coupling constant of .1 ps, with separate coupling of solutes (protein, ligand, and lipids) and solvents (water and ions). The chosen temperature for the simulations was well above the phase transition temperature of DPPC (Tm = 315 K) to ensure that bilayers were in the liquidcrystalline state (Nagle, 1993). The pressure was coupled semiisotropically (coupling constant 1.0 ps and compressibility 4.5 Â 10 À5 bar À1 ), which resulted in the independent coupling of lateral P(x + y) and perpendicular (Pz) pressures. The SETTLE algorithm (Miyamoto & Kollman, 1992) was used to keep water molecule bond lengths at their equilibrium values, and the LINCS algorithm (Hess, 2008) was employed to keep all remaining bonds constrained. The time step for integrating the equations of motion was .002 ps. The particle mesh Ewald (Darden, York, & Pedersen, 1993) method with a cutoff value of 1.2 nm was used to calculate long-range electrostatics.
Periodic boundary conditions were applied in all three spatial directions. Ligand topologies and coordinates were generated using the PRODRG (http://davapc1.bioch.dundee.ac.uk/prodrg/) server (Aalten et al., 1996). In the literature, the topologies for ligands generated from PRODDRG server have been questioned for their correct charge groups and charge assignments (Lemkul, Allen, & Bevan, 2010). To overcome the issue raised in literature about PRODRG, we cross-checked the topologies of all the simulated ligands for the correct charge assignment and charge group. The constant-temperature, constantvolume ensemble (NVT) was adopted at a constant temperature of 323 K, a coupling constant of .1 ps, and duration of 100 ps. Generally, a short NVT simulation is followed by the longer temperature stabilization, the constant-temperature, constant-pressure ensemble (NPT) simulation. Because, raising temperature of a system to a desired level (323 K) using NVT simulation is necessary to take care of the introduced voids that are generated while deleting overlapping lipids. Also NPT simulation should not be carried out initially, because pressure calculation at lower temperature could be inaccurate. However, NPT simulation should be run immediately after NVT to correct the density that is too low due to creation of voids and spaces around the solvent while building system. Thus, a short NVT simulation was carried out to stabilize system for temperature and later on switched to longer NPT simulation to stabilize pressure and equilibrate voids. After NPT was performed at a constant pressure of 1 bar, with a coupling constant of 5.0 ps, and duration of 1000 ps. Final MD simulation for 20,000 ps was carried out using a Nose´-Hoover thermostat (Hoover, 1985) at 323 K and a Parrinello-Rahman barostat at 1 bar (Parrinello & Rahman, 1981).

Sequence and homology model analyses
Sequence analysis showed that the BLAST search retrieved a template structure for CXCR4 (PDB ID: 3ODU) with a sequence identity of 35%. This higher sequence identity showed that CXCR4 was a better template for CCR5 modeling than the previously reported bovine rhodopsin and β 2 adrenergic receptor structures. Expectation (E) value represents the number of different alignments with scores equivalent to or better than scores expected to occur during a random database search. Generally, a lower E-value in BLAST indicates that the alignment is real and does not occur by chance. The E-value for CCR5 was 1e-33. To model the 3D structure of human CCR5, sequence alignment was done with the CXCR4 sequence using CLUSTALW (Figure 1). Alignment obtained by CLUSTALW was checked for correct TM prediction with that of the template sequence. The Modeller program was used to generate 100 3D models of CCR5 with the given alignment, and the developed CCR5 models were checked visually for TMs prediction, and the correct TMs were found to be predicted according to the template (PDB ID: 3ODU) structure.
Developed 3D models of CCR5 were checked for their stereochemical quality, nonbonded atomic interaction, and overall energy, and the results obtained showed the stability of the developed model. A comparison between the high-resolution template structure and the CCR5 model was performed to check for normal and/or unusual distributions of u-ψ torsion angles, nonbonded atomic interactions, and energy profiles. In the CCR5 model, the conserved disulfide bridges between Cys20-Cys269 and Cys101-Cys178 were maintained. A final CCR5 model was selected based on a low MolPdf score and RMSD with that of the template. The selected CCR5 model showed .0179 nm RMSD to the template structure (3ODU) backbone. All hydrogen atoms were added to the CCR5 model using SYBYL (2008). ASN and GLN side chains were reoriented using the amide flip option in SYBYL to ensure hydrogen bonding optimization in CCR5. PROCHECK (Laskowski et al., 1993) analysis ( Figure 2) showed that the selected model had 92.8% of residues in most favored regions, 6.8% residues in additionally allowed regions, .4% residues in generously allowed regions, and no residue in a disallowed region. ERRAT (Colovos & Yeates, 1993) plot analysis (Supplementary Figure S1) showed an overall model quality of 95.32%, which indicated that the model was well within the template structure range. Another structure validation parameter, namely, goodness factor (G-factor) was explored. For a good structure, the overall value of G-factor should be above À.5 and close to zero. Our developed CCR5 model had a G-factor value of zero, indicating that the developed model was of good quality. In addition, the developed model was validated for knowledge-based energy using ProSA (Wiederstein & Sippl, 2007) Figure S2). The developed CCR5 model showed an overall Z-score of À3.03, which is well within the range of the native structure. Based on the analyses of CCR5 model (Table 2), we concluded that the developed model was suitable for docking studies and MD simulation.

Molecular docking analysis
Molecular docking is the most useful means of analyzing ligand-protein interactions at the molecular level. For inhibitor docking, the putative CCR5 binding pocket was defined based on previous experimental results (Castonguay et al., 2003;Dragic et al., 2000;Govaerts et al., 2003;Nishikawa et al., 2005;Tsamis et al., 2003).  The red color indicates residues in the most favored region, yellow indicates residues in an additionally allowed region, faint yellow indicates residues in a generously allowed region, and white indicates the disallowed region. In our derived 3D-CCR5 model, no residue was in the disallowed region. According to mutational data, it was apparent that the binding pocket is located in TM1 (Tyr37), TM2 (Trp86, Tyr89), TM3 (Tyr108, Tyr109, Phe113), TM5 (Ile198), TM6 (Tyr251), and TM7 (Gln280, Glu283) with partial involvement of extracellular loop 2 (ECL2). A 100 docked solutions were generated for the each ligand (14, 25, and 37) in the CCR5 active site using the given Autodock protocol. These poses were separated into clusters based on a RMSD tolerance of .1 nm. The representative binding modes for 14, 25 and 37 in CCR5 were chosen from the first cluster which has the more negative binding free energy. In addition, we used the knowledge of previous report (Castonguay et al., 2003;Dragic et al., 2000;Govaerts et al., 2003;Nishikawa et al., 2005;Tsamis et al., 2003) on important residues (Tyr37, Trp86, Tyr89, Tyr108, Phe109, Phe112, Ile198, Tyr251, Asn252, Leu255, Gln280, and Glu283). As well, the most important salt bridge interaction between basic nitrogen of ligand and acidic Glu283 of CCR5 was used to select a candidate pose. As far as previous mutational report is concerned, there was no contradiction and we felt the chosen docked poses could be acceptable. However, it should be noted that the bioactive conformation is not necessarily identical with the lowest energy conformation. On the other hand, it cannot be a conformation that is so high in energy that it is excluded from the population of conformations in solution (Holtje & Folkers, 2003). Thus, calculated information on the bioactive conformation could serve to verify simulated structures wherever possible. Hydrophobic interactions occurred at Trp86, Tyr89, Leu104, Tyr108, Phe109, Phe112, Phe113, Ile198, Leu203, Trp248, Tyr251, and Leu255. The ligand basic nitrogen makes a salt bridge contact with the acidic Glu283 residue of CCR5 at a distance of .42 nm. The pyridine ring of the ligand interacts through π-π stacking interaction with the phenyl ring of Trp86, and docked into a highly hydrophobic pocket formed by residues Phe85, Tyr89, Trp94, and Leu104. Residues from ECL2 (Cys178, Ser179) also participate in binding pocket formation for the ligand and closely interact with the ligand. The piperidine ring, which contains basic nitrogen, docked in a pocket was composed of residues Tyr108, Tyr251, and Glu283. In addition, ethoxime containing substituents of the ligand are situated deep in the binding site of CCR5, which contains Tyr108, Phe109, Phe112, Phe113, and Ile198 residues. The hydrophobic bromophenyl ligand entity docks into the cavity formed by residues Phe112, Ile198, Trp248, Tyr251, Asn252, and Leu255.
Binding pose of the 37 (Figure 3(B)) was analyzed, which indicates that ligand interacts mainly through the hydrophobic interactions. Basic nitrogen of 37 makes salt bridge contact with the acidic Glu283 at a distance of .38 nm. However, hydrophobic contacts were observed at Phe85, Trp86, Tyr89, Trp94, Leu104, Tyr108, Phe109, Ile198, Val199, Leu203, Tyr251, and Leu255. Pyridine ring of 37 oriented similarly as 14 and interacts through π-π stacking interaction with the phenyl ring of Trp86. At this hydrophobic pocket other interacting residues are  (14) are shown as faint brown color sticks. 14 is shown as a cyan colored stick. Salt bridge contact between 14-CCR5 is shown as a dotted blue line. (B) Binding mode of the medium potent inhibitor (37) in the CCR5 model after docking simulations. A top view of the putative binding pocket of CCR5 is shown along with the seven TM helices. Interacting residue side chains within 4 Å of the ligand (37) are shown as faint cyan color sticks. 37 is shown as a violet colored stick. Salt bridge contact between 37-CCR5 is shown as a dotted cyan line. CCR5 is shown as transparent surfaces.Note: The figure was generated using the Pymol program (http://www.pymol.org).
Phe86, Tyr89, Trp94, and Leu104, and show tight interactions. On the other hand, hydrophobic bromophenyl moiety docked into a hydrophobic pocket which was lined by the residues Ile198, Val199, Leu203, Tyr251, and Leu255. ECL2 (Cys178 and Ser179) also participated in the binding pocket formation for 37 and interacted. In addition, ethoxime was located into a pocket formed by the Tyr109, Phe109, Gln194, Thr195, and Ile198.
Similarly, docked mode of 25 ( Figure 4) was analyzed. It shows that ligand mainly interacts through hydrophobic contacts. Salt bridge contact between basic nitrogen of 25 and Glu283 of CCR5 was detected at a distance of .4 nm. However, bromophenyl ring docks into a highly hydrophobic pocket formed by Tyr108, Phe109, Phe112, Phe113, and Ile198. Whereas, dimethylphenyl ring at the other end also docked in similar fashion as 14 and 37. Aromatic π-π stacking interactions were observed between dimethylphenyl and the phenyl ring of Trp86. Dimethylphenyl ring was surrounded by the hydrophobic Tyr89, Trp94, and Leu104. It was also scrutinized that the residues belong to TM1 and TM4 and do not participate in interactions.

Molecular dynamics simulation analysis
Ligands-proteins (14-CCR5, 25-CCR5, and 37-CCR5) complexes were inserted into lipid bilayers of DPPC containing 128 molecules. The longer axes of 14-CCR5, 25-CCR5, and 37-CCR5 were placed on and parallel to the z axis of DPPC. The Inflategro program was used to pack lipid molecules around the protein structure. It removed two lipid molecules from the upper layer and four from the lower layer; in total six lipid molecules were removed and final system consisted of 122 lipid molecules. Alternate compression followed by minimization of system was performed to remove steric clashes between protein and lipid interface. During these steps, ligands-proteins (14-CCR5, 25-CCR5, and 37-CCR5) complexes were kept constrained by applying strong positional restrains to avoid movement. At the final step, the area occupied per lipid molecule was measured, and found to be close to the experimental value of DPPC lipid. After lipids had been packed around the ligandsproteins (14-CCR5, 25-CCR5, and 37-CCR5), the complex systems were solvated with water and 14 chloride atoms each to produce an electro-neutral systems. The final systems used for simulations consisted of 14-CCR5, 25-CCR5, and 37-CCR5 in 122 DPPC bilayers, 14 chlorides, and 2846 water molecules ( Figure 5). A 100 ps of NVT followed by a 1000 ps of NPT simulation were carried out to stabilize the systems for temperature and pressure. The temperature stabilization plot during NVT simulation and the pressure stabilization plot during NPT simulation are shown in Supplementary Figure S3. A 20,000 ps MD simulations were carried out on the temperature and pressure stabilized systems.
The backbone RMSD of the protein (14-CCR5) in lipids as compared with the initial structure immersed in DPPC bilayers was monitored throughout the 20,000 ps simulation, and presented in Figure 6(A). It was observed that, during first 2000 ps of simulation, the RMSD of the protein continuously increased and that subsequently it fluctuated between .35 and .41 nm, with an average RMSD of .38 ± .03 nm. It was observed that the protein structure was affected by its lipid environment, as indicated in Figure 6(A). This plot also indicated that protein structure equilibrated at 2000 ps and then remained stable throughout the 20,000 ps simulation. However, in case of CCR5 in 25-CCR5 system, the backbone RMSD increased constantly up to first 5000 ps and leveled off for the rest of simulation period. After 5000 ps of equilibration, RMSD constantly fluctuated in between .33 and .39 nm, with an average of .36 ± .03 nm. The backbone RMSD for the CCR5 in  37-CCR5 system was monitored. It was detected that the first 5000 ps CCR5 equilibrated and afterward fluctuated in .35-.41 nm frequency, with an average of .38 ± .03 nm.
Protein structure stability during simulation (14-CCR5) was monitored using a potential energy graph ( Figure 6(B)) as a function of time. This graph indicated that all structures lied within the À2.83e+05 to À2.90e +05 kJ/mol range. All structures showed small deviations in potential energy and during the last 5000 of the 20,000 ps their average potential, the energy was À2.89e +05 kJ/mol. Furthermore, potential energy plots showed that during the last 5000 ps structures were energetically stable. The potential energy for the 25-CCR5 (Figure 6 (B)) system was calculated, which showed that the all structures positioned within the À2.77e+05 to À2.85e +05 kJ/mol energy range. It was also observed that the potential energy decreased gradually. The average potential energy of last 5000 ps of 20,000 ps was À2.83e +05 kJ/mol. Similarly, potential energy for the 37-CCR5 system was calculated which indicates that the gradual decrease in potential energy throughout simulation period. The average energy of last 5000 ps of 20,000 ps was calculated and it was À2.83e+05 kJ/mol. All the structures had energies between À2.77e+05 and À2.85e +05 kJ/mol.
The temperature during 20,000 ps simulation was monitored and the results obtained confirmed that throughout the 20,000 ps simulation systems (14-CCR5, 25-CCR5, and 37-CCR5), temperatures were maintained at 323 K (Figure 7(A)), indicating that the system was simulated under constant temperature conditions. In addition, systems' (14-CCR5, 25-CCR5, and 37-CCR5) pressures were monitored through the 20,000 ps, and found to be equal to the 1 bar (Supplementary Figure S4). The RMSDs of 14, 25, and 37 atoms were calculated as a function of time during 20,000 ps simulation and results are plotted in Figure 7(B). It shows that the RMSD of ligand (14) atoms remained constant ($.065 nm) up to 14,500 ps, then dropped to $.05 nm, and after 16,000 ps rose and leveled off at $.065 nm throughout the remainder of the simulation. The RMSD for 25 (Figure 7(B)) was calculated which shows that up to 14,192 ps the average RMSD was $.065 nm and afterward it elevated to $.116 nm for rest of the simulation period. However, for 37 another trend of RMSD was noted. The RMSD constantly rose up to 5000 ps and afterward gradually decreased throughout 20,000 ps simulation. The average RMSD for 37 was found to be .137 nm. After MD simulations, the complexes (14-CCR5, 25-CCR5, and 37-CCR5) were extracted from trajectories and minimized for 100 steps using the steepest descent method with the aim of removing bad contacts from the crude complex structures. The steepest descent minimized structures were then reevaluated for u-ψ torsion angles (PROCHECK plot) (14-CCR5,  Figure S8). The PROCHECK plot (Figure 8) showed that 99.2% of residues lay in the most favored and additionally allowed regions, whereas only .8% of residues were in the generously allowed region, and no residue was in a disallowed region, thus indicating reasonable model stability after MD simulation. However, PROCHECK plot for 25-CCR5 (Supplementary Figure S5) shows that the 98.5% residues are in the most favored and additionally allowed region and 1.2% residues are in generously allowed and .4% in disallowed region. PROCHECK plot for 37-CCR5 (Supplementary Figure S6) shows 98.1% residues are in most favored and additionally allowed region, 1.5% residues are in generously allowed, and .4% in disallowed region. ERRAT plot (Supplementary Figure S1, 7) analyses indicated an overall model quality of CCR5 after MDS in 14-CCR5 was 92%, in 25-CCR5 was 95.24%, and in 37-CCR5 was 88.35%. However, the ProSA energy plots (Supplementary Figure S2, 8) showed a reliable Z-score for all the models (14-CCR5, À2.37; 25-CCR5, À2.88; 37-CCR5, À2.30) indicating that structures were well within the Z-score of template structure (À2.34) after MD simulations. Parameters evaluated for the extracted CCR5 structures were summarized in Table 2.
Radii of gyration (Rg) for CCR5 were calculated and plotted (Supplementary Figure S9). The Rg of protein in simulated system is a measure of the compactness of protein aggregation. Throughout the CCR5 simulation (14-CCR5), a slow increase in Rg from $2.075 to $2.185 nm was observed, which indicates CCR5 compactness slowly diminished. However, in 25-CCR5, the Rg value was found to be increased from $2.05 to $2.12 nm up to first $4000 ps and afterward it maintained between $2.07 and $2.12 nm, which indicated that CCR5 structure was compact throughout the simulation and supported that obtained structure is of good quality (Supplementary Figure S7A). In 37-CCR5 system, it was observed that the Rg was increased up to first 5000 ps from $2.062 to 2.13 nm, afterward up to 15,000 ps simulation it was maintained and again rose up in last 5000 ps. It indicates that the compactness of CCR5 was maintained between 5000 and 15,000 ps and afterward it slowly diminished the CCR5 packing.

Structural fluctuations
Atomic fluctuations of CCR5 (14-CCR5, 25-CCR5, and 37-CCR5) residues (also called root mean square fluctuations (RMSF)) were examined (Figure 9). RMSFs for protein backbone (C, N and C α ) atoms were generated and plotted (Supplementary Figure S10). The RMSF plot for protein (14-CCR5) residues (Figure 9) showed that the fluctuations were greatest for Lys138 of the intracellular loop (ICL) 2 at $.43 nm, Gln186 of the extracellular loop (ECL) 2 at $.44 nm, Lys228 (ICL3) at $.44 nm, His231 Figure 8. PROCHECK plot of the u-ψ distribution of the CCR5 model after MD simulation. After MD simulation, 89.1% of residues were in the most favored region (A, B, and L), 10.1% were in the additionally allowed region (a, b, l, and p), .8% were in the generously allowed region ($a, $b, $l, and $p), and no residue was in the disallowed region. (ICL3) at $.42 nm, and Phe304 (C-terminal) $.43 nm. However, the fluctuations profile of 25-CCR5 system shows that it reaches maximum of $.32 nm with the average fluctuations which were found to be below .2 nm. It is also observed that the very few fluctuations exceeded .3 nm such as Ile23 (.31 nm), Ile27 (.31 nm), and Leu174 (.32 nm). For the 37-CCR5 system, it was scrutinized that most of the residues fluctuated less than .2 nm and very few residues fluctuated more than .3 nm (Arg60, .31 nm; Tyr184, .43 nm; Arg225, .33 nm; Arg232, .33 nm; Asn267, .34 nm; and Glu302, .38 nm). From Figure 9 it is also evident that the average fluctuations were below .2 nm. A few fluctuations exceeded .3 nm, indicating that these atoms belong to residues in the N-terminal region, ECL1, ECL2, ECL3, ICL1, ICL2, and ICL3. RMSFs for protein backbone (14-CCR5) atoms ( Supplementary Figure S10) showed that fluctuation was greatest (.4 nm) for residue Thr167 of ECL2. Few residues (Ile23, Ile27, Arg60, Leu174, Tyr184, Arg225, Arg232, Asn267, and Glu302) fluctuated more than .3 nm, and most showed fluctuations below .2 nm. These results also indicate that the N-terminal, C-terminal region, ECL2, ICL2, and ICL3 of CCR5 fluctuate considerably, whereas the β sheet and α helical regions fluctuate by less than .2 nm, which indicates these regions are stabilized by intramolecular hydrogen bonding and disulfide bridges. During MD simulation, it was observed that Tyr108, Tyr251, and Glu283 form intramolecular hydrogen bonds and stabilize the TMs architecture. This indicates that β sheet and α helices are more stable than the N-terminal, C-terminal, ECL2, ICL2, and ICL3 regions of CCR5. The average CCR5 structure was calculated from MD trajectories and the b-factors, as shown in Supplementary Figures S11-S13. Red, orange, and yellow color in S10-S12 was identified in N-terminal, C-terminal, ECL1-3, and ICL1-3, which indicated mobile region in MD simulation.
It can also be seen that ligand-protein (14-CCR5, 25-CCR5, and 37-CCR5) complexes embedded in lipid bilayers stay in equilibrium throughout MD simulations, and hence, provide a more relaxed structures (Figure 8, Supplementary Figure S5, 6) for further analyses. Furthermore, it is evident from Figure 8 that after MD simulation no residue was located in a disallowed region, and in addition, ERRAT plot analysis showed a good overall score for CCR5 after MD simulation (Supplementary Figure S1). After careful examination of the Ramachandran plot (Supplementary Figure S5, 6), it was observed that (25-CCR5 and 37-CCR5) after MDS, few residues were found to be in disallowed regions (Lys303 and Ala133), but they are not a part of active site region.

The characterization of loop regions
In CCR5, there are six loop regions, namely, ECL1, ECL2, ECL3, ICL1, ICL2, and ICL3 are present. Nterminal and C-terminal uncharacterized segments whose geometric constraints are also not well known. It is also known that two disulfide bridges are present in CCR5 that is, Cys20-Cys269 and Cys101-Cys178. The disulfide bridges connecting ECLs in chemokine receptors are important for ligand binding (Genoud, Kajumo, Guo, Thompson, & Dragic, 1999). In the developed model, we generated disulfide bridges between the above-mentioned cysteine residues. However, previous study by Kaushik, Mohanty, and Surolia (2012) reported the importance of disulfide bridges in the enhanced structural stability and function. In the 20,000 ps MD simulations, it was observed that disulfide bridges were intact and that they maintained the conformational integrity of CCR5. Figure 10 represents the RMSDs of the two disulfide bridges in differently simulated systems (14-CCR5, 25-CCR5, and 37-CCR5). For Cys20-Cys269 ( Figure 10(A)), the RMSD between disulfide bonds was variable and the highest value reached was $.052 nm, whereas for Cys101-Cys178 (Figure 10(B)) the RMSD was averaged $.015 nm, which indicates that this disulfide bridge maintains the conformational integrity of CCR5. It is worthwhile noting that the disulfide bridges help maintain the structural organization of CCR5, which is in line with a previous report (Blanpain et al., 1999;Genoud et al., 1999). The average RMSD for Cys20-Cys269 in 25-CCR5 was $.015 nm and in 37-CCR5 it was $.02 nm. However, the RMSD for Cys101-Cys178 in 25-CCR5 was found to be $.012 nm and in 37-CCR5 it was $.011 nm.
To calculate backbone structural deviation during MD simulation as a function of time, N-terminal domain (Pro19-Ala30), ECL1 (Ala90-Gln102), ECL2 (Thr167-Ile198), and ECL3 (Gln161-Gln277) were treated as different groups (RMSDs for all are provided in Figure 11). However, it was found that the ECL2 loop (14-CCR5) was more dynamic and fluctuated with maximum amplitude of $0-.5 nm. The average RMSD for ECL1 was $.082 ± .002 nm and for ECL3 was $0-.135 ± .002 nm. In case of 25-CCR5 system, it was observed that the ECL1 has RMSD of $0-.225 nm, whereas in case of ECL2 RMSD reaches from $0 to À.4 nm, and in ECL3 it was $0-.21 nm. It should be noted that, in 37-CCR5 system, the ECL2 loop RMSD was $0-.3 nm and for ECL3 loop it was $0-.275 nm with the average of $.225 nm. However, lowest RMSD was observed for ECL1 in 37-CCR5 system and it was $0-.16 nm. Furthermore, it was observed that ECL2 was the most mobile loop among ECL 1-3. The RMSD of the Nterminal segment was found to be within the range $0-.35 nm. In 25-CCR5 and 37-CCR5 systems, more deviations can be seen for N-terminal RMSDs. Our results indicate that the N-terminal segment and ECL2 are the most mobile regions of CCR5, which is in line with the previous reports (Blanpain et al., 1999;Genoud et al., 1999;Horuk, 1994;Liu, Shi, Liu, & Sun, 2004). The RMSD of backbone atoms of ICL1 (Lys59-Tyr68), ICL2 (Asp125-Thr141), and ICL3 (Lys219-Aeg235) are shown in Supplementary Figure S14.

The 14-CCR5 interaction
After 20,000 ps of MD simulation, 14-CCR5 complex was extracted and minimized by 100 steps using the steepest descent method to remove bad contacts to analyze the interactions. Figure 12 with few changes in interacting residues. In particular, Figure 12(A) shows that ligand reorientation occurred in the CCR5 cavity with minimal changes in the orientations of key interacting residues. Comparison of ligand docking modes in CCR5 (Figure 3(A)) with the extracted complex ( Figure 12(A)) after MD simulation, indicates that 14 was reoriented in CCR5 with changes in the positions of pyridine, ethoxime, and bromophenyl substituents. The 2,6-Dimethyl-3-pyridinyl group of 14 was reoriented into the hydrophobic cavity formed by residues Tyr37, Trp86, Ala90, and Tyr89. Although, docking predicted π-π interaction between this ring and Trp86 disappeared during the course of the MD simulation, after simulation the ring interacted edge-to-edge with Tyr37 and Tyr89. In docked mode, Tyr37 was far from the 2,6-dimethyl-3-pyridinyl group, but after MD simulation the 2,6-dimethyl-3-pyridinyl group appeared to have moved toward Tyr37, which indicates that the importance of Tyr37 for the interaction between 14 and CCR5. Furthermore, it was also found that Ser180 and His181 from ECL2 closely interacted with the ethoxime part of 14. However, this ethoxime was relocated in CCR5 from Tyr108, Phe109, Phe112, Phe113, and Ile198 to an adjacent pocket containing Thr105, Phe109, Ser180, His181, Phe166, Trp190, and Gln194. In addition, close interactions were found between the ethoxime region of 14 and Ser180, His181, Phe166, and Gln194.
The bromophenyl substituent of 14 was moved to the binding pocket formed by residues Tyr108, Phe109, Phe112, Ile198, Tyr251, and Leu255 from its original position in a pocket formed by Phe112, Ile198, Trp248, Tyr251, Asn252, and Leu255. Bromophenyl ring exhibited a π-π stacking interaction with Phe109, and face-to-edge π-π interaction between bromophenyl and Tyr251 was observed. After MD simulation, the bromophenyl part of 14 was stabilized in a hydrophobic part of CCR5. Furthermore, it was found that the bromine group interacted Phe112 via a C-Br-π interaction (distance .36 nm). Salt bridge interaction was also found to be important between the basic nitrogen of 14 and the acidic Glu283 of CCR5. In 14-CCR5 model, after docking this distance was .42 nm, but after MD simulation, 14 moved closer to Glu283 to a salt bridge distance of .38 nm. This observation indicates that this salt bridge interaction plays crucial role in the interaction between 14 and CCR5. After MD simulation, it was observed that 14 was stabilized in CCR5 by a hydrogen bond interaction, which was not detected in docked mode, that is, Gln280 of CCR5 forms a hydrogen bond (distance .17 nm) with the carbonyl group of 14. The angle criteria were set to > 90.0°to select hydrogen bond between A-H-D (Torshin, Weber, & Harrison 2002). The computed hydrogen bond distance between 14 and CCR5 indicates that the ligand is bound tightly in the CCR5 cavity.

The 37-CCR5 interactions
After 20,000 ps of MD simulation, 37-CCR5 complex was extracted and minimized by 100 steps using the steepest descent method to remove bad contacts to analyze the interactions. Figure 12(B) depicts the binding mode of 37 in CCR5 after MD simulation; surrounding .4 nm residues were shown for clarity; 37 interact mainly through hydrophobic interactions. Salt bridge contact between basic nitrogen of 37 and acidic Glu283 was detected at a distance of .36 nm. The bromophenyl ring of 37 shows similar orientation as of 14 and interacts with the residue mostly by hydrophobic contacts. Bromophenyl ring positioned into a pocket lined by the residues Phe109, Phe112, Ile198, Trp248, Tyr251, and Leu255. On the other hand, the hydrophobic ring was docked into a pocket formed by TM2, TM3, TM7, and ECL2, which was lined by the Trp86, Tyr89, Trp94, Leu104, Cys178, and Tyr108 residues. Partial involvement of hydrophilic residue (Gln280) was observed in this hydrophobic pocket. Ethoxime of 37 was located into the pocket formed by the Thr105, Phe109, Ser160, Trp190, and Gln194. It is worthwhile to note that the internal hydrogen bonding was observed between Glu283, Tyr108, and Tyr251 which maintained the structural organization of CCR5.
The 25-CCR5 interactions MD simulated binding mode of the 25 (lower active compound) in CCR5 is shown in Figure 13; 25 mainly interacted through the hydrophobic interactions. Residues surrounding .4 nm to 25 in the binding pocket were shown for clarity. Salt bridge contact was observed between acidic Glu283 and basic nitrogen of 25 at a distance of .38 nm. No hydrogen bond contact was detected. It was observed that CCR5 was tightly packed in TMs region and internal hydrogen bonding was detected between Glu283 and Tyr108, Tyr251 and Gln280, which is evident through the constant Rg (Supplementary Figure S9). The 2,6-dimethylphenyl of 25 was stabilized hydrophobically in pocket formed by Trp86, Tyr89, Trp94, and Leu104. Trp86 shows π-π stacking interaction with the 2,6-dimethylphenyl and holds it tightly in the binding pocket. Face-to-edge π-π interactions were detected between Tyr89, Trp94, and 2,6-dimethylphenyl. On the other hand, bromophenyl substituents of 25 were stabilized into the hydrophobic pocket lined by the residues Tyr108, Phe112, Ile198, Val199, Leu203, Tyr251, Asn252, and Leu255. Bromophenyl shows edge-to-edge hydrophobic interaction with Phe112 and face-to-edge interaction with Tyr251. The bipiperidine of 25 was placed into the central polar region of the CCR5 lined by the residues Tyr108, Gln280, and Glu283. Polar hydroxyl group of 25 is found to be interacting with the Thr195 and Ile198. A 2D schematic plot Ligplot (Wallace, Laskowski, & Thornton, 1995) of the ligandsproteins (14-CCR5, 37-CCR5, and 25-CCR5) interactions were depicted in Figures 14 and 15.

SAR study of oximino-piperidino-piperidine-amide (OPPA) derivatives
The binding modes of 14 (IC 50 = .2 nM), 25 (IC 50 = 76 nM), and 37 (IC 50 = 6.4 nM) after 20,000 ps MD simulations (Figures 12 and 13) were considered as representative bioactive conformations and used to elucidate the SARs of other OPPA derivatives, previously reported by Anandan et al. (Palani et al., 2001(Palani et al., , 2002. Exploitation of SARs studies on other OPPA derivatives was performed using the MD simulated conformer of 14. shown for clarity. The 2,4-dimethyl-3-pyridinyl substitution of 14 was found to be necessary to increase binding affinity, because this substitution enables interactions with the hydrophobic residues Ala29, Tyr37, Ala90, and Tyr89 and the base of this pocket was formed by Trp86. However, the replacement of 2,4-dimethyl substitution with heteroatoms, such as NH 2 , OH, and Cl afforded compounds 9 (IC 50 = .6 nM), 10 (IC 50 = .5 nM), and 11 (IC 50 = 4 nM), which showed slight reductions in CCR5 binding potency, indicating that hydrophobic 2,4dimethyl substitution is crucial for inhibitory potency. Moreover, this finding also emphasizes that in the hydrophobic environment surrounding CCR5, hydrophobic substitution is likely to increase binding. Figure 12(B) depicts the interactions of the medium active compound 37 (IC 50 = 6.4 nM) in CCR5 after MD simulation. It shows that the methyl substituent at 6th position of pyridine ring closely interacts with Trp86 and Tyr89 in the pocket, which shows 6th position is sterically forbidden for the improvement of activity.
For compounds 36 (IC 50 = .9 nM) and 37 (IC 50 = 6.4 nM; Figure 12(B)), the 6-hydroxy-2,4-dimethyl-3-pyridinyl and 2-hydroxy-4,6-dimethyl-3-pyridinyl substituents showed reduced binding potency as compared with 14 (IC 50 = .2 nM), indicating the 6th position is sterically restricted and that any substitution at this region could reduce activity, which was evident in our MD model of compound 14 (Figure 12(A)), and that the 6th position of the ligand closely interacts with Ala90. On the other hand, replacement of the ethoxime of 14 with other linkers, such as the sp 3 carbon in compound 16 (IC 50 = 10 nM), the carbonyl in compound 24 (IC 50 = 12 nM), or alcohol in compound 25 (IC 50 = 76 nM; Figure 13) provided no advantage in binding potency. This may be because the ethoxime is a large bulky hydrophobic entity that could appropriately fill the hydrophobic cavity, whereas but others are comparatively small and less bulky. From Figure 13, it can be seen that the alcohol substituent is interacting with Thr195 and Ile198, but it cannot interact with the other residues (Thr105, Phe109, Ser160, Phe166, His181, Trp190, and Gln194) identified by 14 and 37, which clearly indicates that the cavity filling hydrophobic substitutions (like ethoxime or methoxime) are necessary to improve the binding affinity.
In compounds 1 (IC 50 = 10 nM), 2 (IC 50 = 5.5 nM), and 3 (IC 50 = 13 nM) the linker X=CH 2 , NH and CO showed lower potency compounds compared to 14, suggesting the necessity of bulky group at this position for improved inhibitory activity, which is quiet evident through our MD simulated models ( Figures  12 and 13). It also indicates that the TMs 3-6 enclosed a large hydrophobic cavity and the appropriate bulky substituents filling this cavity can lead to a better inhibitor. On the other hand, the large hydrophobic pocket is formed by the TMs 1-3, and 7, indicates the importance of hydrophobic substituent for improved potency.
However, it can be seen from our MD simulated model (Figure 12(A)), the ethoxime is located in a cavity containing residues Phe109, Ser180, Phe166, His181, Trp190, and Gln194 and closely interacts with Phe109, His181, Phe166, and Gln194. Furthermore, the removal of ethyl from the ethoxime of 14 afforded compound 30 (IC 50 = 10 nM), a poorer inhibitor, which suggests that an appropriate bulky substitution at this position is needed to promote affinity (Figure 12(A)). However, removal of ethoxime moiety afforded compounds 17 (IC 50 = 23 nM), 21 (IC 50 = 38 nM), and 22 (IC 50 = 26 nM) are lower active inhibitors and proposes importance of steric bulk to encourage binding affinity. Replacement of the ethyl group of ethoxime of 14 with the more bulky t-butyl and i-propyl groups afforded compounds 31 (IC 50 = 5.7 nM) and 32 (IC 50 = 1.7 nM), respectively, which showed reduced inhibitory potency, indicating that additional bulky substitution confers no advantage. This situation can be rationalized using the MD simulated model of 14 (Figure 12(A)), which shows that these substituents could interact sterically with the binding pocket formed by residues Phe109, Ser180, His181, Phe166, Trp190, and Gln194. MD simulation results of 14 also explain why E-isomer is less active than the Z-isomer of OPPA derivatives. The E-isomer appears to interact well with a steric pocket lined by the residues Thr105 and Ser180 and to show low affinity for CCR5, whereas the Z-isomer extends in a pocket formed by the residues Figure 15. 2D schematic interactions plot between 25 and CCR5 generated by Ligplot server.
The SARs data also showed that a methyl or ethoxime moiety in addition to the 2,6-dimethylnicotinamide N-oxide moiety of 15, 39-44 is optimal for inhibitory potency. The bromine group of compound 15 (IC 50 = .6 nM) was replaced by CF 3 and OCF 3 in compounds 42 (IC 50 = 1.3 nM) and 43 (IC 50 = 1.8 nM), respectively, which displayed reduced inhibitory potencies, indicating that bulky groups in this region diminish inhibitory activity. Figure 12(A) and (B) provides a view of this region in 14-CCR5 and 37-CCR5, and demonstrates that the bromophenyl ring interacts sterically with the pocket wall containing residues Tyr108, Phe109, Phe112, Ile198, Tyr251, and Leu255. At the bottom of this pocket, Br closely interacts with Phe112 and bulky substituents like CF 3 and OCF 3 are difficult to accommodate, which explains the reduced inhibitory effects of these compounds.
Our interpretation of SARs data and MD simulation results of 14, 25, and 37 with CCR5 in lipid bilayers underscore the importances of the hydrophobicity of CCR5 cavity and R 1 , R 2 region of 14, 25, and 37 to design more potent CCR5 inhibitors.
However, we employed a recently reported CXCR4 crystal structure as template. The main advantages of this template structure over others are its higher sequence identity and expanded binding pocket. CCR5 shares 35.4% sequence identity with CXCR4, which is higher than those of bovine rhodopsin (23.3%) and β 2 adrenergic (25.4%) receptors. Based on a report issued by Kimura et al. (2008), it is evident that models based on bovine rhodopsin have smaller binding pocket volume. In this previous study, a balloon expansion method was used to expand the pocket.
The 3D molecular structures of ligand/protein complexes are considered foundations of structure-based drug design. Frequently, 3D data are available for the protein and drug separately, but infrequently for the protein/drug complex. Molecular docking is a process that enables a ligand and its receptor protein to be arranged in their native 3D conformation. After validation of the developed CCR5 model, to check ligand-CCR5 interactions at the molecular level, automated docking of inhibitors (14, 25 and 37) were performed using a CCR5 model using the Autodock program (Morris et al., 1998). Our docking results (Figures 3 and 4) are in accord with those of previous mutational studies, and indicate that the ligand is docked in a pocket containing TM1, TM2, TM3, TM5, TM6, and TM7 (Castonguay et al., 2003;Dragic et al., 2000;Govaerts et al., 2003;Nishikawa et al., 2005;Tsamis et al., 2003). Compound 14 interacted primarily hydrophobically with the CCR5 cavity. Furthermore, the basic nitrogen of 14 formed a salt bridge with acidic Glu283 of CCR5 at a distance .42 nm. For CCR5, it has been shown that the salt bridge interaction between the basic nitrogen of a ligand and acidic Glu283 of CCR5 plays a crucial role, as demonstrated in a previous report, in which it was confirmed that a mutation of Glu283Ala reduced the binding affinity of TAK779 (Castonguay et al., 2003;Seibert et al., 2006). In the docked models of the 37 ( Figure 3B) and 25 (Figure 4), salt bridge contact was detected which shows the importance of it for the current inhibitors for affinity. Both the inhibitors interacted mainly through the hydrophobic interactions.
Molecular dynamic simulations can provide intimate details of individual particles motions as a function of time, and thus, can be used to tackle specific questions about ligand-receptor interactions without resorting to laboratory-based experimentation. Although experiments play a crucial role in the validation of simulation methodologies, as comparisons of simulation and experimental data serve to confirm the validities of calculated results and provide criteria for improving simulation methods. The elucidation of ligand binding mechanisms is the essential step to the identification of more potent and selective inhibitors of CCR5. In the present study, the docked modes of 14, 25, and 37 complexed with CCR5 were allowed to experience its natural environ-ment in presence of DPPC lipid bilayers for as long as 20,000 ps. It was found that the structures of CCR5 simulated with 14, 25, and 37 were affected by the lipid environments, as shown in Figure 6(A). Furthermore, the potential energy graph indicated that the energy of the CCR5 systems decreased gradually with time ( Figure 6 (B)). Figure 7(A) shows that the systems were simulated under constant temperature and supplemental Figure S4 indicates constant pressure was maintained throughout simulations.
After MD simulations, the structures (14-CCR5, 25-CCR5, and 37-CCR5) were extracted and minimized using the 100 step steepest descent method and validated using different methods, that is, by using PROCHECK (Figure 8, Supplementary Figure S5, 6) and ERRAT (Supplementary Figure S1, 7) which showed good correspondence with the initial model (Table 2). Overall analyses of the CCR5 structures after MD simulations indicated the reliability of it. Structural fluctuations over the last 5000 ps of MD simulations were monitored for all residues (Figure 9) and for backbone atoms (Supplementary Figure S10). It was found that Nterminal residues and residues in the ECLs and ICLs region show most movement.
It was observed that the disulfide bridges (Cys20-Cys269 and Cys101-Cys178) created in the CCR5 model were intact after MD simulations in all the simulated systems (14-CCR5, 25-CCR5, and 37-CCR5). Distances between disulfide bridges are shown in Figure 10. The MD simulation of 14-CCR5, 25-CCR5, and 37-CCR5 in membrane rationalizes previous mutational and modeling data, and shows that disulfide bridges between Cys20-Cys269 and Cys101-Cys178 plays a crucial role in maintaining the structural organization (Blanpain et al., 1999;Genoud et al., 1999). However, Kaushik et al. (2012) also explained the importance of disulfide bridges in maintaining structural organization. Loop region characterization showed that the Nterminal region and ECL2 (among ECL 1 to 3) are the most dynamic (Figure 11), which is in line with previous observations (Blanpain et al., 1999;Genoud et al., 1999;Horuk, 1994;Liu et al., 2004) and shows that the dynamic nature of the N-terminal region and of ECL2 of CCR5 are required for viral entry into host cells.
After MD simulations, complex structures (14-CCR5, 25-CCR5, and 37-CCR5) were analyzed for close interactions (Figures 12 and 13), and it was found that ligands (14, 25, and 37) interact mainly through salt bridge and hydrophobic interactions. The 2D interactions between ligand-CCR5 (14-CCR5, 25-CCR5, and 37-CCR5) were interpreted using the Ligplot program (Figures 14 and 15). MD simulation revealed that compound 14 was reoriented in the CCR5 binding pocket after the introduction of new residues, such as Ala29, Ala90, Ser180, His181, Phe166, and Trp190. Whereas, in case of 25 and 37, after MD simulation, it was examined that the docking identified Phe85 is not an important residue for inhibition and vanished. In addition to hydrophobic interactions, it was found that the 14 was secured in CCR5 by hydrogen bonding and a salt bridge. The distance between the basic nitrogen of compound 14 and Glu283 decreased from (before MD) .42 nm to (after MD) .38 nm, indicating the importance of Glu283 in 14-CCR5 interaction, as has been previously suggested (Castonguay et al., 2003;Seibert et al., 2006). The smallest hydrogen bond distance (.17 nm) between the carbonyl of 14 and Gln280 of CCR5 is the evidence of their strong interaction, and pulling of the 2,4-dimethyl-3-pyridinyl ring of 14 towards Tyr37 it indicates the importance of this residue in the interaction between 14 and CCR5, which is consistent with previous mutational results (Seibert et al., 2006).
However, MD simulated models of 25 and 37 interacted mainly through hydrophobic and salt bridge interactions, and no hydrogen bond interactions were detected. The distance between basic nitrogen of ligands (25 and 37) and acidic Glu283 (before MD; .4 and .38 nm) of CCR5 has been decreased after MD simulation (after MD; .38 and .36 nm), which indicates the importance of this interaction for the CCR5 inhibition. It seems that acidic Glu283 of CCR5 pulls ligands (25, 37) basic nitrogen toward it. From this observation, we could say that the salt bridge contact between ligand-CCR5 might play a major role in the inhibition mechanism. MD simulation also predicted the importance of Ala29 and Ala90 in the vicinity of the 2,4-dimethyl-3-pyridinyl moiety. These residues closely interact with this moiety and further substitution at the 6th position of the pyridinyl ring by hydroxyl or methyl group reduced binding activity, as indicated by our SAR study (Palani et al., 2001(Palani et al., , 2002 and the MD simulated 14-CCR5 complex ( Figure 12(A)). Indeed, a previous mutational study (Seibert et al., 2006) also emphasized that bindings to TAK779, AD101, and SCH-C are altered by mutation of Ala29 and Ala90 with aliphatic residues, which also suggest that these residues line the drug binding pocket. Our MD simulated models (25-CCR5 and 37-CCR5) ( Figures  12(B), and 13) also predicted the importance of Trp86 in inhibition mechanism, which is supported by the previous report (Kondru et al., 2008) which stressed that the mutation of Trp86Ala showed significant reduction in binding affinity for Maraviroc, Vicriviroc, TAK220, and TAK779. However, Figures 12(B) and 13 show that the Trp94 is somewhat distant and shows less interaction which is in line with the previous report (Kondru et al., 2008), indicating that mutation of Trp94Ala had no significant effect on CCR5 inhibition. MD simulation results of 25-CCR5 ( Figure 13) point out that the Thr195 (TM5) is an important residue and supported by the previous report of Kondru et al. (2008), which highlighted that the mutation of Thr195Ala reduced the binding affinity for Aplaviroc, Maraviroc, Vicriviroc, and TAK779 in CCR5.
However, the MD simulated model showed that the bromophenyl and ethoxime moiety changed their positions in CCR5 cavity (Figure 11(B)). Moreover, the bromophenyl ring was placed in a pocket lined by highly hydrophobic residues, such as Tyr108, Phe109, Phe112, Ile198, Tyr251, and Leu255. Previous mutational reports have shown that Ile198 is an important residue in the binding pocket, which is in line with our results, and that mutation of Ile198 affects the binding of TAK779 (Billick et al., 2004;Dragic et al., 2000;Nishikawa et al., 2005;Paterlini, 2002). MD simulated models (14-CCR5, 25-CCR5 and 37-CCR5; Figures 12 and 13) identify that I198 is one of the important residue in the strongest hydrophobic interactions of this series of molecules, which is in line with the prior mutational study (Kondru et al., 2008) indicating that mutation of Ile198Ala markedly reduced binding affinity for the Aplaviroc, Maraviroc, Vicriviroc TAK220, and TAK779.
However, Nishikawa et al. (2005) explained that Leu255 is an important residue for the binding of small molecules like TAK220 and TAK779, which is in accord with our results. In addition, they also explained the importance of hydrophobic residues Ile198 and Asn252 for the inhibition of TAK220 and TAK779. They proposed that the TMs 4, 5, and 6 are involved in the formation of drug binding pocket for TAK220 in addition to previously reported TMs 1-3, and 7. Indeed, the observation of Nishikawa et al. (2005) is in line with our MD simulated model emphasizing that the binding pocket for presently simulated compounds were formed by TMs 1-7. However, TM1 and 4 residues do not participate in the binding pocket formation for 37 and 25. This indicates that the most potent compound (14) interacts with all the TMs, whereas medium potent compound (37) does not interact with the TM1, and the least potent compound (25) is devoid of any interactions from TM1 to TM4. Results analysis indicated that the Phe166 from TM4 (Figure 12(A)) participate in close interaction with the ethoxime moiety of 14, which is in line with Nishikawa et al. (2005) results which stressed that mutation of Phe166Ala had little effect on antiviral activity of TAK220 and TAK779. Met279 from TM7 (Figure 12 (B)) and Thr195 from TM5 ( Figure 13) could have minor contributions for CCR5 antagonism with 37 and 25 and are consistent with the previous study (Nishikawa et al., 2005). Our simulated models showed that the salt bridge distance between basic nitrogen of ligands (14, 25 and 37) and acidic Glu283 of CCR5 has been reduced after MD simulation, which indicates the importance of this interaction. In fact, it seems Glu283 anchor the antagonists into CCR5 cavity (Dragic et al., 2000;Tsamis et al., 2003) and pulls ligands toward it. Nishika-wa et al. (2005) also illustrated that the mutation of Glu283Gln abolishes the activity of TAK220 and TAK779, which is consistent with our results, shows significance of Glu283 for CCR5 antagonism. It should be noted that Glu283 on TM7 is shared by all three inhibitors (14, 25 and 37), which is in accord with the existence of potentially positively charged core in the ligands. Furthermore, consistent with our result, the extensive computer modeling of TAK779 in CCR5 predicted the interaction with Asn252 and Leu255 in TM6 (Paterlini, 2002). However, Castonguay et al. (2003) reported that aromatic Tyr251 is necessary to bind small molecular inhibitors (2-aryl-4-(piperidin-1-yl)butanamines and 1,2,4-trisubstituted pyrrolidines) to CCR5, which is consistent with our results. Tyr251 interacts with the bromophenyl of the ligand in an edge-to-face (T-shaped) aromatic stacked form of interaction, which is critical for CCR5 interactions. Phe109 shows π-π stacking with the bromophenyl ring and stabilizes it in the binding pocket. Furthermore, halogen-pi interaction (C-Br-π) was detected between bromine and Phe112 at a distance of .36 nm, and this group was found to be stabilized in the CCR5 pocket by edge-to-face aromatic interaction. Furthermore, ethoxime reoriented after MD simulation and was located on new set of residues, such as, Phe109, Ser180, His181, Ser160, Phe166, Trp190, and Gln194.
Our MD simulations' identified results are in line with most of the previous mutational reports (Billick et al., 2004;Castonguay et al., 2003;Dragic et al., 2000;Kondru et al., 2008;Nishikawa et al., 2005;Paterlini, 2002;Seibert et al., 2006;Tsamis et al., 2003), but in contrast with the Shahlaei et al. (2011) report. In their report, it was noticed that the CCR5 drifted more compared to original inserted CCR5 structure in POPE bilayers, which is evident through the backbone RMSD plot, which shows CCR5 had $4.5 nm deviations in first phase of simulation. However, in all our investigated systems, backbone RMSDs (Figure 6(A)) was less than .42 nm. Moreover, inhibitors in the current MD study maintained essential salt bridge interactions with Glu283 (TM7), which are supported by the previous mutagenesis reports (Castonguay et al., 2003;Dragic et al., 2000;Nishikawa et al., 2005;Seibert et al., 2006;Tsamis et al., 2003). But, Shahleai et al. (2011) simulated model failed to show essential interaction with the most important Glu283 in CCR5 binding pocket. In addition, current simulated ligands (14, 25, and 37) (Figures 12  and 13) identified essential interaction in binding pocket with the important residues such as Trp86 (TM2), Tyr89 (TM2), Tyr108 (TM3), Phe112 (TM3), Ile198 (TM5), Tyr251 (TM6), Leu255 (TM6), and Gln280 (TM7), which is in line with previous reports (Castonguay et al., 2003;Dragic et al., 2000;Nishikawa et al., 2005;Seibert et al., 2006;Tsamis et al., 2003). The locations of the simulated inhibitors are on the extracellular surface of CCR5 and relevant to understand the mechanism of action. However, Shahleai et al. 2011 report proposed different binding pocket which could be deeper into the TMs rather than at ECL regions of CCR5, and that is why their simulations could not identify any of the previously reported residues.
We have been working on chemokine receptors using in silico modeling methods (Gadhe et al., 2010;Kothandan, Gadhe, Madhavan, & Cho, 2011). In particular, we have developed QSAR models for CCR5 inhibitors using ligand-based methods, such as CoMFA, CoMSIA, and HQSAR. More recently, we compared the binding sites of CCR2 and CCR5 to facilitate the development of dual antagonists targeting both CCR2 and CCR5 (Kothandan et al., 2012). This SAR study of OPPA derivatives concurred well with the MD simulated conformation of ligands (14, 25, and 37), which underscores the importance of the hydrophobicity, appropriate orientations of the 2,4-dimethyl-3-pyridine ring and bromophenyl group, and of the Z-isomer of oxime. Newly introduced residues in vicinity of the ethoxime moiety, such as Ser180, His181, Phe166, and Trp190 might be crucial in the CCR5 binding pocket and additional mutational studies on these residues are necessary to determine their contributions to CCR5 inhibition.
Here, we propose a general binding mode for the OPPA derivatives based on the MD simulation results of higher (14), medium (37), and lower (25) potent inhibitors. Interestingly, we find some trend for these inhibitors such as salt bridge interaction between basic nitrogen of ligand and acidic Glu283, which seems to be necessary for inhibitory activity. Also two aromatic pockets (pocket I, TM1-3 and pocket II, TM3-6) were connected by the central polar region in TM7, and all the simulated inhibitors show important interactions with the Trp86, Tyr89, Tyr108, Phe112, Ile198, Tyr251, Leu255, and Gln280 and Glu283. The results obtained by MD simulations of 14, 25, and 37 in CCR5 can be utilized to design potent inhibitors.

Conclusion
In this paper, we describe a homology model of CCR5 based on a recently reported CXCR4 template structure. Automated molecular docking was performed for the potent (14), medium (37), and lower (25) active CCR5 antagonists in a CCR5 model. The complexes (14-CCR5, 25-CCR5, and 37-CCR5) were simulated in a lipid environment to determine the effect of lipid on these complexes and to explore the structural features and binding mechanism of antagonists as a guide to inhibitor design.
Automated molecular docking simulation allowed us to propose a general binding mode for such inhibitors in the active site of CCR5 and to recognize the resi-dues important for inhibitory activity. Because binding of an inhibitor to a receptor is not a static process in vivo, we performed a long duration MD simulation for these complexes (14-CCR5, 25-CCR5, and 37-CCR5) immersed in lipid bilayers to sample its dynamic behavior. It was detected that ligand position and orientation changed after MD simulation, which indicates that this type of simulation is useful for determining the global conformations of ligands. Furthermore, MD results indicated that energetically ligands (14, 25, and 37) and CCR5 were stable throughout the simulation. In addition, it was observed that residue locations were changed after MD simulation and some new residues were introduced into the vicinity of the ligands. These newly introduced residues were Ala29, Ala90, Ser160, Phe166, Ser180, His181, and Trp190, and they were found to be crucial for ligand binding. It is evident from MD simulation that acidic Glu283 is an important residue for salt bridge contact with ligands. Our results are in line with previous mutational reports. MD modeling provided a detailed interaction mechanism between ligands (14, 25, and 37) and CCR5 and the results of our SARs study on other OPPA derivatives agreed well with the MD simulated conformations of 14, 25, and 37.
Good correlation was observed between 14 and SARs data, indicating that to design a high affinity binder the pyridine ring should be substituted with a hydrophobic group rather than a polar group, because it binds in a hydrophobic cavity. In addition, the bromophenyl ring interacts sterically with the pocket wall and any further bulky substitution at this region would decrease binding affinity. Furthermore, the Z-isomer showed more activity than the E-isomer, because it extended in a pocket lined by the residues Phe109, Ser180, His181, Phe166, Trp190, and Gln194. So far, no site-directed mutagenesis study has been conducted on the Ser180, His181, Phe166, or Trp190 residues, and thus, mutational studies are required to determine their contributions to CCR5 antagonism. It was also detected that after MD simulations, 25 and 37 were devoid of any interactions with the residues from TM1 and 4. The application of MD simulation after docking could be used to identify stable ligand binding mode in protein, and to aid understanding of the process. MD simulations have become a standard tool for the exploration of ligand-protein systems. Simulations are now being performed on large systems containing tens of thousands of atoms using more realistic boundary conditions and better sampling due to the longer simulation times used. Furthermore, realistic simulations of systems as complex as TM proteins embedded in lipid bilayers have become more feasible, and assist our understanding of ligand-protein interactions and provide a dynamic dimension to structural data.