Discovery and Identification of Pyrazolopyramidine Analogs as Novel Potent Androgen Receptor Antagonists

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.

Introduction

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

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.

FIGURE 2

Figure 2. The results of screening process step by step.

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 Y₁ was previously reported from our group (Wang et al., 2015). Other three models Y₂, Y₃, and Y₄ were newly built. Experimental and predicted antiandrogenic activities in these three models are listed in Supporting Information Tables S1–S3. Here the antagonistic activities (IC₅₀) 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, pIC₅₀. 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.

Y₁ = 1.188GATS7v − 1.108BEHp7 − 7.281E2u + 1.453HATS4u − 2.647H6m + 28.810R6u+ − 46.815R8u+ + 12.157

Y₂ = −2.89IC5 + 1.01GATS5e − 3.17DISPp − 12.99HATS3u + 27.13

Y₃ = 5.601IVDE − 1.70C-009 + 1.11BLTF96 + 3.02

Y₄ = 3.18 IC1 − 0.46F05[N-F] + 1.47R4u + 0.409 Depressant-80 + 3.85HATS7m − 7.15

TABLE 1

Table 1. The QSAR MLR models and corresponding statistical parameters.

TABLE 2

Table 2. The descriptors used to build the four linear model and corresponding meanings.

In Table 1, the parameter R² describes the fitting ability and $Q_{LOO}^2$ represents the reliability and stability of a QSAR model. Leave-one-out cross-validated correlation coefficient ($Q_{LOO}^2$) 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 R² greater than 0.75. The predictive abilities were evaluated by means of four external validation parameters (Gramatica and Sangion, 2016), $Q_{F1}^2$ (www.oecd.org/dataoecd/33/37/37849783.pdf), $Q_{F2}^2$ (Schüürmann et al., 2008), $Q_{F3}^3$ (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)¹ program to calculate all molecular descriptors. Then the descriptors involved in these four QSAR models were extracted and submitted to the equations Y₁-Y₄ 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.

FIGURE 3

Figure 3. The skeleton structure of pyrazolopyrimidine analogs.

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.

TABLE 4

Table 4. Structures and experimental activities of identified active compounds.

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 (IC₅₀) 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 IC₅₀ 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.

FIGURE 4

Figure 4. Inhibition of the proliferation for LNCaP cells by compound III, IV and bicalutamide.

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

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).

FIGURE 6

Figure 6. The RMSDs of backbone atoms of the AR in all three trajectories.

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.

TABLE 5

Table 5. The calculated binding free energies (kcal/mol) of the three systems by MM/GBSA method.

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.

TABLE 6

Table 6. The hydrogen bonds between compound III and key residues in ARs.

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

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.

Conclusion

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.

Database Preparation

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

Protein Preparation

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 (Y₁) 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 (Y₂, Y₃, and Y₄) were established and used in the present study. The processes to build models Y₂, Y₃, and Y₄ are as follows.

Data Set

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 Y₂ − Y₄ respectively. The biological activities (IC₅₀) were converted into negative logarithmic unit pIC₅₀ 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)² 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.

Descriptors Generation

The obtained lowest energy confirmations were submitted to DRAGON 5.5 (DRAGON for Windows, version 5.5)³ 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.

Structural Analysis

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

Materials Preparation

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% CO₂. Before being treated with compounds, cells were cultured in 100 mm cell dish at a density of 5 × 10⁵ per well, then sub-cultured when it grew to approximate 90% confluence. Subsequently, the compounds were evaluated in a cellular assay by measuring IC₅₀ 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% CO₂ 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−(RLU_test − RLU_blank)/(RLU_DHT − RLU_blank) × 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)⁴. 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⁻¹. 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 (ΔG_bind) is consisted of electrostatic interaction energy (ΔE_ele), van der Waals interaction energy (ΔE_vdw), polar solvation free energy (ΔG_p) and nonpolar solvation free energy (ΔG_np). Dielectric constants of 1.0 and 80.0 were set for solute and solvent respectively. The polar part of desolvation (ΔG_GB) 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 (ΔG_SA) was determined based on the solvent accessible surface area (SASA) computed by the LCPO algorithm: (Weiser et al., 1999) ΔG_SA = 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).

Author Contributions

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.

Acknowledgments

This research was supported by National Natural Science Foundation of China (81403145) and Fundamental Research Funds for the Central Universities (lzujbky-2018-136).

Supplementary Material

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.

Footnotes

1. ^DRAGON for Windows(Software for Molecular Descriptor Calculations), version 5.5 - Talet srl.http://www.talete.mi.it/

2. ^SYBYL version 6.9. Tripos Inc., St. Louis, M. O.

3. ^DRAGON for Windows(Software for Molecular Descriptor Calculations), version 5.5 - Talet srl. http://www.talete.mi.it/

4. ^RCSB Protein Data Bank. http://www.rcsb.org/pdb

References

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Darden, T. A., Cheatham, T. E., Simmerling, C. L., Wang, J., Duke, R. E., et al. (2012). AMBER 12. San Francisco, CA: University of California.

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Consonni, V., Ballabio, D., and Todeschini, R. (2009). Comments on the definition of the Q2 parameter for QSAR validation. J. Chem. Inf. Model. 49, 1669–1678. doi: 10.1021/ci900115y

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Discovery Studio version 2.5 (2009). Discovery Studio Version 2.5. California, CA: Accelrys Inc.

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.

Evans, R. M. (1988). The steroid and thyroid hormone receptor superfamily. Science 240, 889–895. doi: 10.1126/science.3283939

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Fox, T., and Kollman, P. A. (1998). Application of the RESP methodology in the parametrization of organic solvents. J. Phys. Chem. B 102, 8070–8079. doi: 10.1021/jp9717655

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2009). Gaussian 09. Wallingford, CT: Gaussian Inc.

Gao, W., Bohl, C. E., and Dalton, J. T. (2005). Chemistry and structural biology of and rogen receptor. Chem. Rev. 105, 3352–3370. doi: 10.1021/cr020456u

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, T. J., Li, Y. Y., and Wang, W. (2011). Prediction of peptides binding to the PKA RII alpha subunit using a hierarchical strategy. Bioinformatics 27, 1814–1821. doi: 10.1093/bioinformatics/btr294

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Jernberg, E., Bergh, A., and Wikström, P. (2017). Clinical relevance of androgen receptor alterations in prostate cancer. Endocr. Connect. 6, R146–R161. doi: 10.1530/EC-17-0118

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, L. I. (1989). A concordance correlation coefficient to evaluate reproducibility. Biometrics 45, 255–268. doi: 10.2307/2532051

PubMed Abstract

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Sadar, M. D. (2012). Advances in small molecule inhibitors of androgen receptor for the treatment of advanced prostate cancer. World J. Urol. 30, 311–318. doi: 10.1007/s00345-011-0745-5

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2016). Cancer statistics, 2016. CA. Cancer J. Clin. 66, 7–30. doi: 10.3322/caac.21332

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Svensson, F., Karlén, A., and Sköld, C. (2012). Virtual screening data fusion using both structure- and ligand-based methods. J. Chem. Inf. Model. 52, 225–232. doi: 10.1021/ci2004835

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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.

PubMed Abstract | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Wang, J. M., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi: 10.1002/jcc.20035

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar