ABSTRACT

CC chemokine receptor 5 (CCR5), a member of G protein-coupled receptors (GPCRs), not only plays a significant role in inflammatory responses, but also correlates with HIV-1 infection and cancer progression. Recently, blocking of CCR5 has been considered as an effective strategy in HIV-1/cancers therapy. So far, only Maraviroc has been approved by FDA in 2007, while the other CCR5 inhibitors have failed in their clinical trials. In this study, a highly selective structure-based pharmacophore model was constructed, validated, and applied for virtual screening to retrieve novel CCR5 inhibitors from NCI database. Finally, one potential CCR5 inhibitor candidate, NSC13165, was identified after molecular docking, molecular dynamics (MD) simulations, binding free energy analyses and ADMET prediction. Docking and MD simulation results not only suggested that NSC13165 reserves the common binding mode of the most known CCR5 inhibitors, but also provided important insights toward the allosteric inhibition mechanism of CCR5. The results of binding free energy analyses indicated that the binding affinity of NSC13165 is much better than that of Maraviroc and that van der Waals interaction is the key driving force during the binding process. ADMET prediction suggested that NSC13165 exhibits very low risk of causing lethal side effects. Altogether, our results strongly suggest that NSC13165 has great potential to serve as a novel CCR5 inhibitor, which may be further tested in vitro/in vivo as a drug target for HIV-1/cancers therapy or be used as a lead compound for improving its efficacy through chemical modifications.

Keywords:

CCR5 Inhibitor, HIV-1, Maraviroc, Virtual Screening, Molecular Dynamics Simulation

1. Introduction

CC chemokine receptor 5 (CCR5), a member of the seven-transmembrane domain G protein-coupled receptors (GPCRs) family, is mainly expressed on the surface of white blood cells and plays an important role in the human inflammatory responses to infection. Besides, CCR5 also involves in many autoimmune diseases, such as organ transplant rejection [1], diabetes [2], and rheumatoid arthritis [3]. Recent studies also indicated that CCR5 is crucial for both HIV-1 infection [4 , 5] and different kinds of cancer progressions [6]. During HIV-1 infection, two immune cell surface proteins, cluster of differentiation 4 (CD4) and CCR5, are most commonly adopted for entering host cells [7]. Previous studies have indicated that individuals carrying homozygous CCR5Δ32 deletion are protected from HIV-1 infection [8]. Therefore, blocking the function of CCR5 by CCR5 inhibitors has been considered as an effective and relatively harmless HIV-1 therapeutic strategy [9]. On the other hand, CCR5 and its natural binding substrate CCL5 (known as CCR5-CCL5 axis) have been recently proven to exhibit strong correlation with the progression of several cancers [6], such as breast cancer [10], prostate cancer [11], colon cancer [12], gastric cancer [13] and ovarian cancer [14]. In both breast and prostate cancers, CCR5 are overexpressed [11 , 15]. Previous experimental results indicate that CCR5 inhibitors can suppress the migration, invasive and metastasis of these cancer cells [16 , 17]. Thus, CCR5 inhibitors have become as an effective target for cancer therapy.

The discovery and implication of CCR5 in HIV-1 infection, cancer progression and several inflammatory diseases have triggered the development of various CCR5 inhibitors, such as Maraviroc [18], Aplaviroc [19] and Vicriviroc [20]. So far, Maraviroc, a competitive CCR5 inhibitor and negative allosteric modulator of CCR5, is the only CCR5 inhibitor received full FDA approval for HIV-1 treatment in 2007. Both Aplaviroc and Vicriviroc fail in their clinical trials due to hepatotoxicity [21] or lack of efficacy [22]. After the success of Maraviroc, many pharmaceutical companies and academic researchers have been dedicating to elucidate the inhibition mechanism of CCR5 [23] and to develop novel CCR5 inhibitors [24 , 25]. However, most of these CCR5 inhibitors fail due to poor pharmacokinetic profiles in animal studies, inhibition of the metabolic enzyme CYP2D6, or unsafe level of human Etheràgo-go-Related Gene (hERG) inhibition [26]. Cenicriviroc, a dual CCR2 and CCR5 inhibitors, is the only CCR5 inhibitor currently in clinical trials [25]. Thus, it is highly desired to develop novel CCR5 inhibitors to be used as Maraviroc alternatives in the pharmaceutical market for HIV-1/cancers therapy.

In the past few years, due to the lack of CCR5 crystal structure, the strategies for developing novel CCR5 inhibitors were limited to ligand-based level, such as quantitative-structure activity relationship (QSAR) [27], ligand-based pharmacophore modeling [28], shape-based virtual screening [29], structure-activity relationships (SAR) analysis [30] and lead modification [18]. In 2013, the first and only crystal structure of CCR5-Maraviroc complex was reported [23], which makes the investigation of protein-ligand binding interaction on the molecular level and accurate structure-based drug design become possible. In this study, a highly selective structure-based pharmacophore model was constructed, validated and used to screen potential and novel CCR5 inhibitors with excellent efficacy, good drug-likeness and less toxicity from the National Cancer Institute (NCI) database. Molecular docking was performed to provide detail binding interaction and geometry prediction of the screened compounds. Molecular dynamics (MD) simulations and binding free energy analyses were conducted to gain insights into the binding mode for the selected compounds with excellent binding affinity. Previous reports have evidenced that unexpected toxicity, poor pharmacokinetic properties and metabolism problems are the major problems in the clinical trials [20 , 21 , 26]. To date, computational ADMET (absorption, distribution, metabolism, excretion and toxicity) prediction and drug-likeness filters have become a widely used tool in virtual screening. Therefore, Lipinski rule of five, Veber rule and comprehensive ADMET prediction [31] were also executed in this study to reduce the risk of failure in clinical trials. The selected CCR5 inhibitors may be further tested in vitro/in vivo as a drug target for HIV-1/cancers therapy and provide important information for designing novel CCR5 inhibitors for clinical use.

2. Materials and Methods

2.1. Receptor, Ligands, Database Preparation

The crystal structure of the CCR5-Maraviroc binding complex (PDB ID: 4MBS) [23], at a resolution of 2.71 Å in the inactive state, was downloaded from the RCSB Protein Data Bank. In 4MBS, rubredoxin (residues 1001 - 1054) was inserted into the bottom of CCR5 between residues 223 and 227 (Figure 1(A)), roughly 45 Å away from the binding pocket and had a negligible effect on the protein-ligand binding interaction and binding free energy calculation. Thus, rubredoxin was remained in our simulation processes but excluded from our analysis. Other point mutations of the crystal structure and incomplete residues were fixed and refined using the Prepare Protein protocol embedded in BIOVIA Discovery Studio 2016 (DS 2016).

A group of 41 known CCR5 inhibitors, characterized by their excellent experimental performances (IC₅₀ ≤ 10 nM), including Maraviroc, Ancriviroc, Aplaviroc, INCB-9471 and SCH-350634 (Figure 2), were collected from previous literatures [32 - 38]. A decoy set based on the chemical structure of Maraviroc was generated from the online server Database of Useful Decoys: Enhanced (DUD.E) [39], which provides a set of 393 decoys that are available from the online ZINC database. The publicly available NCI database (>265,000 structures) was downloaded from the NCI website for virtual screening. All of the known inhibitors, decoys and NCI database were prepared under the ionization state at pH 7.0 using the Prepare Ligand protocol of DS 2016.

2.2. Molecular Docking

In this study, both LibDock and CDOCKER embedded in DS 2016 and Schrödinger Maestro 11.3 were employed for molecular docking experiments and data analyses. LibDock was used for rapid rigid docking of the screened compounds while CDOCKER was used for more accurate and precise flexible docking. Before our docking experiments, Maraviroc was removed from the refined 4MBS structure and the binding site was defined as a sphere radius of 10 Å from the position of Maraviroc. Most of the docking parameters were set to default for both LibDock and CDOCKER, except for the conformation parameter, which was set to Fast for LibDock.

2.3. Structure-Based Pharmacophore Model Generation and Validation

To obtain comprehensive information of the ligand-receptor binding complexes, each of the 40 known inhibitors (except Maraviroc) was docked into the CCR5 binding site and a maximum of 10 binding poses was generated using CDOCKER. Among these docking poses, only 64 out of 400 poses that contained similar CCR5 binding interactions [23 , 40 - 45] along with the refined 4 MBS were selected for the generation of structure-based pharmacophore models using the Common Feature Pharmacophore Generation of DS 2016. For each binding complex, a maximum of 10 pharmacophore models with 4 - 6 pharmacophore features were automatically generated using default parameters. Finally, a total of 38 structure-based pharmacophore models were constructed.

To select the best pharmacophore model among these 38 structure-based pharmacophore models, a test 3D database comprising both 41 known CCR5 inhibitors and 393 decoys was constructed using Build 3D Database (maximum 255 conformations for each compound) of DS 2016. Then, these models were used for screening the test 3D database using Screen Library with default parameters except that energy threshold was set to 5 in DS 2016. The screening results of each model were evaluated using the Güner-Henry scoring method [46]. Finally, the best pharmacophore model (named as Model_1), which was constructed from the CCR5-inhibitor A complex and yielded the highest GH score of 0.866, was selected for virtual screening.

2.4. Virtual Screening

To enhance the drug-likeness of our screening result, the prepared NCI database was filtered by the Lipinski rule of five and Veber rule using DS 2016. Then the filtered NCI database was built into 3D database (maximum 255 conformations for each compound) for virtual screening by Model_1 using Screen Library protocol of DS 2016. Compounds that successfully mapped all the pharmacophore features of Model_1 were considered as potential hits, which were further submitted to molecular docking studies with LibDock and CDOCKER. Finally, only compounds that map well with Model_1, have higher docking scores comparing to inhibitor A, and contain similar binding mode and critical interactions with the key residues were considered as potential CCR5 inhibitor candidates.

2.5. Maintaining the Integrity of the Specifications

After virtual screening, the CDOCKER binding complex of each potential CCR5 inhibitor candidate was used as the initial conformation for MD simulations to ensure the optimal molecular interactions with CCR5. All MD simulations were performed using GROMACS 2016.3 [47] and combined with NVIDIA CUDA 8.0, Antechamber in AMBER tool 2017 [48], ACPYPE [49]. Each complex was solvated and neutralized in a 10.0 Å dodecahedron periodic box with TIP3P water molecules. Then additional sodium and chloride ions were added to reach 0.15 M salt concentration as physiology state at saline. For ligand and protein in each complex, the general AMBER force field (GAFF) [50] and AMBER ff99SB-ILDN force field [51] were applied. Temperature was set to 310 K at 1.0 atm pressure. For each solvated system, minimization was performed before two steps of restrained equilibrations. Restrained NVT equilibration was carried out for 200 ps and then restrained NPT equilibration was applied for another 200 ps. Finally, the non-restrained MD production run was performed for 50 ns in MD I and extended to 200 ns in MD II.

2.6. Binding Free Energy Calculations

The molecular mechanics energies combined with the Poisson-Boltzmann and surface area continuum solvation (MM/PBSA) method was applied to calculate the binding free energies through g_mmpbsa for GROMACS-5.1.X [52]. Although entropy is omitted in MM/PBSA calculation, it is still a reliable tool to compare the binding free energy of different ligands in the same protein binding pocket. For each complex, one trajectory was extracted from MD production for g_mmpbsa calculation. The solvent accessible volume only (SAV-only) nonpolar model was applied for nonpolar solvation energy calculation. Other parameters were set to default except the temperature was set to 310 K. To perform a rapid binding energy comparison in MD I, a total of 100 snapshots were collected from the last 10 ns MD simulation trajectory. For a more precise and reliable binding energy analysis in MD II, a total of 2000 snapshots were collected from the last 20 ns MD simulation trajectory.

2.7. ADMET Prediction

Many drug candidates have failed in the early stage of clinical trials due to poor ADME or toxicity profile [26]. Therefore, the online server admetSAR [31], a reliable and widely used free tool for evaluating chemical ADMET properties, was adopted to evaluate the ADMET properties of the screened CCR5 inhibitor candidates in this study. In admetSAR, ADMET properties are predicted by the simplified molecular-input line-entry systems (SMILES) for each selected compound. Compounds with better predicted ADMET properties comparing to Maraviroc were finally selected as potential CCR5 inhibitor candidates.

3. Results and Discussion

During the early stage of CCR5 inhibitor development, lethal side effects, bad ADMET properties and lack of efficacy are the major problems in clinical trials [26]. Previously, the strategies for developing novel CCR5 inhibitors are limited to ligand-based level due to the lack of CCR5 crystal structure. Since the first and only X-ray crystal structure of CCR5-Maraviroc complex, 4MBS, was reported in 2013 [23], precise molecular level studies through molecular docking [42 , 53] and MD simulation [40 , 41 , 43 , 44] have become possible. However, there is still a lack of structure-based virtual screening study for the identification of novel CCR5 inhibitors. To the best of our knowledge, this study is the first attempt to discover potential CCR5 inhibitors with excellent efficacy, good drug-likeness and less toxicity by a combination of structure-based pharmacophore modeling, virtual screening, molecular docking, MD simulations, binding free energy analyses and ADMET prediction.

3.1. Structure-Based Pharmacophore Model Generation and Validation

Structure-based pharmacophore models are derived from the structure of the target protein-ligand complex by investigating important interaction and excluded volumes in the binding site and translated into pharmacophore features. Comparing to 3D QSAR and ligand-based approaches, structure-based approach has advantages of more precise prediction, good filter to remove decoys and finding structurally novel ligands [54]. In this study, the best structural-based pharmacophore model, Model_1, consisting of 4 pharmacophore features, including 2 hydrophobic features and 2 hydrogen bond acceptors (Figure 1(B)), was constructed from the CCR5-inhibitor A (IC₅₀ = 0.15 nM) binding complex (Figure 1(A) and Figure 1(C)) [36]. As shown in Table 1, Model_1 can successfully recognize 40 out of 41 known CCR5 inhibitors (sensitivity = 0.976) and correctly exclude 386 out of 393 decoys (specificity = 0.982). To comprehensively evaluate the statistical quality and overall fitness of Model_1, rigorous validation was carried out using Güner-Henry (GH) scoring method. GH score higher than 0.7 suggests a very good and reliable model [46]. As can be seen in Table 1, Model_1 possess an excellent GH score of 0.866. The above validation results confirmed that Model_1 exhibits great quality to be used as a 3D query for the further virtual screening process.

* $GH=\left(\frac{Ha\left(3A+Ht\right)}{4HtA}\right)\left(1-\frac{Ht-Ha}{D-A}\right)$ . GH score higher than 0.7 suggests a very good pharmacophore model.

3.2. Virtual Screening

Virtual screening is a useful technique to discover potential candidates that are able to block or trigger the activity of the selected drug target. Comparing to docking-based virtual screening, pharmacophore-based virtual screening considers more about protein flexibility and reduces the impact of insufficiently designed or optimized docking scoring functions [55]. Another advantage is that it can identify compounds with different scaffolds but sharing similar biological activities, also known as “scaffold hopping” [56], which can bring the advantages of improving physicochemical and ADMET properties, enhancing binding affinity and selectivity, and/or discovering patentable analogues. Here, the well validated Model_1 was adopted to screen potential CCR5 inhibitors from the NCI 3D database. Finally, a total of 293 compounds that fitted well with all the features of Model_1 were selected and considered as hits. Then, these compounds were submitted to the further molecular docking studies.

3.3. Molecular Docking

Molecular docking was performed to investigate the binding conformation of each selected compound. To obtain the hits with better binding affinity, the LibDock score (118.396) and CDOCKER score (CDOCKER interaction energy, −50.6996 kcal・mol⁻¹) of inhibitor A were set as threshold values (Figure 3). LibDock was firstly applied to reduce the number of hit compounds, resulting in a total of 70 compounds out of the 293 selected compounds from the previous virtual screening process. Comparing to LibDock, CDOCKER, a MD simulated-annealing-based docking algorithm, provides a more accuracy docking result. To enhance the precision of binding structures and reduce false positives, these 70 compounds were docked into the binding pocket using CDOCKER. Finally, 18 compounds with higher CDOCKER scores comparing to that of inhibitor A were selected.

