In Silico Docking of Rhodanine Derivatives and 3D-QSAR Study to Identify Potent Prostate Cancer Inhibitors

Koffi Alexis Respect Kouassi¹, Adenidji Ganiyou^2*, Diomande Gbe Gondo Didier², Anoubilé Benié^1,3, Ziao Nahossé¹

¹Laboratoire de Thermodynamique et de Physico-Chimie du Milieu, UFR SFA, Université Nangui Abrogoua, Abidjan, Côte d’Ivoire.
²Laboratoire de Constitution et de Réaction de la matière, UFR-SSMT, Université Félix Houphouët-Boigny, Abidjan, Côte d’Ivoire.
³Laboratoire de Chimie Bio-Organique et de Substances Naturelles, UFR-SFA, Université Nangui Abrogoua, Abidjan, Côte d’Ivoire.

DOI: 10.4236/cc.2022.102002                

Abstract

The 3VHE protein is considered as a potential target for the treatment of prostate cancer. In order to find new 3VHE inhibitors, pharmacophore models based on the molecular structure of rhodanine derivatives and a three-dimensional quantitative structure-activity relationship model (3D-QSAR) have been developed and validated by different methods. The 3D-QSAR model was evaluated for its predictive performance on a diverse test set containing 18 prostate cancer inhibitors. It presents very interesting internal and external statistical validation parameters (SD = 0.081; R2 = 0.903; Q2 = 0.869; ; F = 247.2). This result suggests that the 3D-QSAR combinatorial model can be used to search for new 3VHE inhibitors and predict their potential activity. Based on the combinatorial pharmacophore model, a virtual screening of the Enamine database was performed. Compounds selected after virtual screening were subjected to molecular docking protocols (HTVS, SP, XP and IFD). Twenty new active compounds have been identified and their absorption, distribution, metabolism and excretion (ADME) property calculated using Schr?dinger’s Qikprop module. These results suggest that these new compounds could constitute new chemical starting points for further structural optimization of 3VHE inhibitors.

Keywords

Rhodanine, Pharmacophores, 3D-QSAR, Molecular Docking

1. Introduction

Prostate cancer accounts for 40% of all cancers affecting men. The incidence and mortality of prostate cancer continue to rise and is a major health problem in the world. It is the most common newly diagnosed cancer in men and has a very high death rate. This cancer death rate in men is second behind lung cancer [1]. Many factors of genetic, toxicological and diet-related origin seem to be involved in the development of this cancer. The treatment suggested to the patient depends on the size of the tumor in cancer. Doctors also take into account the evolutionary nature of it and the general state of health of the person. The treatment depends on the stage of advancement. The main treatment methods for prostate cancer are surgery, radiotherapy (external radiotherapy), active monitoring (which makes it possible to postpone the start of treatment), hormone therapy and chemotherapy. However, fundamental questions remain about the prevention and treatment of prostate cancer. Treatment for metastatic disease is strictly palliative and there is still no treatment for the disease. The thiazolidin-2,4-diones, rhodanine derivatives, five-membered heterocyclic compounds, have a wide range of biological activities that include antioxidants [2], anti-in-flammatory [3] [4], antibacterial [5] [6], antifungals [7] [8], and the most important being the anticancer activity [9] [10] [11]. The presence of rhodanine in a very wide range of compounds possessing very varied biological properties makes it an important compound in the search for new drugs. Computer-aided design methods which include a series of techniques for discovering, designing and improving chemicals in silico are of considerable support in this research. [12] [13]. Five-membered heterocycles, which exhibit interesting biological activities as well as important industrial applications, are such potential targets [14]. This is the context in which our work, which aims to identify, using in silico methods, new prostate cancer inhibitors with rhodanine and its derivatives as base molecules.

2. Materials and Methods

2.1. Selection of Biological Dataset

A data set of 74 rhodanine derivatives with their anticancer activity IC50 (µM) as inhibitors of prostate cancer (PC3) was extracted from the work of Coulibaly et al. [15]. All molecular structures and activity data used for pharmacophore modeling, 3D-QSAR study, virtual screening and molecular docking are shown in Table 1. For these studies, the inhibitory activities (IC50 values in M) for each compound were changed to a negative logarithm of IC₅₀ (pIC₅₀). The pIC50 values were used as dependent variable for the development of the model in the 3D QSAR. All compounds have a similar structure and bioassay method.

2.2. Ligand Preparation

The 3D structures of the ligands were generated using the construction panel in Maestro and optimized using the LigPrep module [16]. Partial atomic charges

Table 1. Structures and pIC50 of PC3 inhibitors.

were assigned and possible ionization states were generated at a pH of 7.0 ± 2.0. The OPLS_2005 force field was used for optimization of the production of the low energy ligand conformer [17]. Energy minimization was performed for each ligand until it reached a root mean square deviation threshold of 0.01 Å.

2.3. Pharmacophore Model Generation

Schrödinger’s Phase module for ligand-based drug design was used to develop pharmacophore hypotheses [18]. Before generating the pharmacophoric model, Phase first recruits the Glide XP tool to anchor the ligand already bound to the active site of the protein under study and selects the best-ranked pose for the generation of the pharmacophores model. The chemical characteristics of all ligands were defined by six characteristics of pharmacophores: H-bond acceptor (A), H-bond donor (D), hydrophobic group (H), negatively charged group (N), positively charged group (P) and aromatic ring (R). An active analogue approach was used to identify common pharmacophore hypotheses, in which common pharmacophores were selected from the conformations of the active ligand set using a hierarchical partitioning technique that groups similar pharmacophores according to their inter-site distances [19]. The resulting pharmacophores were then noted and classified. Scoring was done to identify the best candidate hypothesis, which provided an overall ranking of all hypotheses. The scoring algorithm included contributions from site point and vector alignment, volume overlap, selectivity, number of paired ligands, relative conformational energy, and biological activity [19].

2.4. Pharmacophore Validation

The validation of a pharmacophore model is a fundamental first step that must be performed to prove its accuracy and specificity in the selection of active molecules, while guiding the virtual screening of ligands from a database. In the present study, a set of 324 decoy molecules extracted from Directory of Useful Decoys (http://dude.docking.org/) [20] [21] enriched with 9 active molecules was used. Before validation, the preprocessing of the active and decoy datasets was performed using the LigPrep module. All possible ionizable states as well as tautomeric forms at a pH range of 7.0 ± 2.0 were generated using LigPrep module [16]. For each compound, at most 32 conformers were generated by default and low energy stereoisomers with correct chirality were engaged for further analysis. The hypothesis validation tool of the Phase module [22] has been used. This tool uses the hypothesis file and the decoy and active dataset as input to calculate the performance parameters mentioned above. Various statistical parameters including Enrichment Factors (EF), Robust Initial Improvement (RIE), Boltzmann Enhanced Discrimination of Receiver Operating Characteristic (BEDROC), Area Under Accumulation Curve (AUC) and receiver operating characteristics (ROC) were calculated to validate the hypothesis [23].

2.5. Guner-Henry Score Validation

Güner-Henry scoring method is applied for quantification of model selectivity and evaluation of model effectiveness of similarity search. This scoring evokes the actives from a molecule dataset consisting of known active and inactive molecules. This scoring system ranges from 0 to 1, where 0 specifies a null model and 1 specifies an ideal model. The score is expected to be greater than 0.7 [24]. The formulas used for calculating GH score are given below:

$GH=\left(\frac{H_a\left(3A+H_t\right)}{4H_{t}A}\right)\left(1-\frac{H_t-H_a}{D-A}\right)$

$\%A=\frac{H_a}{A}\times 100;\%Y=\frac{H_a}{H_t}\times 100;EF=\frac{H_a/H_t}{A/D}$

where H_a is the number of actives in the hits list (true positives), A is the number of active compounds in the database, H_t is the number of hits retrieved, D is the number of compounds in the database, %A is the percentage of known active compounds obtained from the database, %Y is the percentage of known actives in the hits list, EF is the enrichment of the concentration of actives by the model relative to random screening without a pharmacophoric approach. Güner-Henry score is considered as a relevant metric, it takes into account both percent yield of actives in a database (%Y) and the percent ratio of actives in the hit list (%A).

2.6. 3D-QSAR

QSAR modeling was carried out by dividing the data set (74 molecules) into a training set (75%) and a test set (25%) randomly. PHASE presents two options for the alignment of the 3D structure of molecules: the pharmacophore-based alignment and the atom-based alignment [19] [25]. In this study, we used an atom-based QSAR model, which is more useful in explaining the structure-activity relationship. In atom-based QSAR, a molecule is treated as a set of overlapping van der Waal spheres. Each atom (and therefore each sphere) is classified into one of six categories according to a set of simple rules: hydrogens attached to polar atoms are classified as donors of hydrogen bonds (D); C-H carbons, halogens and hydrogen are classified as hydrophobic/nonpolar (H); atoms with an explicit negative ionic charge are classified as negative ionic (N); atoms with an explicit positive ionic charge are classified as positive ionic (P); nonionic nitrogen and oxygen are classified as electron attractors (W); and all other types of atoms are classified as miscellaneous (X).

For the purposes of QSAR development, the van der Waal models of the aligned molecules of training set were placed in a regular grid of cubes, each cube being assigned zero or more “bits” to account for the different types of atoms in learning set that occupy the cube. This representation results in binary-valued occupancy models that can be used as independent variables to create partial least squares (PLS) QSAR models. A five-component model (PLS factor) with good statistics was obtained for the dataset. The statistical quality of the generated QSAR models was judged by parameters such as the regression coefficient (R²), cross validation ( $Q_{CV}^2$ ), the variance (F), the confidence interval (P), the mean squared error (RMSE) and Pearson’s coefficient r [26].

2.7. High Throughput Virtual Screening and Molecular Docking

The molecules obtained after pharmacophore screening were filtered by High Throughput Virtual Screening (HTVS), followed by Glide docking SP (standard precision) and XP (extra precision) at the crystal structure binding sites with Glide. The co-crystallized ligand has been centralized for grid generation using the grid generation tools in Glide. MM-GBSA post-docking minimization (Molecular Mechanical Energies Combined with Generalized Born and Surface) was performed to optimize the geometries of the molecules recovered, and the 10% of the molecules recovered from each step were selected for the next level. Finally, all non-peptide molecules (peptide compounds are orally degradable) were run on the Glide XP molecular docking system, using the 3VHE crystal structure to estimate the docking scores of the resulting molecules after screening.

2.8. Induced Fit Docking

A molecular docking method, known as induced fit docking (IFD) [27], where the receptor is flexible in the docking study, has been performed. The energy minimisation of the protein structure was performed using the OPLS_2005 force field. The prepared molecules were docked to the rigid protein using Glide with the default parameters. Energy minimisation was applied on the crystal structure of PDB code: 3VHE. XP molecular docking was used for initial docking and 21 ligand poses were maintained for protein structure refinement. Schrödinger’s 2017-4 Prime module, was used to refine residues within 5.0 Å of the ligand poses and induced fit protein-ligand complexes were developed. After these refinements [28], the ranking of each of the 21 complexes was performed by Prime energy. Complex structures with an energy less than 30 kcal·mol⁻¹ were re-docked for the final step of scoring. Each ligand was docked into each refined low energy receptor structure developed in the refinement step. The binding affinity of each complex was estimated by the docking score. The lowest negative docking score and IFD score was considered as a more favourable binding condition with the active site of 3VHE.

2.9. ADME Prediction

QikProp tool of Schrodinger [29] was applied to predict the druggable property of ten best hits by assessing the ADME profile. During this, the Lipinski rule of five and various descriptors like QPlogHERG, QPPCaco, QPlogBB, and % human oral absorption were calculated.

3. Results and Discussion

3.1. Generation of Pharmacophoric Models

For the modeling of pharmacophores, a set of 74 molecules in the activity range on the logarithmic scale (3.665 - 4.959) were selected. Compounds with a pIC50 activity > 4.4814 were selected as active while those with pIC50 < 4.0409 were considered inactive. The selected active and inactive molecules were used to test the specificity of the pharmacophore hypothesis and to define the excluded volumes. Models of pharmacophores containing 5 sites were generated using three features: hydrogen bond acceptor (A), hydrophobic (H) and aromatic ring (R). The survival score parameter gives a general ranking to all the pharmacophoric models generated and was used as selection criteria for the pharmacophores models. The best pharmacophores models obtained based on the site score, vector score and volume score parameters which lead to the survival score and BEDROC score parameters have been reported in Table 2. A good pharmacophores model is characterized by a high value of the parameters. Survival score. The survival score value of all pharmacophoric models generated varies between 4.911 and 5.129. Model AAHHR-3 had the highest value while model AAAHR-8 gave the lowest value for survival score.

3.2. Analysis of Pharmacophoric Models

The pharmacophoric models constructed were validated using a set of parameters (Table 3) in a database containing 324 decoy molecules, collected using 9 inhibitors of prostate cancer cells (PC3) which were not used in the construction

Table 2. Pharmacophores model and scoring parameters.

Table 3. Characteristics of the best pharmacophores models obtained after validation.

of pharmacophoric models. The pharmacophore obtained were used for the screening of the database using the phase software screening tool implemented in the Schrödinger suite version 17-4. The ability of a pharmacophore to differentiate between active agents and decoys is given by the enrichment factor EF. The EF scores (1%) of the models obtained are in the range 20.81 - 33.30, indicating that these models are all able to identify assets from a large dataset of compounds [30].

The position of the actives in relation to the compounds classified in an ordered manner corresponds to the ROC value. This value varied between 0 and 1 where a value greater than or equal to 0.7 is considered an appropriate performance measurement value [31]. In the present study, ROC value of between 0.98 and 1 indicates that the pharmacophores models obtained have a strong capacity for selecting active molecules. In addition, the % screen graph and the ROC graph (Figure 1) revealed the sensitivity and specificity of recognizing active molecules. BEDROC measurements measure early recognition of database assets and range from 0 to 1 [32]. We considered α = 20.0 for the BEDROC metric, which means that 80% of BEDROC results come from the first 8% of the molecules classified [31]. Therefore, a substantial BEDROC value (α = 20.0) between 0.80 and 0.99 suggested the early detection of active compounds in the database. The selection of the best hypotheses was carried out using the Güner-Henry (GH) scoring method. The analysis of GH notation was carried out by calculating the following variables: %A is the percentage of known active compounds extracted from the database (precision), H_a is the number of assets in the list of results (true positives), A is the number of active compounds in the database, %Y, the percentage of known assets in the results list, H_t is the number of results retrieved, D is the number of compounds in the database, EF is the enrichment of the active concentration by the model compared to the random screening without any pharmacophore approach and GH is the Güner-Henry score. Models 1 and 5 succeeded in recovering 90% of the active compounds and a GH score of 0.92 indicated the good quality of these models (Table 3). All these enrichment results suggest that the pharmacophores models generated are satisfactory for the recovery of the assets of a large database of molecules.

Based on the GH parameter, the models with the greatest power of selectivity for active molecules are AAHHR_5 and AAAHR_8 (GH = 0.922). The AAHHR_5 hypothesis is characterized by two attractor sites (A), two hydrophobic sites (H) and one aromatic site (R). The other model, AAAHR_8 is represented by three attractor sites (A), a hydrophobic site (H) and an aromatic site (R). The distances and angles between the different sites of the models are given in Table 4 and Table 5, respectively and illustrated in Figure 2. Interactions less than 3.1 Å are considered to be strong. Those between 3.1 Å and 3.55 Å are assumed to be average while those greater than 3.55 Å are low [33].

Figure 1. Roc curve of the pharmacophores models AAHHR_5 and AAAHR_8.

Figure 2. Distances and angles between pharmacophoric sites of models AAHHR_5 and AAAHR_8.

Table 4. Distance between the different pharmacophoric sites of models AAHHR_5 and AAAHR_8.

Table 5. Angles between the different pharmacophoric sites of the AAHHR_5 and AAAHR_8 models.

3.3. Development and Validation of 3D-QSAR Models

Using the phase software of Schrödinger suite, we generated a set of 3D-QSAR models using the “atom-based QSAR” module which gradually incorporates PLS (Partial Least Square) regression factors. To build the model the grid spacing parameter was set to 1 Å, the PLS factors were set to 5 and the variables with |t-values| ≺ 2.00 have been eliminated. A 3D-QSAR model was generated with 74 ligands which were randomly divided into the test set (56 ligands) and validation set (18 ligands). The Phase module generated a total of 5 models with the different statistical parameters which are presented in Table 6. The model with a PLS factor equal to 2 was considered for further analysis because a small difference between biological activity and biological activity was obtained. Predicted activity with a standard deviation of 0.10005 (Table 6). The curve validating the accuracy of the 3D-QSAR model is presented in Figure 3 which illustrates the evolution of the predicted biological activity of the molecules as a function of their experimental biological activities. We observe a good linear correlation between the activities predicted by the model and the experimental activities.

Table 6. Experimental and predicted values of the molecules of the test set (training set) and the validation set (test set) based on the PLS factor 2.

Figure 3. Line scatter plot illustrating the correlation of actual activity versus predicted activity for the test (1) and validation (2) set using an atom-based 3D-QSAR model.

The residual scale which is the difference between the experimental activity and the predicted activity was used to rank the predictions. Residues less than 0.8 were considered good predictions while residues between 0.8 and 1.6 were considered weak predictions. Residues greater than 1.6 were considered bad predictions [34]. The results obtained show that all the molecules have good predictions as a function of their residual scale (Table 6).

To provide reliable predictability, a good model must undergo internal and external validation. The predictive power of the generated 3D-QSAR model was analysed using a set of 18 compounds and the statistical significance of the model was obtained using a PLS factor 2. The robustness of the model to predict active molecules was considered as a function of various internal and PLS parameters, including cross-validation coefficient (Q²), the correlation coefficient between predicted and experimental biological activity (R²), the standard deviation (SD), the mean squared error (RMSE), the variance (F), the significance level of variance ratio (P), Pearson’s correlation coefficient (Pearson-r) and model stability (Table 7). For the PLS factor 2, the regression coefficient (R² = 0.9032) on the one hand and cross-validation coefficient (Q² = 0.9077) on the other hand, indicate that the model has good internal predictive power. In addition, the very large value of the variance (F = 247.2) with a small value of P = 1.348e−27, low SD (0.0813), low RMSE (0.08) and high value of Pearson-r 0.9330 confirmed the importance of the selected model (Table 7). There is little difference between predicted and biological activities of the test and validation sets. This small difference highlights the effectiveness of 3D-QSAR model. The value R² − Q² = 0.0336 < 0.3, all internal statistical parameters are well within the defined range,

Table 7. Result of the quantitative structure-activity relationship (QSAR) of the developed model.

SD. regression standard deviation; R² is for the regression coefficient; F is the ratio of the model variance to the observed activity variance (variance ratio); P: significance level of the variance ratio; RMSE: mean square error. Q² directly analogous to R² but based on the predictions of the test set is a validation metric calculated using observations (y-axis) and predictions (x-axis) and the Pearson-r value for the correlation between predicted and observed activity for the test set; RMSE RMS error in test set predictions.

but external statistical validation is also essential for greater reliability of the selected model.

The atomic contributions of the molecules were analysed in order to understand 3D-QSAR model obtained and to clarify the design requirements for developing more potent human liver cancer inhibitors. The result is shown in Table 8. According to the results of the constructed QSAR model, the hydrophobic contributions/nonpolar substituents. Electron-withdrawing groups favorably contribute to the biological activity of the compounds studied.

3.4. Analysis of External Statistical Validation

In order to assess the external predictive capacity of 3D-QSAR model obtained. The Golbraikh-Tropsha parameters and Roy’s metric were used to calculate the external predictive correlation coefficients $r_{\text{pred}}^2$, CCC, $Q_{F_1}^2$ and $Q_{F_2}^2$ [35]. In Table 9, the values of $r_{\text{pred}}^2$ of 3D-QSAR models and high values of CCC, $Q_{F_1}^2$ and $Q_{F_2}^2$ shows the robustness and efficiency of the external prediction capability of 3D-QSAR model. This external prediction capability indicates the reliability of 3D-QSAR model in predicting the biological activity of new compounds. As can be seen in Table 8, the established models also meet the criteria of Golbraikh-Tropsha and Roy. According to the data in Table 10, these criteria are present in the established 3D-QSAR model. Also the result showed that the established 3D-QSAR model was free from systematic error and can therefore be applied to the prediction of the biological activity of the ligands of the test set.

Table 9 shows the performance parameters based on the MAE criteria [36] for the external validation tests of the 3D-QSAR model. If a QSAR model follows the criteria: MAE ≤ 0.1 × range of the test set and MAE ± 3σ ≤ 0.2 × range of the test set then the model can be considered a good predictor. According to the data in Table 9, these criteria are present in the established 3D-QSAR model. Also the result showed that the established 3D-QSAR model was free from systematic error and can therefore be applied to the prediction of the biological activity of ligands of the validation set.

Table 8. Contribution of the atom-type fraction.

Table 9. Statistical properties calculated for the external validation of the developed QSAR models.

Table 10. Results of external statistical validation of the 3D-QSAR model by the MAE method.

3.5. 3D-QSAR Visualization MapAnalysis

Color maps provide information on positions which require a particular physicochemical property to enhance the cytotoxic activity of a ligand. The QSAR model displays the 3D features as cubes. The blue cubes show positive coefficients which are favorable while the red cubes show negative coefficients which are unfavorable characteristics for the activity. These maps give a clue of which functional group is desirable or undesirable at certain positions in a ligand. Their positive contribution is represented by blue cubes and their negative contribution is represented by red cubes. The compounds F64, F38 and F19 were selected as model molecules for the study of the resulting 3D QSAR model.

v Hydrogen Bond Donor interaction

For more active compound F64, the red region is around the oxygen atoms and also around the carbon between two nitrogen atoms of benzimidazole. Cyclohexane and the sulfur atoms of rhodanine are the areas of compound F38 around which the red cubes are concentrated. This indicates that substitutions at these positions by groups with more hydrogen bond donor property do not promote the inhibitory activity of the ligand in the active site of the 3VHE receptor. The blues regions are seen around the nitrogen atoms of the F64 ligand. For compound F38, the blue area is concentrated around the oxygen atoms and also around the carbon-carbon double bond close to rhodanine. Concerning the F19 molecule, the blue region is located around the two carbon atoms (C27, C30) located between the two rhodanine rings. These areas indicate that substitution at these positions with groups with greater hydrogen bond donor properties promotes the inhibitory activity of the ligand in the active site of the 3VHE receptor (Figure 4).

v Hydrophobic interactions

For more active compound F64, the red region is around the carbonyl group. For compound F38, the red cubes are located around the C2 and C4 carbon atoms of cyclohexane. Concerning compound F19, the red zone is located around the C30 carbon between the two rhodanine rings and also around the C=S double bond of rhodanine. Substitution by groups having more hydrophobic properties at these positions therefore it does not promote the inhibitory activity of the ligand in the active site of the 3VHE receptor. The blue regions of the most active compound F64 are observed around the cyclohexane and the piperazine ring including all the atoms located between these two rings. For compound F38, the blue region is observed around the sulfur atoms of rhodanine and also around the oxygen atoms bound to cyclohexane. For the F19 molecule, the blue cubes are located around two C=O double bonds and also the two nitrogen atoms of the rhodanine sites. These areas indicate that substitution at these positions by groups with strong hydrophobic properties promotes the inhibitory activity of the ligand in the active site of the 3VHE receptor (Figure 5).

v Electro-attractor interaction

Visual analysis of compound F64 indicates the presence of blue cubes at the

Figure 4. Overview of contour maps for 3D-QSAR hydrogen bond donor models for F64, F38 and F19.

Figure 5. Overview of contour maps for 3D-QSAR models of hydrophobic interactions for F64, F38 and F19.

benzimidazole and two nitrogen atoms of the piperazine ring. For the F38 molecule, these blue cubes are seen around the oxygen and sulfur atoms. Regarding F19, the blue zone is observed around the carbon-carbon double bond linked to rhodanine. These blue regions indicate that substitution by compounds with more electron-withdrawing properties at these positions promotes the inhibitory activity of the ligand against the 3VHE receptor. The areas represented by the red cubes indicate that substitution with compounds with greater electron-withdrawing properties does not promote the inhibitory activity of the ligand in the active site of 3VHE receptor. For the F64 compounds this zone is located around the carbon atoms of the piperazine ring and the oxygen atoms. For compound F38, it is localized around the two cyclohexane rings and around rhodanine. Regarding F19, we observe the red cubes around all the oxygen atoms of the molecule but also around the C=O double bonds of rhodanines (Figure 6).

3.6. Virtual Screening of the Enamine Chemical Library

The increasing numbers of genomic targets of therapeutic interest [37] and macromolecules (proteins, nucleic acids) for which a three-dimensional (3D) structure is available [36] make virtual screening techniques increasingly attractive for bioactive molecule identification projects [38]. Virtual screening is any computer search process in molecular databases that allows the selection of molecules. In this work, docking studies were performed using the Schrödinger Suite Glide grid-based method [39]. Virtual stepwise screening was performed using the Glide HTVS, SP and XP docking methodologies described above. All docking poses were scored with the MM-GBSA approach, as implemented in the Prime program of the Schrödinger software suite. After the Prime MM-GBSA analysis we retained the molecules with energies below −50 kcal/mol. To take into account the flexibility of the protein, the resulting set of molecules was subjected to the IFD protocol. The sequential virtual screening including the HTVS, SP, XP prime MM-GBSA and IFD protocols allowed us to select a total of 20 new molecules. Table 11 and Table 12 present some parameters of these obtained compounds as well as their 3D structures.

All the molecules obtained have a fitness score greater than 1.5, which means that they are well aligned with the pharmacophore hypotheses that enabled them to be found. The higher the fitness score, the better the structural alignment. The fitness scores of the leads are between 1.507 and 1.775. The compound PV-002558797812 found by the AAAHR_6 pharmacophore model with a fitness score of 1.772 is the highest ranked. The compound Z1694049401 found by the AAAHR_5 model with a fitness score of 1.716 follows. The compound Z1684534882 found by the pharmacophore model AAAHR_7 has the lowest fitness score with 1.507.

Figure 6. Overview of contour maps for 3D-QSAR models of electron-withdrawing interactions for F64, F38 and F19.

Table 11. Molecules retained after virtual screening by pharmacophore models and by molecular docking.

Table 12. 2D structures of the 20 hits obtained.

3.7. HTVS, SP and XP Analysis of Hits

The compounds were subjected to a three-stage docking strategy based on Glide in which all compounds were docked through three stages of the docking protocol, high-throughput virtual screening (HTVS), standard precision (SP) and extra precision (XP). In the first step, the Glide high-throughput virtual screening mode was used and 10% of the best performing ligands were used for the next step, Glide SP. Again, 10% of the best performing leads from Glide SP were retained and docked with Glide XP to refine the correct ligands. The glide energy, glide emodel and docking score parameters of HTVS, SP and XP of the selected hits are given in Table 13.

The compounds obtained by XP screening are classified according to the calculated docking score. These values range from −12.031 kcal/mol to −10.959 kcal/mol. These values are all higher than those of the reference molecule. This means that these 20 hits have an affinity comparable to that of the reference molecule, which confers a better stability in the active site of the 3HVE protein.

3.8. Prime MM-GBSA Hits Analysis

The Prime/MM-GBSA method [40] based on the complex obtained after Docking XP was used to calculate the free enthalpy of binding ΔG_bind of ligands in the active site of the target protein and the results obtained are summarised in Table 14. The free enthalpy of binding ΔG_bind of the hits ranged from −94.520 to −48.197 kcal/mol. All the obtained hits except the compounds Z1683225090 and PV-002569833477 have a binding energy higher than that of the reference ligand.

Van der Waals interactions ΔG_vdW, (ranging from −65.383 to −38.016 kcal/mol), electrostatic or coulomb interactions ΔG_coulomb (ranging from -30.244 to −3.18 kcal/mol) and lipophilic interactions ΔG_lipo (ranging from −34.648 to −19.472 kcal/mol) are the main energetic factors favorable to ligand binding.

Table 13. Glide energy, glide emodel and docking score of the 20 poses after HTVS, SP and XP screening.

Table 14. Calculated free binding energies (kcal/mol) of the 20 selected molecules.

All reported energies are in kilocalories per mole (kcal/mol). ΔG_bind free energy of protein−ligand binding; ΔG_Coul the columbic binding free energy; ΔG_Cov the covalent binding free energy; ΔG_vdW van der Waals binding free energy; ΔG_SolvSA solvation binding free energy of the surface area; ΔG_SolvGB generalized Born solvation binding free energy.

The contribution of hydrogen bonding energy ΔG_H-bond (ranging from −2.37 to −0.566 kcal/mol) and packing energy ΔG_packing (ranging from −10.973 to −0.305 kcal/mol) is small in the free enthalpy of binding.

The unfavourable energy contributions to ligand binding are the covalent interaction energy ΔG_covalent (ranging from 2.566 to 18.012 kcal/mol) and solvation energy ΔG_solvGB (ranging from 7.002 to 42.109 kcal/mol).

These results confirm that the molecules obtained have a higher affinity and therefore a possibly higher inhibition rate than the available reference molecule.

3.9. IFD Hits Analysis

The parameters glide energy, glide Emodel, docking score and IFD score as well as the different types of interactions between the hits obtained and the residues of the active site of the 3VHE protein are summarised in Table 15. The docking score, IFD score, glide Emodel, glide energy parameters of the hits obtained are between −14.07 and −10.813 kcal/mol; −645.847 and −617.2 kcal/mol; −123.445 and −71.983 kcal; −73.135 and −50.605 kcal/mol respectively. The binding modes of the new molecules with the lowest IFDscore (Z1694049401) as well as the one with the highest IFDscore (Z1683053210, PV-001434258746) were discussed in detail.

● Binding Mode of Compound Z1683053210

Of all the compounds selected, the compound Z1683053210 was the best inhibitor in this screening with an IFDscore of −632.286 kcal/mol. Visual analysis shows that this compound forms 2 hydrogen bonds with the active site, one of which is with the Asp 1046 residue and the nitrogen of the 1,2,4-oxadiazole ring with a distance of 2.20 Å and an angle ∠NHN = 164.7˚. The second bond is established between the hydrogen of residue CYS 919 and the oxygen of the morpholine ring with a distance of 1.70 Å and an angle ∠NHO = 164˚. The LYS 868

Table 15. IFD analysis of the 20 selected molecules.

residue engages in a pi-cation bond with the morpholine ring. The benzene group and benzofuran of the compound establish π-π bonds with the Phe 1047 and Hie 1026 residue respectively. Finally, this compound establishes a pi-cation bond with the residue Lys 868 (Figure 7).

Numerous hydrophobic interactions also ensure the stability of the 3VHEPV-001434258746 complex. These interactions are found between the inhibitor and residues Leu 840, Val 848, Ala 866, Leu 889, Val 916, Phe 918, Cys 919, Leu 1019, Leu 1035, Ile 1044, Cys 1045 and Phe 1047 (Figure 8).

● Binding mode of the compound PV-001434258746

The compound PV-001434258746 is considered the second best inhibitor that came out of our virtual screening with an IFDscore of −633.864 kcal/mol. It establishes three hydrogen bonds with the 3VHE receptor. A first bond between the hydrogen of residue Asp1046 and the oxygen of the carbonyl group with a distance of 1.98 Å and an angle ∠NHO = 149.7˚. A second bond is established between the hydrogen of residue Cys919 and the 1,7-naphthyridine nitrogen with a distance of 2.12 Å and an angle ∠NHN = 147.8˚. The last bond is observed between the hydrogen of the nitrogen near the carbonyl group and the oxygen of residue VAL 916 with a distance of 1.78 Å and an angle ∠NHO = 141.5˚. This compound establishes two π-π bonds between its benzene ring and the Hie 1026 residue and then between the 1,7-naphthyridine and the Phe 1047 residue (Figure 9).

Numerous hydrophobic-type interactions also ensure the stability of the 3VHEPV-001434258746 complex. These interactions are observed between the inhibitor and the residues Leu 840, Val 848, Ala 866, Leu 889, Val 916, Phe 918, Cys 919, Leu 1019, Leu 1035, Ile 1044, Cys 1045 and Phe 1047 (Figure 10).

Figure 7. Binding mode of compound Z1683053210 in the active site cavity of the 3VHE receptor.

Figure 8. Hydrophobic interactions (green color) of the Z1683053210 compound with the active site residue of 3VHE.

Figure 9. Binding mode of compound PV-001434258746 in the active site cavity of the 3VHE receptor.

● Binding mode of compound Z1694049401

Z1694049401 establishes six hydrogen bonds with the 3VHE receptor. Through its two oxygen atoms in its two carbonyl groups, it establishes two hydrogen bonds with residue ASP 1046 with distances of 1.97 Å and 2.52 Å making angles ∠NHO equal to 128.8˚ and 158.5˚. The hydrogen of residue CYS 1045 establishes a third hydrogen bond with the oxygen of one of the carboxylic functions with a distance d (H...N) = 2.64 Å and an angle ∠SHO = 123.3˚. Cys residue 919 also establishes two hydrogen bonds, one of which is with the nitrogen of the 1,2,4-thiadiazole and the other with the hydrogen of the nitrogen close to this group with distances of 2.09 Å and 2 Å, respectively. These bonds are made with

Figure 10. Hydrophobic interactions (green color) of PV-001434258746 with residues of the 3VHE active site.

angles ∠NHN = 151.6˚ and ∠NHO = 151.8˚, respectively. The last bond involves the oxygen of the Glu 885 residue and the hydrogen of the nitrogen located between the two carboxylic functions with a distance of 1.69 Å and an angle ∠NHN = 162.5˚ (Figure 11).

It is also necessary to underline the intervention of residues Val 848, Ala 866, Leu 889, Ile 892, Val 898, Val 899, Val 916, Leu 1035, Ile 1044, Cys 1045 and Phe 1047 in the stability of the 3VHE-Z1694049401 complex by the formation of numerous hydrophobic interactions between these residues and the compound Z1694049401 (Figure 12).

3.10. Prediction of Hit Activity

The 3D-QSAR model constructed above was used to predict the inhibitory activity of the twenty molecules (20) selected after virtual screening. The predicted IC50 activity (µM) of the twenty molecules (20) varies between 0.408 µM and 0.817 µM (Table 16). These theoretical activities calculated for the new molecules obtained after virtual screening are significantly lower than that obtained by Coulibaly et al. [15] (between 0.587 µM and 0.756 µM). These new compounds show better or similar biological activities to those of the molecules synthesized by Coulibaly et al. [15]. These molecules could constitute a solid basis in the search for new inhibitors of 3VHE prostate cancer cells.

3.11. ADMET Properties of New Compounds

The ADMET properties of the 20 newly identified hits were evaluated using Qikprop. The results obtained are reported in Table 17. The above 20 hits meet the drug properties based on Lipinski’s rule five. The pharmacokinetic parameters (ADME-Tox) for all 20 hits are within an acceptable range for human use, revealing their potential drug-like properties.

QPlogPo/w is the expected octanol/water partition coefficient; QPlogS predicts the solubility in water. log S. S in mol/l is the concentration of the solute in

Figure 11. Binding mode of compound Z1694049401 in the active site cavity of the 3VHE receptor.

Figure 12. Hydrophobic interactions of Z1694049401 with residues of the 3VHE active site.

Table 16. Predicted activities of new molecules.

Table 17. ADMET properties of the 20 new compounds.

a saturated solution which is in equilibrium with the crystalline solid; QPPCaco predicted the apparent permeability of Caco-2 cells in nm/sec; Caco-2 cells are a model for the intestinal blood barrier; QPPMDCK Expected apparent permeability of MDCK cells in nm/sec; QPlogBB predicted brain/blood partition coefficient.

[Recommended values]: PSA = 7 to 200; QPPCaCo < 25 Poor, >500 Good; QPPMDCK < 25 Poor; >500 Excellent; QlogBB between −3.0 and −1.2; QPlogPo/w between −2 and 6.5; QPlogS between −6.5 and 0.5; QPlogKhsa between −1.5 and 1.5; MW between 130 and 725.

4. Conclusion

In this study, a pharmacophore modeling and virtual screening protocol was performed to identify novel prostate cancer inhibitors. Using in vitro data of prostate cancer inhibitors, we constructed different pharmacophore hypotheses using PHASE software. In addition, the discriminatory ability of the pharmacophore models was validated by enrichment studies using the DECOY method. The developed models show significant value of validation parameters AUC, ROC, BEDROC and RIE. The obtained pharmacophore models were used to screen a set of molecules from the Enamine database. We also developed an atom-based 3D-QSAR model using 56 molecules in the test set and 18 molecules in the validation set. The generated model exhibited a high coefficient of determination (R² = 0.9032) and cross-validation coefficient (Q² = 0.8695) with low RMSE (0.08) and SD (0.0813). The high predictive ability of this model was validated by an external validation test on all molecules in the validation set. Analysis of the atom-based 3D-QSAR contour representation revealed the structural requirements to enhance the biological activity of the inhibitors against the prostate cancer enzyme. Compounds selected by the VSW protocol were subjected to free enthalpy of binding calculation by the Prime-MMGBSA protocol, IFD: induced fit docking protocol, and pharmacological properties calculation (ADMET). The results suggest that all the 20 hits obtained satisfy good in silico criteria such as Docking score, ADME properties, presence of hydrogen bonding interaction and also favorable and lower binding energy and IFD score than the reference molecule.

Conflicts of Interest

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

References

[1]	Parkin, D., Bray, F. and Devesa, S. (2001) Cancer Burden in the Year 2000. The Global Picture. European Journal of Cancer, 37, 466. https://doi.org/10.1016/S0959-8049(01)00267-2
[2]	Jeong, T.-S., Kim, J.-R., Kim, K.S., Cho, K.-H., Bae, K.-H. and Lee, W.S. (2004) Inhibitory Effects of Multi-Substituted Benzylidenethiazolidine-2, 4-Diones on LDL Oxidation. Bioorganic & Medicinal Chemistry, 12, 4017-4023. https://doi.org/10.1016/j.bmc.2004.06.001
[3]	Koppireddi, S., Komsani, J.R., Avula, S., Pombala, S., Vasamsetti, S., Kotamraju, S. and Yadla, R. (2013) Novel 2-(2,4-dioxo-1,3-thiazolidin-5-yl) Acetamides as Antioxidant and/or Anti-Inflammatory Compounds. European Journal of Medicinal chemistry, 66, 305-313. https://doi.org/10.1016/j.ejmech.2013.06.005
[4]	Barros, C.D., Amato, A.A., de Oliveira, T.B., Iannini, K.B.R., da Silva, A.L., da Silva, T.G., Leite, E.S., Hernandes, M.Z., de Lima, M.d.C.A. and Galdino, S.L. (2010) Synthesis and Anti-Inflammatory Activity of New Arylidene-Thiazolidine-2,4-Diones as PPARγ Ligands. Bioorganic & Medicinal Chemistry, 18, 3805-3811. https://doi.org/10.1016/j.bmc.2010.04.045
[5]	Heerding, D.A., Christmann, L.T., Clark, T.J., Holmes, D.J., Rittenhouse, S.F., Takata, D.T. and Venslavsky, J.W. (2003) New Benzylidenethiazolidinediones as Antibacterial Agents. Bioorganic & Medicinal Chemistry Letters, 13, 3771-3773. https://doi.org/10.1016/j.bmcl.2003.07.010
[6]	Trotsko, N., Kosikowska, U., Paneth, A., Wujec, M. and Malm, A. (2018) Synthesis and Antibacterial Activity of New (2,4-Dioxothiazolidin-5-yl/ylidene) Acetic Acid Derivatives with Thiazolidine-2,4-Dione, Rhodanine and 2-Thiohydantoin Moieties. Saudi Pharmaceutical Journal, 26, 568-577. https://doi.org/10.1016/j.jsps.2018.01.016
[7]	Marc, G., Ionuț, I., Pirnau, A., Vlase, L., Vodnar, D.C., Duma, M., Tiperciuc, B. and Oniga, O. (2017) Microwave Assisted Synthesis of 3,5-Disubstituted Thiazolidine-2,4-Diones with Antifungal Activity. Design, Synthesis, Virtual and in Vitro Antifungal Screening. Farmacia, 65, 414-422.
[8]	Tuncbilek, M. and Altanlar, N. (2006) Synthesis of New 3-(Substituted Phenacyl)-5-[3’-(4H-4-oxo-1-benzopyran-2-yl)-benzylidene]-2,4-Thiazolidinediones and Their Antimicrobial Activity. Archiv der Pharmazie, 339, 213-216. https://doi.org/10.1002/ardp.200500180
[9]	Liu, X.-F., Zheng, C.-J., Sun, L.-P., Liu, X.-K. and Piao, H.-R. (2011) Synthesis of New Chalcone Derivatives Bearing 2, 4-Thiazolidinedione and Benzoic Acid Moieties as Potential Anti-Bacterial Agents. European Journal of Medicinal Chemistry, 46, 3469-3473. https://doi.org/10.1016/j.ejmech.2011.05.012
[10]	Ibrahim, M.A., Abdel-Hamed, M.A.-M. and El-Gohary, N.M. (2011) A New Approach for the Synthesis of Bioactive Heteroaryl Thiazolidine-2, 4-Diones. Journal of the Brazilian Chemical Society, 22, 1130-1139. https://doi.org/10.1590/S0103-50532011000600019
[11]	Mousavi, S.M., Zarei, M., Hashemi, S.A., Babapoor, A. and Amani, A.M. (2019) A Conceptual Review of Rhodanine. Artificial Cells, Nanomedicine, and Biotechnology, 47, 1132-1148. https://doi.org/10.1080/21691401.2019.1573824
[12]	Almeida, V.L.d., Lopes, J.C.D., Oliveira, S.R., Donnici, C.L. and Montanari, C.A. (2010) Estudos de relações estrutura-atividade quantitativas (QSAR) de bis-benzamidinas com atividade antifúngica. Química Nova, 33, 1482-1489. https://doi.org/10.1590/S0100-40422010000700011
[13]	Andrade, C.H., Pasqualoto, K.F.M., Ferreira, E.I. and Hopfinger, A.J. (2010) 4DQSAR. Molecules, 15, 3281-3294. https://doi.org/10.3390/molecules15053281
[14]	Cantello, B.C.C., Cawthorne, M.A., Haigh, D., Hindley, R.M., Smith, S.A. and Thurlby, P.L. (1994) The Synthesis of BRL 49653—A Novel and Potent Antihyperglycaemic Agent. Bioorganic & Medicinal Chemistry Letters, 4, 1181-1184. https://doi.org/10.1016/S0960-894X(01)80325-5
[15]	Coulibaly, W.K., Paquin, L., Bénié, A., Bekro, Y.-A., Durieux, E., Meijer, L., Le Guével, R., Corlu, A. and Bazureau, J.-P. (2012) Synthesis of New N,N’-Bis (5-Arylidene-4-Oxo-4,5-Dihydrothiazolin-2-Yl) Piperazine Derivatives under Microwave Irradiation and Preliminary Biological Evaluation. Scientia Pharmaceutica, 80, 825-836. https://doi.org/10.3797/scipharm.1206-04
[16]	Schrödinger Release 2017-4: LigPrep (2017). Schrödinger, LLC, New York.
[17]	Shivakumar, D., Williams, J., Wu, Y., Damm, W., Shelley, J. and Sherman, W. (2010) Prediction of Absolute Solvation Free Energies Using Molecular Dynamics Free Energy Perturbation and the OPLS Force Field. Journal of Chemical Theory and Computation, 6, 1509-1519. https://doi.org/10.1021/ct900587b
[18]	Dixon, S.L., Smondyrev, A.M., Knoll, E.H., Rao, S.N., Shaw, D.E. and Friesner, R.A. (2006) PHASE: A New Engine for Pharmacophore Perception, 3D QSAR Model Development, and 3D Database Screening: 1. Methodology and Preliminary Results. Journal of Computer-Aided Molecular Design, 20, 647-671. https://doi.org/10.1007/s10822-006-9087-6
[19]	Dixon, S.L., Smondyrev, A.M. and Rao, S.N. (2006) PHASE: A Novel Approach to Pharmacophore Modeling and 3D Database Searching. Chemical Biology & Drug Design, 67, 370-372. https://doi.org/10.1111/j.1747-0285.2006.00384.x
[20]	Huang, N., Shoichet, B.K. and Irwin, J.J. (2006) Benchmarking Sets for Molecular Docking. Journal of Medicinal Chemistry, 49, 6789-6801. https://doi.org/10.1021/jm0608356
[21]	Mysinger, M.M., Carchia, M., Irwin, J.J. and Shoichet, B.K. (2012) Directory of Useful Decoys, Enhanced (DUD-E). Journal of Medicinal Chemistry, 55, 6582-6594. https://doi.org/10.1021/jm300687e
[22]	Schrödinger Release 2017-4: Phase (2017). Schrödinger, LLC, New York.
[23]	Braga, R.C. and Andrade, C.H. (2013) Assessing the Performance of 3D Pharmacophore Models in Virtual Screening. Current Topics in Medicinal Chemistry, 13, 1127-1138. https://doi.org/10.2174/1568026611313090010
[24]	Yang, J.M. and Shen, T.W. (2005) A Pharmacophore Based Evolutionary Approach for Screening Selective. Proteins: Structure, Function, and Bioinformatics, 59, 205-220. https://doi.org/10.1002/prot.20387
[25]	Teli, M.K. and Rajanikant, G.K. (2012) Pharmacophore Generation and Atom-Based 3D-QSAR of Novel Quinoline-3-Carbonitrile Derivatives as Tpl2 Kinase Inhibitors. Journal of Enzyme Inhibition and Medicinal Chemistry, 27, 558-570. https://doi.org/10.3109/14756366.2011.603128
[26]	Prakash Tanwar, O., Karthikeyan, C., Hari Narayana Moorthy, N.S. and Trivedi, P. (2010) 3D QSAR of Aminophenyl Benzamide Derivatives as Histone Deacetylase Inhibitors. Medicinal Chemistry, 6, 277-285. https://doi.org/10.2174/157340610793358846
[27]	Wang, H., Aslanian, R. and Madison, V.S. (2008) Induced-Fit Docking of Mometasone Furoate and Further Evidence for Glucocorticoid Receptor 17α Pocket Flexibility. Journal of Molecular Graphics and Modelling, 27, 512-521. https://doi.org/10.1016/j.jmgm.2008.09.002
[28]	Saubern, S., Guha, R. and Baell, J.B. (2011) KNIME Workflow to Assess PAINS Filters in SMARTS Format. Comparison of RDKit and Indigo Cheminformatics Libraries. Molecular Informatics, 30, 847-850. https://doi.org/10.1002/minf.201100076
[29]	Schrödinger Release 2017-4: QikProp (2017). Schrödinger, LLC, New York.
[30]	Samal, S.K., Routray, S., Veeramachaneni, G.K., Dash, R. and Botlagunta, M. (2015) Ketorolac Salt Is a Newly Discovered DDX3 Inhibitor to Treat Oral Cancer. Scientific Reports, 5, Article No. 9982. https://doi.org/10.1038/srep09982
[31]	Truchon, J.-F. and Bayly, C.I. (2007) Evaluating Virtual Screening Methods. Journal of Chemical Information and Modeling, 47, 488-508. https://doi.org/10.1021/ci600426e
[32]	Kumar, A., Roy, S., Tripathi, S. and Sharma, A. (2016) Molecular Docking Based Virtual Screening of Natural Compounds as Potential BACE1 Inhibitors. Journal of Biomolecular Structure and Dynamics, 34, 239-249. https://doi.org/10.1080/07391102.2015.1022603
[33]	Imberty, A., Hardman, K.D., Carver, J.P. and Perez, S. (1991) Molecular Modelling of Protein-Carbohydrate Interactions. Docking of Monosaccharides in the Binding Site of Concanavalin A. Glycobiology, 1, 631-642. https://doi.org/10.1093/glycob/1.6.631
[34]	Sindhu, T. and Srinivasan, P. (2014) Pharmacophore Modeling, 3D-QSAR and Molecular Docking Studies of Benzimidazole Derivatives as Potential FXR Agonists. Journal of Receptor and Signal Transduction Research, 34, 241-253. https://doi.org/10.3109/10799893.2014.885048
[35]	Veerasamy, R., Rajak, H., Jain, A., Sivadasan, S., Varghese, C.P. and Agrawal, R.K. (2011) Validation of QSAR Models-Strategies and Importance. International Journal of Drug Design and Discovery, 2, 511-519.
[36]	Berman, E. and Machin, S. (2000) Skill-Biased Technology Transfer around the World. Oxford Review of Economic Policy, 16, 12-22. https://doi.org/10.1093/oxrep/16.3.12
[37]	Rask-Andersen, M., Masuram, S. and Schiöth, H.B. (2014) The Druggable Genome. Annual Review of Pharmacology and Toxicology, 54, 9-26. https://doi.org/10.1146/annurev-pharmtox-011613-135943
[38]	Lengauer, T., Lemmen, C., Rarey, M. and Zimmermann, M. (2004) Novel Technologies for Virtual Screening. Drug discovery today, 9, 27-34. https://doi.org/10.1016/S1359-6446(04)02939-3
[39]	Schrödinger Release 2017-4: Glide (2017). Schrödinger, LLC, New York.
[40]	Li, H. (2011) A statistical Framework for SNP Calling, Mutation Discovery, Association Mapping and Population Genetical Parameter Estimation from Sequencing Data. Bioinformatics, 27, 2987-2993. https://doi.org/10.1093/bioinformatics/btr509