Original Research ARTICLE
Discovery and Identification of Pyrazolopyramidine Analogs as Novel Potent Androgen Receptor Antagonists
- School of Pharmacy, Lanzhou University, Lanzhou, China
Androgen receptor (AR), an important target in the current androgen derivation therapy, plays a critical role in the development and progress of prostate cancer (PCa). Nonsteroidal antiandrogens, such as enzalutamide and bicalutamide, are commonly used in clinic to treat PCa. Though they are very effective at the beginning, drug resistance problem appears after about 18 months. One of the reasons is that these antiandrogens share similar structure skeleton. Therefore, it is urgent to discover novel antiandrogens with different skeletons for resistance problem. Herein, we combined structure- and ligand-based methodologies for virtual screening chemical databases to identify potent AR antagonists. Then the cytotoxic activities of the screened hit samples were evaluated by using LNCaP prostate cancer cells. Virtual screening and biological evaluation assay results suggest that several chemicals with novel pyrazolopyrimidine skeleton can inhibit the proliferation of prostate cancer cells with similar, or even higher, bioactivities to bicalutamide. AR reporter gene assay experiments proved that Compound III showed potential antagonistic effects. In addition, molecular dynamics simulations results proved that Compound III can properly bind to AR and prevent helix 12 (H12) from closing to distort the formation of activation function 2 (AF2) site, resulting in the invalid transcription. Hence, pyrazolopyrimidine was discovered as a novel, potent and promising antiandrogen skeleton deserved to be further studied.
According to the latest cancer statistics 2016 (Siegel et al., 2016), prostate cancer (PCa) is the second most common cancer among males around the world. The estimated death rate of PCa is 8%, which is just lower than the leading lung cancer, but the estimated new cases of prostate cancer become highest among the diagnosed cancers. Androgen receptor (AR) plays a critical role not only in the development of prostate cancer but also in the progress of the advanced castration states (Tilley et al., 1994; Taplin et al., 2003; Jernberg et al., 2017). Therefore, effective suppression of AR activity remains mainstream therapeutic schemes to the treatment of advanced, recurrent and metastatic prostate cancer.
AR, a class of ligand-activated transcription factor, is a member of the steroid and nuclear receptor (NR) superfamily (Evans, 1988). The AR structure consists of four basic elements: N-terminal domain (NTD), DNA-binding domain (DBD), hinge region and ligand-binding domain (LBD), that is a highly structurally conserved throughout the NR superfamily (Gao et al., 2005). Natural hormone testosterone (T) and dihydrotestosterone (DHT), binding to AR LBD, are the endogenous ligands of AR. The main mechanism of androgen action is to regulate the gene expression, change the level of specific proteins in cells and control cells behavior (Fix et al., 2004). Once bound to AR, androgens play pivotal roles in the sexual development, function, musculoskeletal growth of male and the progress of prostate cancer.
The standard treatment of prostate cancer involves androgen derivation therapy in conjunction with small molecule antiandrogens that block AR signaling (Hodgson et al., 2007). Antiandrogens compete with DHT for the binding to AR, inhibiting AR transactivation through a variety of mechanisms, including disruption of nuclear localization, interruption of DNA binding and interference with co-activator recruitment (Tran et al., 2009; Sadar, 2012). Unfortunately, most patients receiving antiandrogen therapy eventually develop drug resistance problem indicated by the rising level of serum prostate-specific antigen (PSA), leading to the lethal disease state termed castration-resistant prostate cancer (CRPC) (Denmeade et al., 2003). To deal with this situation, more efforts have been concentrated on design and discovery of new AR antagonists with novel skeleton to overcome drug resistance in CRPC.
In the current study, a combined structure- and ligand-based screening strategy was applied to virtually screen a big commercial chemical database to fish potential AR antagonists. Considering the published AR protein crystal structures are all in the agonistic manner (Bohl et al., 2005a,b), we firstly constructed its antagonistic conformation by using homology modeling according to the reported process (Liu et al., 2012). Structure-based docking method was used to find out chemicals that could bind to AR with high affinities. Then, a set of strictly validated QSAR models were used to predict the in silico antiandrogen abilities of the screened chemicals. Subsequently, the cell proliferation and AR reporter gene assay were used to evaluate the biological activities of the filtered molecules. Finally, molecular dynamics (MD) simulations and molecular mechanics Generalized Born (GB) surface area (MM-GBSA) calculations were applied to explore the interaction mechanics between pyrazolopyrimidine and ARs including wild type (WT), F876L and H874Y mutant types.
Results and Discussion
In order to effectively dig out chemicals with high antiandrogen potency, a combined structure- and ligand-based strategy was applied in this work. The flowchart of the combined virtual screening procedure was shown in Figure 1.
Figure 1. The flowchart of the combination of structure- and ligand- based virtual screening strategy.
Structure-Based Virtual Screening
The AR ligand binding domain represents a hydrophobic pocket, which can accommodate various chemicals of different sizes and chemotypes. Taking into account that the crystal structure of AR with antagonistic manner (Anti-AR) is unclear, homology modeling was applied to construct its Anti-AR conformation using the same method as previously reported (Liu et al., 2012). Then the AR crystal structure (PDB ID:1T65) and Anti-AR were both used to virtually screen the ChemDiv database by using Glide standard precision (SP) docking program. Glide SP, a semi flexible docking method, performs exhaustive sampling and is the recommended balance between speed and accuracy.
ChemDiv is now the recognized global leader in discovery chemistry with the industry's largest, most diverse, and most pharmacologically-relevant commercial collection of 1.6 million individually crafted, lead-like, drug-like small molecules etc. All chemicals that can be docked into the ligand binding domains of both 1T65 and Anti-AR were reserved. Subsequently, these compounds were further filtered following Lipinski's rules of five, including less than five hydrogen-bond donors, ten hydrogen-bond acceptors, 500 Da of molecular weight and AlogP of 5 in Discovery Studio 2.5, to keep compounds with potential bioavailability. Finally, 3,149 compounds were reserved in Database 2 for further screening, as shown in Figure 2.
Ligand-Based Virtual Screening
Quantitative structure-activity relationship (QSAR) is an analytical methodology that can be used to interpret the quantitative relationship between molecular structures and corresponding biological activities. Most importantly, after thorough validation, QSAR models can be used to quantitatively predict the activities of new chemicals. In this study, four QSAR models constructed by using multiple linear regression (MLR) method were employed to in silico evaluate the bioactivities of the filtered ChemDiv chemicals contained in database 2.
The Development and Validation of QSAR Models
The linear MLR model Y1 was previously reported from our group (Wang et al., 2015). Other three models Y2, Y3, and Y4 were newly built. Experimental and predicted antiandrogenic activities in these three models are listed in Supporting Information Tables S1–S3. Here the antagonistic activities (IC50) were experimentally tested from the AR reporter gene assays. The Y in these QSAR models stands for the negative logarithmic unit of the biological activities, pIC50. These four linear equations were list as follows. The corresponding parameters, demonstrating the model performances, were listed in Table 1. The descriptors used to build the four linear model and corresponding meanings were listed in Table 2.
Y1 = 1.188GATS7v − 1.108BEHp7 − 7.281E2u + 1.453HATS4u − 2.647H6m + 28.810R6u+ − 46.815R8u+ + 12.157
Y2 = −2.89IC5 + 1.01GATS5e − 3.17DISPp − 12.99HATS3u + 27.13
Y3 = 5.601IVDE − 1.70C-009 + 1.11BLTF96 + 3.02
Y4 = 3.18 IC1 − 0.46F05[N-F] + 1.47R4u + 0.409 Depressant-80 + 3.85HATS7m − 7.15
In Table 1, the parameter R2 describes the fitting ability and represents the reliability and stability of a QSAR model. Leave-one-out cross-validated correlation coefficient () is the one of internal validation methods which can avoid the risk of overfitting and the possibility of overestimating predictability. From the above parameters listed in Table 1, it can be seen that the fitting abilities of these MLR models are all high enough with R2 greater than 0.75. The predictive abilities were evaluated by means of four external validation parameters (Gramatica and Sangion, 2016), (www.oecd.org/dataoecd/33/37/37849783.pdf), (Schüürmann et al., 2008), (Consonni et al., 2009) and CCC (Lin, 1989; Chirico and Gramatica, 2011, 2012), and the formulas of parameters were shown in the Table S4. All these parameters are high enough to guarantee the predictive ability of these models. Especially the parameter CCC values are all beyond 0.85 as suggested in the literature (Chirico and Gramatica, 2012), which is a highly reliable parameter to guarantee the external predictability of a model. All the root mean square error (RMSE) values for the training set and prediction set are similarly low in each model. After thorough and rigorous validation, all these four QSAR models were proved to be reliable, stable and predictive and can be used to predict the bioactivities of new chemicals.
QSAR-Based Virtual Screening
After models construction and validation, these four models were used to calculate the in silico antiandrogenic bioactivities of the filtered compounds in database 2. Firstly, these 3,149 molecules were imported to DRAGON (DRAGON for Windows, version 5.5)1 program to calculate all molecular descriptors. Then the descriptors involved in these four QSAR models were extracted and submitted to the equations Y1-Y4 to calculate the biological activity of each molecule. Subsequently, all these chemicals were ranked by the calculated biological activities, considering simultaneously the results from four models. Only molecules with high potency, predicted bioactivities from four models greater than 7, can be retained for the succeeding process. Finally, 396 compounds were reserved in database 3 for further structural analysis and experiential analysis.
The 396 different compounds can be divided into 8 groups according to their skeleton structures. But some of them contain poisonous substructures from the point view of drug, some of them are easy to be hydrolyzed, or some of them have too complex structures etc. Finally, we focused on pyrazolopyrimidine skeleton (Figure 3) showing high in silico for antiandrogenic potency further study. The representative pyrazolopyrimidine analogs and corresponding predicted activities were listed in the Table 3.
Table 3. The details of docking score and calculated bioactivities using four models of representative pyrazolopyrimidine analogs.
Screening of Pyrazolopyrimidine Analogs
To find out more molecules with pyrazolopyrimidine skeleton, we searched the ZINC database, a free database of commercially-available compounds for virtual screening. On the ZINC website, the user can choose to search the database according to the structure similarity or substructure. In this study, our goal is to find out more molecules with pyrazolopyrimidine skeleton, so we used pyrazolopyrimidine as a substructure query to search ZINC database. Totally 5076 samples were found and downloaded to build a new chemical Database 4.
Then the four QSAR models were used again to calculate the in silico antiandrogenic potency of the chemicals contained in Database 4. We assumed that the chemicals with predictive value less than 7 were less bioactive, which were removed in the following study. As a result, 715 chemicals were remained for the following screening process. To find out compounds having high binding affinities with AR, the remained molecules were then imported to LigandFit and Libdock modules to do docking simulations. Taking into sufficient consideration of all the above results, we selected 45 molecules with high potency and submitted them to J&K Chemicals and SPECS. Finally we successfully purchased four compounds (Table 4), and the predicted biological activity from the Y1-Y4 models and the docking affinities were listed in Table S5.
In vitro Evaluation of Hit Compounds
In order to determine whether the screened hit compounds can affect the proliferation of prostate cancer cell, we carried out the cell proliferation assay by using the LNCaP cell line. The corresponding structures and biological activities (IC50) were listed in Table 4. It can be seen that all these compounds can suppress the proliferation of the LNCaP cells to varying degree, which suggests that pyrazolopyrimidine analogs can act as potential hits against CRPC. Compound IV, with the IC50 value of 45.8 μM, is 2-fold less active than bicalutamide. Compound III were demonstrated to be more effective than the first-generation AR antagonist bicalutamide and exhibited a dose dependent manner as shown in Figure 4.
To examine whether these two compounds have antagonistic activities toward AR, we performed transient transfection assay with androgen reporter pMMTV-Luc, which contains the natural AR target promoters, mouse mammary tumor virus (MMTV) long terminal repeat promoters, and luciferase gene bound at the downstream of an AR promoter. The concentration of DHT is 10 nM in all AR luciferase assays. As the result, compound IV have no significantly antagonistic activity against AR, while compound III can moderately inhibit the DHT-induced transcriptional activation of AR at the concentration of 10 μM (29.8%) (Figure 5). We assumed that pyrazolopyrimidine analogs with bulky groups may play a vital role in pushing away H12 to form open conformation (antagonistic form) and hence show the antagonistic bioactivity to AR, which could be further proved by molecular dynamics simulations.
Figure 5. AR antagonistic activity. The COS-7 cells were treated with positive compound Bic, compound III and IV (each of 10 μM) in the presence of 10 nM DHT. **p < 0.01 compared to cells treated only with DHT. All experiments were repeated at least three times.
Molecular Dynamics Simulation
In order to investigate how Compound III produces antagonistic effect on androgen receptor, molecular dynamics simulations and MM/GBSA methods were employed to study the interaction between Compound III and AR, including wild type AR, F876L and H874Y mutant types, which are commonly occurred in the AR LBD and often caused enzalutamide resistance problem (Liu et al., 2017). Then free energy calculation and decomposition analysis were carried out to analyze the MD results. Figure 6 shows the root mean square deviations (RMSDs) of all the backbone atoms of the protein during all 100 ns MD, from where it can be seen that after 80 ns, the RMSD values of the protein backbone atoms and binding pocket atoms have small fluctuations, so the last 20 ns trajectories were used to do all the analysis. We analyzed the results with earlier work aimed to enzalutamide and AR (Liu et al., 2017).
In Table 5, the binding free energies of Compound III/WT AR, Compound III/F876L AR, and Compound III/H874Y AR systems, are −44.37 kcal/mol, −57.04 kcal/mol and −42.02 kcal/mol, close to the enzalutamide' s data (Liu et al., 2017). It can be seen from Table 5 that van der Waal provides the main driving force for the binding between Compound III and WT/mutant ARs. Electrostatic energies are also very important for binding except Compound III/F876L AR complex. However, the polar solvation energies are unfavorable to the binding free energy.
Hydrogen bond interaction (Table 6) between Compound III and ARs show that all three complexes have relatively stable hydrogen bonds during last 20 ns in simulations. For the Compound III/WT AR system, Compound III formed two hydrogen bonds with T877 (98.86%) and Q711 (51.25%), while for the Compound III/F876L AR system, Compound III formed only one stable hydrogen bond with N705 (97.94%). Two hydrogen bonds are formed between R752 (85.41%), Q711 (65.50%) and Compound III respectively. Hydrogen bonds are beneficial to stabilize the position of the benzene of Compound III to prevent H12 from closing, as shown in Figure 7.
Figure 7. The superimposition of WT AR's initial structure with optimized structures of three systems in last 20 ns. Green cartoon and sticks represent the initial crystal structure. (A) The complexes of initial structure and WT (deep salmon cartoon and sticks). (B) The complexes of initial structure and Compound III/F876L AR (deep salmon cartoon and sticks). (C) The complexes of initial structure and Compound III/H874Y AR (deep salmon cartoon and sticks).
To get the optimized structure of each system, we performed a cluster analysis for the last 20 ns trajectories using Kclust algorithm (https://mmtsb.org/workshops/sean-bin_workshop_2012/Tutorials/MMTSB_EnsembleAnalysis/MMTSBEnsembleAnalysis.html). The systems of the Compound III/WT AR, Compound III/F876L AR and Compound III/H874Y AR were clustered into 6, 7, and 7 classes and the representative conformation accounted for 44.0, 29.7, and 35.0%, respectively. From the largest number of cluster, the conformation with the lowest RMSD to the cluster centers was selected. These optimized structures were used to analyze the conformation difference between corresponding initial structure and structure after 80 ns MD simulation, shown in Figure 6. We can see the superimposition of the wild type AR and the three optimized structures respectively in Figure 7. This figure shows that for all three systems, benzene ring of Compound III is close to the H12, meanwhile, the binding free energy decomposition analysis (Figure 8) also shows that residue L/F876 has weaker interactions with enzalutamide or Compound III than enza/F876L AR and enza/H874Y AR complexes. It has been experimentally confirmed (Korpal et al., 2013) or computational predicted (Liu et al., 2017) that F876L and H874Y mutation could switch enzalutamide from AR antagonist to AR agonist. The benzene ring linked to pyrazole of Compound III is close to H12 in all these systems, pushing H12 far away from the ligand binding pocket, which distorts the AF2 site resulting the inactivation of transcription. Based on the MD simulations results above, F876L and H874Y mutations in AR cannot convert compound III from AR antagonist to AR agonist.
Figure 8. Contribution of the important residues for ligand binding. All structures are average conformations generated from the last 20 ns snapshots of each MD system.
In this study, structure-based docking and ligand-based QSAR methodology were combined to virtually screen chemical databases to discover new androgen receptor antagonists with novel chemical skeleton. After virtual screening and the cell proliferation assay validation, the inhibition rate of prostate cancer cell LNCaP and DHT-induced transcriptional activation of AR were detected on a series of pyrazolopyrimidine analogs in vitro. Importantly, Compound III was proved to be the one of the most potent of these non-nitrophenyl and non-cyanophenyl type nonsteroidal AR antagonists and exhibited potent antiandrogenic activity. MD simulations and MM/GBSA methods were employed to investigate the antagonist mechanism of representative Compound III to AR, which further proved that this compound can prevent H12 in AR LBD from closing to distort the formation of AF2, resulting the invalid transcription. Hence, Pyrazolopyrimidine, as a novel antiandrogenic skeleton serves as a core structure of AR antagonist, deserves to be further optimized to discover more AR antagonists with high potency.
Materials and Methods
It was reported that virtual screening strategy combining structure- and ligand-based methods improves the performance and consistency of virtual screening compared to the single methods alone (Svensson et al., 2012). In the current study, molecular docking and QSAR models were used together to search a chemical database for potential androgen receptor antagonists.
Considering the diversity of the molecular structure, the commercially available ChemDiv structural library (http://www.chemdiv.com/) was used for screening, which includes 1.6 million compounds. All these molecular structures were protonated/deprotonated by generating ionization states with significant population at the specific 7.0 ± 2.0 in Maestro (Friesner et al., 2004). After adding partial charges, the molecular structures were minimized with MMFFs force field keeping its original stereochemistry.
Structure-Based Virtual Screening
The androgen receptor crystal structures were downloaded from Protein Data Bank (PDB code: 1T65) and prepared by using the Protein Preparation Wizard in Maestro. The treatments on AR structure include deleting the involved solvents, adding the missing hydrogen atoms (missing loops) and residues via Prime module, protonation protein crystal structure to neutral (pH = 7), minimizing side chains using OPLS 2005 force field, and definition receptor grid using a 12 Å box centered on the ligand. All other adjustable settings were set as default. Moreover, we constructed an antagonistic AR (named Anti-AR) model according to the reported homology modeling methodology (Liu et al., 2012).
Docking-Based Virtual Screening
First, all the chemical structures were collected and docked into the ligand binding domain of 1T65 and Anti-AR in Glide standard precision (SP) model, respectively. This program is generally used to approximate a systematic search of the conformational, oriented and positional space for the docked compound. The involved active site was defined from the coordination of the built-in ligand, using the default settings. The chemicals that could accommodate into the ligand binding domains of 1T65 and Anti-AR formed database 1. Then, the Lipinski's rule of five (Lipinski et al., 2001) was further utilized as a benchmark to filter the compounds to form database 2, which were submitted to the next process.
Ligand-Based Virtual Screening
Our group has published a predictive QSAR models (Y1) on androgen receptor antagonists (Wang et al., 2015), which can be used here to evaluate the in silico bioactivities of the screened chemicals. Furthermore, we found more AR antagonists with different structures from the literatures, and three new QSAR models (Y2, Y3, and Y4) were established and used in the present study. The processes to build models Y2, Y3, and Y4 are as follows.
The chemical structures together with corresponding bioactivities taken from literatures (Hamann et al., 1998; Zhi et al., 1998, 1999; Kong et al., 2000; Tachibana et al., 2007a,b; Zhao et al., 2008) were collected to develop models Y2 − Y4 respectively. The biological activities (IC50) were converted into negative logarithmic unit pIC50 and used as dependent variables. The 3D structures of these molecules were imported into SYBYL 6.9 (SYBYL, version 6.9, Tripos Inc., St. Louis, MO)2 to calculate Gasteiger-Hückel partial charges. The energy minimization was performed using the Tripos force field with convergence criterion of 0.01 kcal/mol·Å. Then the multi-search routine was performed to do conformation search to obtain the lowest energy conformation for each molecule. Subsequently, the compounds were randomly split into training sets and prediction sets.
The obtained lowest energy confirmations were submitted to DRAGON 5.5 (DRAGON for Windows, version 5.5)3 to calculate all the 20 types of descriptors. To decrease redundant information, the calculated descriptors were pretreated to exclude the constant or near constant value (variance < 1) and high correlated pairwise (R > 0.99).
Modeling Method and Validation
Here in this study, genetic algorithm was employed to select important descriptors highly related to the bioactivities, using the MLR method to build the QSAR models, executed in QSARINS (Gramatica et al., 2013). To thoroughly validate the built QSAR models, several internal and external validation standards were combined to evaluate the robustness and predictive ability of the built models (Gramatica and Sangion, 2016). After strict validation, these models were used to predict the bioactivities of the screened chemicals. The chemicals with high in silico bioactivities formed database 3.
After the structure- and ligand-based virtual screening, chemicals with high binding affinities and high in silico bioactivities were filtered out for the succeeding analysis. Molecules with similar skeleton were extracted.
Then, according to our experience on medicinal chemistry and combining with the docking and QSAR results, the novel skeleton of pyrazolopyrimidine was eventually proposed for further research. To thoroughly explore the pyrazolopyrimidine analogs, the specific skeleton was used as a substructure to search the ZINC online database (http://zinc15.docking.org/), and the resulted chemicals were downloaded to form a new pyrazolopyrimidine database 4. Subsequently, LigandFit and Libdock modules were used to dock all these pyrazolopyrimidine analogs to the ligand binding domain of 1T65 and Anti-AR protein, and the four QSAR models were employed to predict corresponding bioactivities. At the end, the hit chemicals with excellent performance both in the docking and QSAR model predictions were selected for further experimental validation.
In vitro Assay and Biological Evaluation
The selected chemicals were purchased from the established suppliers, including the J&K Chemicals and SPECS.
Cell Proliferation Assay
The LNCaP cells were cultured in RPMI 1640 medium containing 10% FBS, 100 unit/ml penicillin and 100 unit/ml streptomycin at 37°C in a humidified atmosphere with 5% CO2. Before being treated with compounds, cells were cultured in 100 mm cell dish at a density of 5 × 105 per well, then sub-cultured when it grew to approximate 90% confluence. Subsequently, the compounds were evaluated in a cellular assay by measuring IC50 values. All experiments were performed during the logarithmic phase of cell growth. For analysis of the effect of all chemicals on LNCaP cell proliferation, the cells were seeded in 96 well plates at a density of 5,000 LNCaP cells per well and treated with chemicals for 120 h. A 20 μL aliquot of tetrazole (MTT) was added to each well. After 4 h incubation, 150 μL of isopropanol was added to dissolve the crystal. Finally, the absorbance was measured at a wavelength 570 nm.
AR Reporter Gene Assay
Luciferase reporter assay based on androgen response elements (ARE) was used to examine chemicals for androgenic activity. COS-7 cells were seeded in a 24 well plate, and cultured in medium (DMEM medium containing 10% charcoal-stripped FBS (CSS), 2 mM glutamine) for 24 h in a 24 well plate. The plasmids of pcDNA3.1-AR, pMMTV-Luc and pRL-SV40 were transfected in COS-7 cells by using transfection reagent Lipofectamine 3000 according to the instructions of manufacturer. After culturing at 37°C in a 5% CO2 atmosphere for 24 h, cells were treated with the chemicals (final concentration, 10 μM) in the presence of 10 nM DHT for 24 h. Then, cells were harvested with 150 μL of cell passive lysis buffer (Promega). The firefly and renilla luciferase activities were determined with a Dual-Glo Luciferase Assay Kit (Promega). The data were obtained in triplicate and expressed as inhibition rate over the DHT control. Inhibition% = 1−(RLUtest − RLUblank)/(RLUDHT − RLUblank) × 100%. RLU = relative light unit.
Molecular Dynamics Simulations
The crystal structures of WT, H874Y mutant AR (PDB ID 2Q7L) were obtained from the Protein Data Bank (http://www.rcsb.org/pdb)4. Then Pymol (v1.7) program was used to generate the 3D structure of F876L mutant AR by mutated wild type AR 876F to 876L. Then the CDOCKER module of Discovery Studio 2.5 [(Discovery Studio version 2.5, 2009) Accelrys Inc. CA] was used to dock Compound III to all the three kinds of ARs. The lowest-energy models were selected as the docking results.
MD simulations were performed with AMBER12 (Case et al., 2012) package. The geometry optimization and partial charges calculation of small-molecule ligand were performed in Gaussian09 program using HF/6-31G*(Frisch et al., 2009) basis set. Then the restrained electrostatic protential (RESP) (Bayly et al., 1993; Cieplak et al., 1995; Fox and Kollman, 1998) of Compound III was assigned using the general AMBER force filed (GAFF) (Wang et al., 2004). The ff99SB force field (Hornak et al., 2006) was used to parameterize the systems. Finally, all the systems were neutralized with chloride ions and immersed in a cubic TIP3P water box (Ryckaert et al., 1977). The solvent boundary was at least 10Å away from the proteins. Subsequently, all systems were minimized using a steepest decent method followed by conjugate gradient method, and then heated from 0 to 310 K. The temperature was controlled by Langevin thermostat with a coupling coefficient of 2.0 ps−1. SHAKE algorithm (Eaamann et al., 1995) was used to constrained bond lengths involving hydrogen atoms. All equilibration and subsequent MD stages were carried out in the isothermal isobaric (NTP) ensemble with a target pressure of 1 bar. The time step was set as 2 fs and the production time for all the systems was 100 ns.
The binding free energy of each system was calculated using molecular mechanics/generalized born surface area (MM/GBSA) methodology (Massova and Kollman, 2000; Xu et al., 2013; Sun et al., 2014a,b). The binding free energy (ΔGbind) is consisted of electrostatic interaction energy (ΔEele), van der Waals interaction energy (ΔEvdw), polar solvation free energy (ΔGp) and nonpolar solvation free energy (ΔGnp). Dielectric constants of 1.0 and 80.0 were set for solute and solvent respectively. The polar part of desolvation (ΔGGB) was computed with the modified GB model developed by Onufriev et al. (referred to as igb = 2 in Amber) (Onufriev et al., 2004), and the nonpolar part of desolvation (ΔGSA) was determined based on the solvent accessible surface area (SASA) computed by the LCPO algorithm: (Weiser et al., 1999) ΔGSA = 0.0072ΔSASA. The change of the conformational entropy (−TΔS) was not considered due to high computational cost and low prediction accuracy (Wang et al., 2006; Hou et al., 2011). In total, 2,000 snapshots evenly extracted from 80 to 100 ns were used to calculate the energy terms. For each systems, the interaction spectrum between Compound III and ARs on a per-residue basis was calculated by MM/GBSA decomposition analysis supported by the mmpbsa.py module in AMBER (Gohlke et al., 2003; Hou et al., 2008, 2012).
XW and JL conceived and coordinated the study. LW did the structure-based virtual screening, screening of pyrazolopyrimidine analogs, in vitro evaluation of hit compounds, TS did the molecular docking, molecular dynamics simulations, and they wrote the paper. All authors analyzed the results and approved the final version of the manuscript.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This research was supported by National Natural Science Foundation of China (81403145) and Fundamental Research Funds for the Central Universities (lzujbky-2018-136).
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2018.00864/full#supplementary-material
Molecular structures and corresponding biological activities used in this study to build QSAR models Y2-Y4 were listed. Definitions of external validation metrics were listed in Table S4. The predicted biological activity from the Y1-Y4 models and the docking affinities were listed in Table S5.
2. ^SYBYL version 6.9. Tripos Inc., St. Louis, M. O.
Bayly, C. I., Cieplak, P., Cornell, W. D., and Kollman, P. A. (1993). A well behaved electrostatic potential based method using charge restraints for deriving atomic charges—the RESP model. J. Phys. Chem. B 97, 10269–10280. doi: 10.1021/j100142a004
Bohl, C. E., Gao, W., Miller, D. D., Bell, C. E., and Dalton, J. T. (2005a). Structural basis for antagonism and resistance of bicalutamide in prostate cancer. Proc. Natl. Acad. Sci. U.S.A. 102, 6201–6206. doi: 10.1073/pnas.0500381102
Bohl, C. E., Miller, D. D., Chen, J., Bell, C. E., and Dalton, J. T. (2005b). Structural basis for accommodation of nonsteroidal ligands in the androgen receptor. J. Biol. Chem. 280, 37747–37754. doi: 10.1074/jbc.M507464200
Chirico, N., and Gramatica, P. (2011). Real external predictivity of QSAR models: how to evaluate it? Comparison of different validation criteria and proposal of using the concordance correlation coefficient. J. Chem. Inf. Model. 51, 2320–2335. doi: 10.1021/ci200211n
Chirico, N., and Gramatica, P. (2012). Real external predictivity of QSAR models. Part 2. New intercomparable thresholds for different validation criteria and the need for scatter plot inspection. J. Chem. Inf. Model. 52, 2044–2058. doi: 10.1021/ci300084j
Cieplak, P., Cornell, W. D., Bayly, C., and AKollman, P. (1995). Application of the multimolecule and multiconformational RESP methodology to biopolymers: charge derivation for DNA, RNA, and proteins. J. Comput. Chem. 16, 1357–1377. doi: 10.1002/jcc.540161106
Denmeade, S. R., Jakobsen, C. M., Janssen, S., Khan, S. R., Garrett, E. S., Lilja, H., et al. (2003). Prostate-specific antigen-activated thapsigargin prodrug as targeted therapy for prostate cancer. J. Natl. Cancer Inst. 95, 990–1110. doi: 10.1093/jnci/95.13.990
Eaamann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). Properties of organic liquids when simulated with simulated with long-range Lennard-Jones interactions. J. Chem. Phys. 103, 8577–8593.
Fix, C., Jordan, C., Cano, P., and Walker, W. H. (2004). Testosterone activates mitogen-activated protein kinase and the cAMP response element binding protein transcription factor in Sertoli cells. Proc. Natl. Acad. Sci. U.S.A. 101, 10919–10924. doi: 10.1073/pnas.0404278101
Friesner, R. A., Banks, J. L., Murphy, R. B., Halgren, T. A., Klicic, J. J., Mainz, D. T., et al. (2004). Glide: a new approach for rapid, accurate docking and scoring. 1. method and assessment of docking accuracy. J. Med. Chem. 47, 1739–1749. doi: 10.1021/jm0306430
Gohlke, H., Kiel, C., and Case, D. A. (2003). Insights into protein-protein binding by binding free energy calculation and free energy decomposition for the Ras-Raf and Ras-RaIGDS complexes. J. Mol. Biol. 330, 891–913. doi: 10.1016/S0022-2836(03)00610-7
Gramatica, P., Chirico, N., Papa, E., Cassani, S., and Kovarich, S. (2013). QSARINS: a new software for the development, analysis, and validation of QSAR MLR models. J. Comput. Chem. 34, 2121–2132. doi: 10.1002/jcc.23361
Gramatica, P., and Sangion, A. (2016). A historical excursus on the statistical validation parameters for QSAR models: a clarification concerning metrics and terminology. J. Chem. Inf. Model. 56, 1127–1131. doi: 10.1021/acs.jcim.6b00088
Hamann, L. G., Higuchi, R. I., Zhi, L., Edwards, J. P., Wang, X. N., Marschke, K. B., et al. (1998). Synthesis and biological activity of a novel series of nonsteroidal, peripherally selective androgen receptor antagonists derived from 1,2-dihydropyridono[5,6-g] quinolones. J. Med. Chem. 41, 623–639. doi: 10.1021/jm970699s
Hodgson, M. C., Astapoval, I., Hollenberg, A. N., and Balk, S. P. (2007). Activity of androgen receptor antagonist bicalutamide in prostate cancer cells is independent of NCoR and SMRT corepressors. Cancer Res. 67, 8388–8395. doi: 10.1158/0008-5472.CAN-07-0617
Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C. (2006). Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 65, 712–725. doi: 10.1002/prot.21123
Hou, T. J., Li, N., Li, Y. Y., and Wang, W. (2012). Characterization of domain-peptide interaction interface: prediction of SH3 domain-mediated protein-protein interaction network in yeast by generic structure-based models. J. Proteome Res. 11, 2982–2995. doi: 10.1021/pr3000688
Hou, T. J., Zhang, W., Case, D. A., and Wang, W. (2008). Characterization of domain-peptide interaction interface: a case study on the amphiphysin-1 SH3 domain. J. Mol. Biol. 376, 1201–1214. doi: 10.1016/j.jmb.2007.12.054
Kong, J. W., Hamann, L. G., Ruppar, D. A., Edwards, J. P., Marschke, K. B., and Jones, T. K. (2000). Effects of isosteric pyridone replacements in androgen receptor antagonists based on 1,2-dihydro-and 1,2,3,4-tetrahydro-2,2-dimethyl-6 -trifluoromethyl-8-pyridono[5,6-g] quinolones. Bioorg. Med. Chem. Lett. 10, 411–414. doi: 10.1016/S0960-894X(00)00010-X
Korpal, M., Korn, J. M., Gao, X. L., Rakiec, D. P., Ruddy, D. A., Doshi, S., et al. (2013). An F876L mutation in androgen receptor confers genetic and phenotypic resistance to MDV3100 (Enzalutamide). Cancer Discov. 3, 1030–1043. doi: 10.1158/2159-8290.CD-13-0142
Lipinski, C. A., Lombardo, F., Dominy, B. W., and Feeney, P. J. (2001). Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv. Drug Deliv. Rev. 46, 3–26. doi: 10.1016/S0169-409X(00)00129-0
Liu, B., Geng, G., Lin, R., Ren, C., and Wu, J. H. (2012). Learning from estrogen receptor antagonism: structure-based identification of novel antiandrogens effective against multiple clinically relevant androgen receptor mutants. Chem. Biol. Drug Des. 79, 300–312. doi: 10.1111/j.1747-0285.2011.01290.x
Liu, H. L., Wang, L. Y., Tian, J. Q., and Li, J. Z. (2017). Molecular dynamics studies on the enzalutamide resistance mechanisms induced by androgen receptor mutations. J. Cell. Biochem. 118, 2792–2801. doi: 10.1002/jcb.25928
Massova, I., and Kollman, P. A. (2000). Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding. Perspect. Drug. Discov. 18, 113–135. doi: 10.1023/A:1008763014207
Onufriev, A., Bashford, D., and Case, D. A. (2004). Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 55, 383–394. doi: 10.1002/prot0.20033
Ryckaert, J. P., Ciccotti, G., and Berendsen, H. J. C. (1977). Numerical integration of cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. doi: 10.1016/0021-9991(77)90098-5
Schüürmann, G., Ebert, R. U., Chen, J., Wang, B., and Kühne, R. (2008). External validation and prediction employing the predictive squared correlation coefficient test set activity mean vs training set activity mean. J. Chem. Inf. Model. 48, 2140–2145. doi: 10.1021/ci800253u
Sun, H., Li, Y., Shen, M., Tian, S., Xu, L., Pan, P., et al. (2014a). Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys. 16, 22035–22045. doi: 10.1039/c4cp03179b
Sun, H., Li, Y., Tian, S., Xu, L., and Hou, T. J. (2014b). Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys. Chem. Chem. Phys. 16, 16719–16729. doi: 10.1039/c4cp01388c
Tachibana, K., Imaoka, I., Yoshino, H., Emura, T., Kodama, H., Furuta, Y., et al. (2007a). Discovery of 7alpha-substituted dihydrotestosterones as androgen receptor pure antagonists and their structure-activity relationships. Bioorg. Med. Chem. 15, 174–185. doi: 10.1016/j.bmc.2006.09.072
Tachibana, K., Imaoka, I., Yoshino, H., Kato, N., Nakamura, N., Ohta, M., et al. (2007b). Discovery and structure–activity relationships of new steroidal compounds bearing a carboxy-terminal side chain as androgen receptor pure antagonists. Bioorg. Med. Chem. Lett. 17, 5573–5576. doi: 10.1016/j.bmcl.2007.07.090
Taplin, M. E., Rajeshkumar, B., Halabi, S., Werner, C. P., Woda, B. A., Picus, J., et al. (2003). Androgen receptor mutations in androgen-independent prostate cancer: cancer and leukemia group B study 9663. J. Clin. Oncol. 21, 2673–2678. doi: 10.1200/JCO.2003.11.102
Tilley, W. D., Lim-Tio, S. S., Horsfall, D. J., Aspinall, J. O., Marshall, V. R., and Skinner, J. M. (1994). Detection of discrete androgen receptor epitopes in prostate cancer by immunostaining: measurement by color video image analysis. Cancer Res. 54, 4096–4102.
Tran, C., Ouk, S., Clegg, N. J., Chen, Y., Watson, P. A., Arora, V., et al. (2009). Development of a second-generation antiandrogen for treatment of advanced prostate cancer. Science 324, 787–790. doi: 10.1126/science.1168175
Wang, J. M., Hou, T. J., and Xu, X. J. (2006). Recent advances in free energy calculations with a combination of molecular mechanics and continuum models. Curr. Comput. Aid. Drug Des. 2, 287–306. doi: 10.2174/157340906778226454
Wang, Y. W., Bai, F., Cao, H., Li, J. Z., Liu, H. X., and Gramatica, P. (2015). A combined quantitative structure-activity relationship research of quinolinone derivatives as androgen receptor antagonists. Comb. Chem. High. T. Scr. 18, 834–845. doi: 10.2174/1386207318666150831125750
Weiser, J., Shenkin, P. S., and Still, W. C. (1999). Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 20, 217–230. doi: 10.1002/(SICI)1096-987X(19990130)20:2<217::AID-JCC4>3.0.CO;2-A
Xu, L., Sun, H., Li, Y., Wang, J., and Hou, T. J. (2013). Assessing the performance of MM/PBSA and MM/GBSA methods. 3. The impact of force fields and ligand charge models. J. Phys. Chem. B 117, 8408–8421. doi: 10.1021/jp404160y
Zhao, S., Shen, Y., van Oeveren, A., Marschke, K. B., and Zhi, L. (2008). Discovery of a novel series of nonsteroidal androgen receptor modulators: 5- or 6-oxachrysen-2-ones. Bioorg. Med. Chem. Lett. 18, 3431–3435. doi: 10.1016/j.bmcl.2008.03.085
Zhi, L., Tegley, C. M., Kallel, E. A., Marschke, K. B., Mais, D. E., Gottardis, M. M., et al. (1998). 5-Aryl-1,2-dihydrochromeno[3,4-f] quinolines: a novel class of nonsteroidal human progesterone receptor agonists. J. Med. Chem. 41, 291–302. doi: 10.1021/jm9705768
Zhi, L., Tegley, C. M., Marschke, K. B., and Jones, T. K. (1999). Switching androgen receptor antagonists to agonists by modifying C-ring substituents on piperidino[3,2-g] quinolinone. Bioorg. Med. Chem. Lett. 9, 1009–1012. doi: 10.1016/S0960-894X(99)00119-5
Keywords: prostate cancer, androgen receptor, virtual screening, pyrazolopyrimidine, AR reporter gene assay, molecular dynamics
Citation: Wang L, Song T, Wang X and Li J (2018) Discovery and Identification of Pyrazolopyramidine Analogs as Novel Potent Androgen Receptor Antagonists. Front. Pharmacol. 9:864. doi: 10.3389/fphar.2018.00864
Received: 09 May 2018; Accepted: 17 July 2018;
Published: 28 August 2018.
Edited by:Salvatore Salomone, Università degli Studi di Catania, Italy
Reviewed by:Elias Castanas, University of Crete, Greece
Greta Varchi, Consiglio Nazionale delle Ricerche (Bologna), Italy
Copyright © 2018 Wang, Song, Wang and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.