Recent studies have indicated that CCR5 inhibitors are allosteric inhibitors, sharing a similar binding site on CCR5 [45 , 57]. The known CCR5 inhibitors all establish an extensive binding interaction network with several CCR5 key residues in the binding pocket, such as Leu33, Tyr37, Trp86, Tyr108, Phe109, Phe112, Gln194, Thr195, Ile198, Trp248, Tyr251, Glu283, Thr284 and Met287 [23 , 40 , 44 , 45 , 57 - 60]. Through the site-specific mutagenesis studies, Trp86, Tyr108, Ile198, Tyr251 and Glu283 have been identified as critical residues for CCR5 inhibition [45 , 57]. Among these critical residues, Glu283 is the most important residue for CCR5 inhibition. Besides, the π-π stacking interactions between Trp86, Tyr108 and Phe109 residues of CCR5 and ligands also help to stabilize the binding process [43]. Based on the above information, only 7 compounds (Figure 3), which all exhibit similar binding interactions with these CCR5 key residues, were selected as potential CCR5 inhibitors for further analysis. These 7 compounds were named as compounds 1-7 according to their CDOCKER scores (from the highest to the lowest).

3.4. The Analyses of Compounds 1-7

3.4.1. MD Simulation I (MD I)

Molecular dynamics (MD) simulations are able to predict the atomic motion of protein-ligand complex, such as identification of the binding site, validation of the docking results and calculation of the binding free energy, and thus play an important role in drug discovery. To determine and compare the binding stabilities of compounds 1-7, inhibitor A and Maraviroc, their CDOCKER docking complexes were subjected to 50 ns MD simulations (MD I) in an explicit hydration environment.

The dynamic stabilities and MD simulation trajectories of these complexes were explored by their backbone root-mean-square deviations (RMSD). Backbone RMSD with small fluctuations and constant value suggests a stable binding system. As shown in Figure 4(A), compounds 1, 2, 7, inhibitor A and Maraviroc all exhibited more stable MD trajectories in the last 10 ns, with average RMSD values of 2.737 Å, 2.726 Å, 2.880 Å, 2.759 Å and 2.786 Å, respectively. In contrast, compounds 3, 4, 5 and 6, exhibited more unstable MD trajectories with higher average RMSD values of 3.320 Å, 3.547 Å, 3.147 Å and 3.033 Å, respectively, in the last 10 ns of MD. These MD results suggest that the stabilities of the dynamic equilibriums of compounds 1, 2 and 7 were better than those of compounds 3, 4, 5 and 6.

To explore the protein residues flexibility throughout the MD simulation, the root-mean-squared fluctuations (RMSF) of individual CCR5 residues were calculated. As shown in Figure 4(B), compounds 1-7 exhibited similar RMSF distributions comparing to inhibitor A and Maraviroc. In addition, the conformational changes of the above mentioned CCR5 key residues (Leu33, Tyr37, Trp86, Tyr108, Phe109, Phe112, Gln194, Thr195, Ile198, Trp248, Tyr251, Glu283, Thr284 and Met287) were all very small (Figure 4(C)). Among these 7 selected compounds, compound 1 showed the lowest RMSF distribution, indicating that it can sufficiently fix most of the CCR5 key residues comparing to the allosteric inhibitor Maraviroc (Figure 4(C)) [23 , 45 , 57]. The results of RMSF analyses suggest that compound 1 exhibits better binding stability with CCR5 key residues.

3.4.2. Binding Free Energy Calculation

Binding affinity describes the interactions of the ligands with the binding pocket, which can be evaluated through binding free energy calculation. In general, compounds with high binding affinities have a higher degree of ligand occupancy in the binding site than compounds with low binding affinities. Thus, binding free energies (ΔG_bind) of compounds 1-7, inhibitor A and Maraviroc were further calculated to determine their binding affinities towards CCR5 using g_mmpbsa in this study and the results are shown in Figure 3. Comparing to Maraviroc (−72.21 kcal・mol⁻¹) and inhibitor A (−69.41 kcal・mol⁻¹), compounds 1 (−83.43 kcal・mol⁻¹), 3 (−93.77 kcal・mol⁻¹), 4 (−93.56 kcal・mol⁻¹) and 6 (−102.11 kcal mol⁻¹) all showed significant lower binding free energies, suggesting that these compounds all exhibit better binding affinities than those of Maraviroc and inhibitor A.

3.4.3. ADMET Prediction

Prediction of the compound ADMET properties describes the disposition of a pharmaceutical compound within the human body. Considering that most of previously developed CCR5 inhibitors failed in their clinical trials due to bad ADMET properties [26], filtering and optimization of ADMET properties in the early stage of the drug discovery are helpful to avoid failure. In this study, the ADMET properties of compounds 1-7 and Maraviroc were predicted and the results are presented in Table 2. For absorption, compounds 1, 3, 4, 5, 6, 7 and Maraviroc are able to possess good human intestinal absorption (HIA), suggesting oral delivery possibility. For distribution, compounds 1, 3, 4, 5 and Maraviroc can penetrate the blood-brain barrier (BBB), which is necessary for treating HIV-1 infection in the central nervous system (CNS) [61]. For toxicity, compounds 1, 2, 5, 6 and 7 are both non-mutagenic (AMES toxicity test) and non-carcinogenic (carcinogenic profile). For solubility, compounds 1, 3, 4, 6 and 7 had better aqueous solubility than that of Maraviroc.

Since the US Food and Drug Administration (FDA) published the first in vitro/in vivo drug interaction guidance documents in 1997 and 1999, respectively, off-target transporter interactions and drug-drug interactions (DDIs) have become the major challenges for drug development [62]. In human heart, hERG potassium channels are essential for regulating electrical activity [63]. The blockage of hERG channel by drug causes lethal long QT syndrome (LQTS), which leads to withdraw from the clinical trials or markets [64]. The results of hERG inhibitor prediction I/II showed that both compounds 1 and 2 are hERG non-inhibitors, involving lower LQTS risk than the others. CYP450 enzymes-based unanticipated DDIs and drug metabolism problems are also a common cause of adverse drug events (ADE) [62]. Interacting with CYP450 enzymes may lead to DDIs, hepatotoxicity, reactive metabolite formation, idiosyncratic adverse drug interactions, and/or even loss the efficacy of drugs [65]. Only compounds 1 and 2 are non-inhibitors to all of the CYP450 superfamily of enzymes and have low CYP450 inhibitory promiscuity. All of the screened compounds are non-inhibitors of renal organic cation transporter (OCT), decreasing the possibility of transporter-mediated renal DDIs [66].

Based on the ADMET prediction, compound 1 not only exhibits better ADMET properties than that of Maraviroc, but also shows great improvement in lowering the risk of causing both lethal CYP450 superfamily inhibition and hERG inhibition. Combining with the results of MD I, binding free energy calculation and ADMET prediction, compound 1 is considered as the most potential CCR5 inhibitor candidate.

^a+ for positive and − for negative. ^b+ for penetrable and − for non-penetrable. ^c+ for inhibitor and − for non-inhibitor. ^d+ for toxic and − for non-toxic. ^e+ for carcinogens and − for non-carcinogens.

3.5. The Comparison of Compound 1 and Maraviroc

3.5.1. The Analysis of Binding Interactions

To gain insights into the binding interaction between CCR5 and its inhibitors, the initial binding complexes of compound 1 and Maraviroc were investigated through their 3D structures and 2D-interaction maps (Figure 5). According to Figure 4(A) and Figure 4(C), both compound 1 and Maraviroc occupied the similar location in the binding pocket. The aromatic rings of compound 1 and Maraviroc reached deep into the bottom of the binding pocket and formed hydrophobic interactions with several key CCR5 residues, which is consistent with previous studies [23].

The 2D-interaction maps showed that Maraviroc formed one hydrogen bond with Glu283 (Figure 5(B)); while compound 1 formed three hydrogen bonds with Tyr251, Glu283 and Thr284 (Figure 5(D)). For Maraviroc, π-π stacking interactions were occurred with Trp86, Tyr108, Phe109, Phe112, Tyr251 (Figure 5(B)); while for compound 1, π-π stacking interactions were formed with Tyr108 and Phe112 (Figure 5(D)). In addition, both compound 1 and Maraviroc formed hydrophobic interaction with Trp248, which may prevent CCR5 from structural changing into active conformation [23].

The results of binding interaction analyses indicate that compound 1 and Maraviroc share similar binding interactions with CCR5 although compound 1 exhibits a totally different chemical structure than Maraviroc (see Table 2), which is known as “scaffold hopping” [56]. With such a totally different chemical structure, compound 1 may bring the advantage of having better ADMET properties than Maraviroc (see Table 2). In order to gain insights into the effects of scaffold hopping on the binding stability towards CCR5, compound 1 and Maraviroc were subjected to further MD II and binding free energy calculation and decomposition analyses.

3.5.2. MD Simulation II (MD II)

To ensure the reliability and precision of MD simulations and binding free energy calculations, the MD production of both compound 1 and Maraviroc were extended from 50 ns in MD I to 200 ns in MD II. In order to explore the dynamic stability of CCR5-inhibitor complexes, the RMSD values of the entire systems and RMSF values of each CCR5 residue were investigated. The total binding free energy were also decomposed into four individual free energy components and break down into the contribution from each residue to elucidate the binding mechanism in detail. Comparing to molecular docking, these MD-based analyses allow to provide the important insights of the allosteric mechanism of CCR5 inhibition and validate the binding affinity of the selected CCR5 inhibitor.

1) RMSD Analysis

To monitor the overall stability of the binding complexes, the RMSD values of the entire protein-ligand backbone atoms were calculated. As shown in Figure 6(A), both compound 1 and Maraviroc retained excellent equilibrium in the extended MD II production runs. The average RMSD value of compound 1 and Maraviroc in the last 10 ns were 2.980 Å and 2.942 Å, respectively. Both compound 1 and Maraviroc had reached an equilibrium phase with very small fluctuations over the entire 200 ns MD simulations, suggesting excellent binding stability with CCR5.

2) RMSF Analysis

The RMSF provides a measure of CCR5 local flexibility after binding with allosteric CCR5 inhibitors. The results of RMSF analyses for compound 1 and Maraviroc are given in Figure 6(B). It shows that both compound 1 and Maraviroc had similar RMSF distribution for most of the residues, suggesting similar conformation changes toward CCR5. The average RMSF values of compound 1 and Maraviroc were 1.95 Å and 2.03 Å, respectively. To investigate the flexibility of each CCR5 residue in the allosteric binding mechanism, the CCR5 regions with RMSF values less than 3 Å were investigated and defined as “stable regions” (Figure 6(B)) [67]. Interestingly, the residues located in the stable regions do not necessarily correspond to the residues near CCR5 binding site, implying that the binding of CCR5 inhibitors not only stabilize CCR5 binding site, but also propagate to the other regions of CCR5, including 7 transmembrane domains (TM 1 - 7) and 2 extracellular loops (ECL1 and ECL2). Most of the residues located out of the stable regions belong to the 2 terminal helices and extracellular loop 3 (ECL3).

The RMSF values of CCR5 key residues (Leu33, Tyr37, Trp86, Tyr108, Phe109, Phe112, Gln194, Thr195, Ile198, Trp248, Tyr251, Glu283, Thr284 and Met287) were also investigated. As shown in Figure 6(C), the key residue RMSF values of compound 1 were all smaller than those of Maraviroc, including the critical residues Trp86, Tyr108, Ile198, Tyr251 and Glu283 for CCR5 inhibition [45 , 57]. For the most important residue Glu283, the RMSF value of compound 1 (0.725 Å) was significantly smaller than that of Maraviroc (1.376 Å). Thus, the results of RMSF analyses suggested that compound 1 exhibits similar conformation changes toward CCR5 binding in a more stable binding effect in the binding site comparing to Maraviroc.

*all energies were calculated by g_mmpbsa and in kcal・mol⁻¹.

3) Binding Free Energy Calculation

Binding free energy was calculated by the MM/PBSA method, which provides more detailed information on the interactions of compound 1 and Maraviroc towards CCR5. Negative values for binding free energies indicate high thermodynamic stability and good binding affinity. Higher binding affinity of a drug target indicates lower concentration required during treatment, which can reduce the risk of side effects caused by binding to off-target proteins. The results of the predicted binding free energies and energy components of compound 1 and Maraviroc from MD II are shown in Table 3. It shows that compound 1 had lower binding free energy (−88.67 kcal・mol⁻¹) comparing to that of Maraviroc (−75.14 kcal・mol⁻¹), indicating that compound 1 exhibits higher binding affinity than Maraviroc.

4) Binding Free Energy Decomposition Analysis

Binding free energy can be decomposed into four individual components (ΔG_vdw, ΔG_elec, ΔG_polar and ΔG_nonpolar) to investigate the effect of these energy terms toward binding. Table 3 showed that van der Waals interactions (ΔG_vdw), electrostatic interaction (ΔG_elec) and nonpolar solvation energy terms (ΔG_nonpolar) provided favorable contributions to the binding, whereas polar solvation energy terms (ΔG_polar) opposed the binding. Table 3 also showed that van der Waals interactions are the major driving force for compound 1 and Maraviroc towards the binding processes. The formation of the strong favorable van der Waals interactions with CCR5 inhibitors towards the binding process is well correlated with the fact that most of CCR5 key residues are hydrophobic [40]. The reason for the formation of unfavorable polar solvation energy terms is that blocking the interaction of negatively charged key residue Glu283 from solvent would enhance the destabilization polar solvation term [43]. Thus, increasing the favorable van der Waals interactions and reducing the unfavorable polar solvation energy terms of CCR5 inhibitors through chemical modifications would be a promising strategy for designing novel CCR5 inhibitors with better binding affinity.

4. Conclusion

Although CCR5 inhibitors have been proven to exhibit great potential in the treatment of HIV-1 infection and cancer progressions, the only FDA-approved Maraviroc still has a risk of causing life-threatening side effects, including heart attack, skin reaction, liver damage, and allergic reaction [68]. Thus, there is a highly desire for discovering novel CCR5 inhibitors with new chemical skeleton, better binding affinity and relatively harmless side effects. In this study, the constructed structure-based pharmacophore model Model_1 was shown to be excellent with its high sensitivity and specificity to identify potential CCR5 inhibitors from the NCI database. NSC13165 (compound 1), which exhibits a different chemical scaffold than Maraviroc, was identified to possess a higher binding affinity toward CCR5 comparing to Maraviroc through a series of computational approaches. Besides, NSC13165 also shows very low risk in several common lethal side effects through the ADMET prediction, such as off-target transporter interactions, drug-drug interactions and cell toxicities, thus it can be further tested through in vitro/in vivo studies or used as a lead compound for improving its efficacy, selectivity, or other pharmacokinetic parameters. Our binding free energy calculation results also suggest that enhancing van der Waals interaction, electrostatic interaction, nonpolar solvation energy and/or decreasing polar solvation energy term can significantly improve the binding affinity of CCR5 inhibitors, leading to a prospective strategy for further chemical modification, SAR, QSAR studies. Finally, the highly selective pharmacophore model can be further adopted to screen more potential CCR5 inhibitors from the other molecular databases, such as ZINC, MCULE and ChEMBL.

Acknowledgements

The authors thank the Ministry of Science and Technology (MOST 104-2221-E-027 -073 -MY3) and National Taipei University of Technology and Taipei Medical University (NTUT-TMU-101-10 and NTUT-TMU-102-10) for their financial supports.

CONFLICTS OF INTEREST

The authors declare no conflicts of interest regarding the publication of this paper.

References

NOTES

*Equal contribution.