Abstract

Rhodanine and its derivatives exhibit interesting biological activities as well as a wide range of biological applications. In this study, a dataset of seventy-four molecules with anticancer activities against human cancer cell line Huh-7D12, were chosen for the modeling of pharmacophores and Quantitative Structure Activity (3D-QSAR) relationship. Pharmacophoric models containing five sites were generated from three characteristics: hydrogen bond acceptor (A), hydrophobic (H) and aromatic ring (R). After the validation, eight hypotheses which presented a good power of selectivity of the active agents were selected (GH > 0.5). Internal and external validation parameters indicated that the generated 3D-QSAR model exhibits good predictive capabilities and significant statistical reliability (R² = 0.9606, Q² = 0.955, = 0.952). Pharmacophoric models and contour maps provided significant information on the main structural features of rhodanine derivatives. Twenty-one molecules were returned from the Enamine chemical database after molecular docking studies (HTVS, SP, XP, and IFD). These provided an estimate of ligand-protein binding interactions essential for anticancer activity. The ADMET prediction of these 21 compounds suggested that their pharmacophoric properties lie within an acceptable range. This result indicates that these new compounds provide an effective basis for the methodical development of potent inhibitors of the protein kinase C-kit.

1. Introduction

Cancer (or malignant tumor) is a disease characterized by an abnormally large cell proliferation (tumor) within a normal tissue of the body, in such a way that the survival of the latter is threatened ¹. The main form of treatment is chemotherapy, where anti-tumor chemicals are given to patients. This treatment is considered effective, especially in the early stages of the disease, but it does not always cure the patient or completely destroy cancer. Various factors are linked to the failure of the treatment, among which we can mention the stage of the disease, the resistance of tumor cells to the drugs and the side effects of the treatment as the drugs used kill both the cancer cells and the normal cells, often becoming resistant to treatment ². It is therefore important to develop effective anticancer therapeutic agents with well-defined pharmacokinetic properties.

Rhodanine derivatives have proven to be attractive compounds because they exhibit interesting biological activities as well as important industrial applications. These activities justify the investigation of rhodanine derivatives for the development of new therapeutic agents. The presence of rhodanine backbone in a very wide range of compounds with very varied biological properties makes it an important compound in the search for new drugs. Computer-aided drug design is an ingenious process of finding new drugs based on knowledge of ligands and biological target. It is a complex, lengthy and very expensive process with a high chance of failure at any stage of drug's development. In contrast, computational strategies such as high throughput in silico chemical library screening help prioritize and identify small molecules for organic synthesis and subsequent bioassays while saving time and money.

The objective of this study is to provide new molecules able of effectively inhibiting human liver cancer cells using molecular modeling methods such as pharmacophore modeling, 3D-QSAR and molecular docking. The new molecules obtained were subjected to an Abosrption, Distribution, Metabolism, Elimination and Toxicity (ADMET) in order to study their physiological, physicochemical and pharmacokinetic properties.

2. Materials and Methods

A set of 74 rhodanine derivatives with their anticancer activity (IC₅₀), was synthesized by Coulibaly et al ³ for the development of new inhibitors of human liver cancer. The negative logarithm of the measured IC₅₀ value (pIC₅₀) was used in this study. For 3D-QSAR studies, these 74 compounds were divided into a training set (54 compounds) and a test set (18 compounds). The molecules in the training set were selected at random so that they contained information both on their structural characteristics and on their ranges of biological activity. The molecules in the dataset have been divided into active, inactive, and moderately active molecules.

Compounds with pIC₅₀ activity > 4.886 were considered active while those with pIC₅₀ activity < 4.356 were considered inactive. Compounds with a pIC₅₀ between 4.886 and 4.356 were considered moderately active. All chemical structures and biological activities of all molecules are listed in Table 1.

All calculations were performed on an Intel Core i5-4570 processor with microprocessors running at 3.20 GHz and 6 GB RAM memory running Windows 8.1 operating system. Access to all Schrödinger modules as well as data organization and analysis were performed exploiting Schrödinger software and the Maestro graphical user interface (GUI) ⁴.

3D ligand structures were generated using the construction panel in Maestro and optimized using the LigPrep module ⁵ Partial atomic charges were assigned and possible ionization states were generated at a pH of 7.0 ± 2.0. The OPLS_2005 force field was used to optimize production of low energy ligand conformer Energy minimization was performed for each ligand until it reached a root mean square deviation threshold of 0.01 Å.

For target proteins whose crystal structures were available with their crystal PDB ligands with good resolution were used for the studies. Maestro's Protein Preparation Wizard ⁶ was used to prepare the protein. If the protein structure is a multimer with double binding sites, the redundant site has been removed by choosing molecules or chains. If the binding interaction required the two identical units to form an active site pocket, neither was removed. Water molecules have been removed. A brief relaxation was performed using constrained minimization of all atoms made with Impact Refinement (Impref) module ⁷ using the OPLS-2005 force field to mitigate steric conflicts that may exist in original PDB structures. Minimization ended when energy converged or the RMS value reached a maximum cutoff of 0.30 Å.

Common pharmacophore hypotheses are a spatial arrangement of chemical characteristics common to two or more active ligands that are proposed to explain the key interactions involved in binding of a ligand with its receptor. The PHASE module ⁸ was used for the generation of pharmacophoric models for anticancer agents ⁹. The prepared ligands were imported to develop a Phase pharmacophore model panel with their respective biological activity values. The default parameters were used with a maximum of 10 conformations per rotary connection and with that a maximum of 1000 conformers were generated by the ConfGen module ¹⁰.

The PHASE module provides a standard set of six pharmacophoric characteristics, namely, a hydrogen bond acceptor (A), a hydrogen bond donor (D), a hydrophobic group (H) and a negatively ionizable (N), positively ionizable (P) and aromatic ring (R) to define chemical characteristics of ligands. A total of 6 variants were generated keeping 7 and 4 as the maximum and minimum number of sites respectively. The resulting hypotheses were scored and ranked based on their vector, volume, site scores, survival scores, and survival assets. ⁹. In hypotheses generated, five sites were found to be common for all the selected compounds. The pIC50 activities ranged from 3.866 to 5.699. Compounds with pIC50 activity > 4.886 were selected as active while those with pIC50 activity < 4.356 were considered inactive. Compounds with a pIC50 of between 4.886 and 4.356 were designated as moderately active.

Many metrics are currently used to assess the performance of virtual screening ranking methods. For example, enrichment factor (EF), area under receiver operating characteristic curve (ROC), Boltzmann enhanced discrimination of receiver operating characteristics (BEDROC), area under the accumulation curve (AUAC), the robust initial enhancement (RIE). These metrics were used to determine the robustness of the hypothesis. Evaluating the performance of virtual screening methods are useful for selecting the method that works best in a given context each time, recovering active compounds from a mixed set of active compounds and decoys (presumably inactive compounds against the examined target). The area under Receiver Operating Characteristic (ROC) curve is used by influential groups to measure virtual screening performance in part because it does not exhibit desirable statistical behavior. A value of 0.5 shows that ranking method does not do better than Randomicking. It can be interpreted as the probability that an active will be ranked before an inactive. It has a value between 0 (worst possible performance) and 1 (best performance) ¹¹. Enrichment Factor (EF) is defined as ratio of odds of finding an active compound in the top X% of the data set ¹¹. A maximum enrichment is therefore 100 if all of the actives (A=1) are found within the top 1% of the decoys (D=0.01).

The quality of pharmacophore models was assessed using Güner-Henry (GH) scoring method ¹². The results of GH score assessment are strongly influenced by the composition of active and inactive compounds in database. A GH value greater than 0.5 was used for desirable pharmacophore models. The following criteria were applied to find several pharmacophore models:

Where, Ha: Number of actives in the hit list (true positives); Ht: Number of the hits retrieved; A: Number of active molecules in the database; D: Total number of database compounds; % A: The ratio of the actives retrieved in the hit list (precision); % Y: The yield of actives (recall); EF: Enrichment factor (i.e. enrichment of the concentration of the actives by the model relative to random screening without any pharmacophoric approach); GH: The Guner-Henry score.

The PHASE module presents two options for the alignment of the 3D structure of molecules; pharmacophore-based alignment and atom-based alignment ¹³. In this study, we used an atom-based QSAR model, which is more useful in explaining relationship between activity and structure. In atom-based QSAR, a molecule is treated as a set of overlapping van der Waals spheres. Each atom (and therefore each sphere) is placed in one of six categories according to a simple set of rules: hydrogens attached to polar atoms are classified as donors of H (D) bonds; C - H carbons, halogens and hydrogens 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 atoms are classified as electron attractors (W); and all other types of atoms are classified as miscellaneous (X).

For purpose of QSAR development, the van der Waals models of the molecules in the aligned learning set were placed in a regular grid of cubes, with each cube assigned zero or more ‘’ bits ‘’ to account for different types atoms in training set that occupy the cube. This representation results in binary-valued occupancy models that can be used as independent variables to create QSAR partial least squares (PLS) models. Atom-based QSAR models were generated using the 54 (75%) compound training set with 1.0 Å gate spacing. The best QSAR model was validated by the predicting activities of the 18 (25%) compounds in test set. A five component model (PLS factor) with good statistics was obtained for dataset while maximum number of PLS factors in each model can be 1/5 of total number of training set molecules. QSAR models containing few PLS factors were generated. The statistical quality of the generated QSAR models was judged by parameters such as regression coefficient (R²), cross validation variance (F), confidence interval (P), mean square error (RMSE) and Pearson correlation coefficient ¹⁴.

In order to better assess predictive power and the quality of model developed, the external validation of the model was carried out with the following validation criteria: the metric introduced by Roy and al. ¹⁵ et parameters ¹⁶ and CCC ¹⁷, these parameters are calculated using the equations below:

Y_{obs (test)} and Y_{pred (test)} respectively represent the observed and predicted activities of the validation set and are the correlation coefficient values between the values of the observed and predicted biological activities of the compounds of validation set with a straight line that passes through the origin or not Changing the axes gives the value of is the constant of the correlation line (at the origin) (predicted values according to the experimental values)is the constant of the (original) correlation line (experimental values versus predicted values), is the correlation coefficient of concordance. For a model with good external predictability, the values of and must be greater than 0.5 and must be less than 0.2, As for the values of and CCC, they must be greater than 0.5, 0.5 and 0.85 respectively.

Tropsha ¹⁸ proposed the following criteria for model validation:

Besides the metric and the Tropsha criteria, the criteria based on MAE analysis (Mean Absolute Error) ^{15, 19} were used to thoroughly test external predictability of the QSAR model. In this approach, a QSAR model has good predictive power if it meets these two conditions:

• MAETest ≤ 0.1 × training set range and MAETest + 3 × σ ≤ 0.2 × training set range].

Where, MAETset is error on the prediction of biological activity of 95% of the compounds in validation set (excluding 5% of compounds whose error on prediction of biological activity is high) and σ is standard deviation of absolute value of errors for the validation set. Criteria based on MAE analysis and other validation measures were estimated using XternalValidationPlus 1.2 software developed by Roy et al, ¹⁵. This program also checks if a systematic error in prediction has occurred.

An Enamine database ²⁰ which includes more than 1323686 commercial compounds, was used as a database for the screening. The database was processed by the phasedb_manage and phasedb_confsites functionalities of PHASE, in order to generate a 3D PHASE database in which multiple conformers, as well as their corresponding pharmacophoric sites, were created and stored for each molecule. A so-called thorough sampling method was used, thus generating sets of conformations for the basic structures of the compounds, while all other parameters were left by default. The PHASE 3D database was screened using the pharmacophoric models generated previously. Then, it was imposed that only compounds corresponding to at least four of the total characteristics of the pharmacophoric models are recovered.

XP ligand docking was performed rather than SP docking because XP is better than SP in scoring function and it also predicts false positive results. This docking was performed in Glide of Schrödinger-Maestro ²¹. The end result of the docking can be found as an energy minimization slip score. For docking, the van der Waals scale factor was set at 0.85 and 0.15 for the ligand compounds and the partial charge cut-off value was set at -10.0 kcal / mole. The lowest slip score containing compounds was then subjected to MM-GBSA analysis for the calculation of free binding energy and the best poses were recorded for each ligand compound.

Calculation of binding free energy was also performed for protein ligand complexes. MM-GBSA is a combined method for the calculation of the binding free energy that was used in this experiment which accumulates energies of OPLSAA molecular mechanics (EMM), an SGB solvation model for polar solvation (GSGB) and a non-polar solvation term (GNP) composed of the surface accessible to non-polar solvents and the van der Waals interactions ²². The best poses of Glide score were used for calculation of free binding energy. The total free energy of the bond:

ΔE_MM: is energy difference between the protein-ligand complex and the sum of energies of the protein (without ligand) and of the ligand. ΔG_solv: is the difference in GBSA solvation energy of the complex and the sum of solvation energies of the protein and the ligand. ΔG_SA: is the difference in surface energy of the complex and the sum of surface energies of the protein and the ligand.

The best leads which have a free enthalpy of ΔGbind binding greater than 50 kcal / mol after MM-GBSA analysis were coupled to 1t46 receptor ²³ in order to study the affinity of protein-ligand complex using Induced Fit Docking (IFD) protocol) ²¹ implemented in Schrödinger. The docking modulus taking into account the IFD flexibility of the protein is one of the most robust and accurate methods which takes into account the flexibility of ligands and receptors.

In this method, the ligands were first docked to the protein with a van der Waals radius scale of 0.5 for the protein and ligand atoms. The protein molecule has been subjected to energy minimization with the OPLS-2005 force field ²⁴ and an implicit solvation model. Initial docking was performed using SP mode and the number of poses generated was set at 20 for protein refinements. The Prime module ²⁵ was used to develop the refined protein-ligand complexes with each of the 20 structures from the previous step. All residues having at least one atom located within 3.5Å of each pose of the corresponding ligand were included in the Prime refinement. This leads to a structure and a conformation of the ligands which are adapted to each pose of the structure of the receptor. The complexes thus obtained are classified with Prime energy parameter and complexes having an energy located at less than 30 kcal / mol of the minimum energy structure were used for Docking and Scoring XP course.

In the last step, each ligand was re-docked into each low energy receptor structure produced in the second step using Glide XP with default settings. An IFD score that takes into account both protein-ligand interaction energy and total system energy was calculated and used to rank IFD poses.

For the analysis of physiological descriptors of a compound such as adsorption, distribution, metabolism and excretion behavior of ligand compounds, ADME analysis was performed in Schrodinger's QikProp module ²⁶. It also predicts the physico-chemical nature of compounds as well as their pharmacokinetic properties. In this study, we used Schrodinger's Qikprop 3.2 module. There are also several other descriptors also analyzed such as Predicted IC50 to block HERG K + channel in vitro, predicted octanol or water partition coefficient [log P (o / w)], number of hydrogen bond acceptors (HBA), number of hydrogen bond donors (HBD), estimated aqueous solubility (log s), skin permeability (log Kp), MDCK cell permeability (MDCK), binding to human serum albumin (log Khsa), blood brain partition coefficient (logBB), percentage of human oral absorption rate.

3. Results and Discussion

The main aim of this study was to identify novel inhibitors of human hepatoma cancer. Schrödinger's PHASE module was used to identify pharmacophoric models, while 3D-QSAR based on the atoms of aligned ligands helped to decipher substitution effects. ADME computation and molecular docking were also performed to enhance the scope of the study and investigate ligand-receptor interactions.

For the modeling of pharmacophores, a set of 74 rhodanine derivatives in the activity range on logarithmic scale (3.866-5.699) was selected. This set was subdivided into two subsets, a first consisting of 56 compounds which made it possible to establish the 3D-QSAR models and a second consisting of 18 compounds which made it possible to validate the models obtained. Compounds with pIC50 activity > 4.886 were selected as active while those with pIC50 activity < 4.356 were considered inactive. Compounds with a pIC50 of between 4.886 and 4.356 were designated as moderately active.

The selected active and inactive molecules were used to test the specificity of the pharmacophore hypothesis and to define the excluded volumes. Pharmacophoric models containing five sites were generated using three characteristics: hydrogen bond acceptor (A), hydrophobic (H) and aromatic ring (R). The Survival score parameter, giving a general classification to all pharmacophoric models generated, was used as selection criteria for pharmacophoric models. The best pharmacophoric models obtained based on the site score, vector score and volume score parameters which lead to Survival score and BEDROC score parameters were reported in Table 2. A good pharmacophoric model is characterized by a high value of Survival score parameters. The Survival score parameter of all pharmacophoric models generated varied between 5.123 and 4.631.

Pharmacophore validation is essential to verify the effectiveness of the model in asset recovery and early classification ¹⁵. For this purpose, different enrichment parameters were calculated and presented in Table 3.

The enrichment factor (EF) is the quality that a pharmacophore has in distinguishing active ingredients from decoys. It is calculated based on the number of known assets identified from a fraction of the database examined (Dror et al., 2009). Emphasis was placed on the EF (1%) because it measures the enrichment score for the top 1% of the lures detected. The EF scores (1%) of the models obtained are between 19.59 and 33.30, indicating that these models are all able to identify assets from a large dataset of compounds ³⁴. The ROC value corresponds to the position of the assets with respect to the compounds classified in an ordered manner which are arranged linearly among the defined internal library ²⁷. The ROC value varying between 0 and 1 where a value greater than or equal to 0.7 is considered as the appropriate execution measurement value ²⁸. Therefore, the ROC value between 0.98 and 0.99 in this study, reflects the good capacity of pharmacophore models for the selection of 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 ¹⁵. We considered α = 20.0 for the BEDROC metric, which means that 80% of the BEDROC results come from the first 8% of the molecules classified ²⁸. Therefore, a substantial BEDROC value (α = 20.0) between 0.788 and 0.99 suggested the early detection of active compounds from the database. The selection of compounds with the greatest power of selectivity was carried out using the Güner-Henry (GH) scoring method. The GH values of the selected models between 0.631 and 0.922, therefore greater than 0.5 indicate good reliability of the models, which suggests that they could be used successfully in virtual screening studies to find new chemical entities against human liver cancer. Models 3 and 8 succeeded in recovering 90% of the active compounds with a GH score of 0.922. All these enrichment results suggest that the pharmacophoric models generated are very satisfactory.

Based on the GH parameter, the models with the greatest power of selectivity for active molecules are AAAHR_3 and AAHHR_4 (GH = 0.922). The AAAHR_3 hypothesis is characterized by three attractor sites (A), a hydrophobic site (H) and an aromatic site (R). The other model, AAHHR_4 is represented by two attractor sites (A), two hydrophobic site (H) and one aromatic site (R).

A good 3D-QSAR model should exhibit reliable predictability that can be validated by validation methods, i.e. internal and external validation. The predictive power of the generated 3D-QSAR model was analyzed using a set of 18 compounds and the statistical significance of the model was obtained using a PLS factor of 3. The robustness of the model to predict active molecules was considered as a function of different internal parameters and PLS, including regression coefficient for training set (Q²), regression coefficient for test set (R²), standard deviation (SD), square error mean (RMSE), variance (F), significance level of variance ratio (P), Pearson correlation coefficient (Pearson-r) and model stability (Table 4). For the PLS factor 3, the regression coefficient for the test set between the predicted biological activity and the experimental biological activity (R² = 0.9606) on the one hand and the regression coefficient of the set of training (Q² = 0.9556) on the other hand, indicate that the model has good internal predictive capacity. In addition, the large value of the variance (F = 422.3) with a smaller value of P = 2.618e-034, the low SD (0.0125), the low RMSE (0.07) and the high value of Pearson-r 0.9777 confirmed the importance of the selected model (Table 4). The difference R²-Q² equal to 0.005 is much less than 0.3. From the results of the internal validation it can be concluded that this model is stable and has very good explanatory power, which is very acceptable ²⁹.

The Phase module generated total of 5 models with the different statistical parameters which are presented in Table 4. The model with a PLS factor equal to 3 was considered for further analysis because a small difference was obtained between biological activity and predicted activity with a standard deviation of 0.005. (Table 4). The curve validating the accuracy of 3D-QSAR model is shown in Figure 2 and Figure 3, in which the predicted biological activity of molecules as a function of biological activity is shown. The scatter plots of observed activity versus predicted activity for generated 3D-QSAR model indicate correlation coefficients of R² = 0.91 for the training set (Figure 2) and R² = 0, 96 for validation set (Figure 3). This result testifies to the good predictive power of the considered 3D-QSAR model.

The predictions were ranked according to the residual scale which is the difference between the experimental activity and the predicted activity. 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 ³⁰. In the present study, all the molecules showed good predictions according to their residual scale (Table 5).

In order to understand 3D-QSAR model obtained, the atomic contributions of the molecules were also analyzed to elucidate the design requirements to develop more potent human liver cancer inhibitors and the result is tabulated in Table 6. According to the results of the constructed QSAR model, the hydrophobic contributions / nonpolar substituents and electron withdrawers favorably contribute to the activity. The presence of electron-withdrawing groups also plays an important role.

We cannot judge the predictive capacity of the model developed from internal validation for any new set of compounds. The external validation was therefore carried out to assess the predictive capacity of the model developed for the prediction of untested molecules. Table 7 through Table 9 show external validation metrics obtained using xternal validation software.

The following statistical parameters obtained by the validation set verify the conditions and acceptability conditions of Golbraikh and Tropsha, thus establishing the predictive capacity of the model obtained:

In order to verify the proximity between the observed and predicted data, the parameter (test), was developed by Roy et al. ³¹. The value of (test) uses the square correlation coefficients between the observed property and the prediction of tested compounds. For an acceptable forecast, the value of (test) should preferably be less than 0.2 provided that the test value or greater than 0.5. The correlation coefficient of concordance (CCC) proposed by Lin ³² which measures both the precision and the trueness of the model was also evaluated. With a CCC value = 0.9771 < 1, the model obtained is very satisfactory.

Table 9 shows performance parameters based on MAE criteria (Roy et al) (citation) for external validation testing of 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 as 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 for the prediction of the biological activity of the ligands of the test set.

The color maps around the ligands obtained by PHASE indicate that the ligand-based 3D QSAR analysis can define the characteristics and the important effects of various substitutions on the rhodanine derivatives for their anticancer action. The graph represents the positive and negative activity coefficients of different properties, namely (a) hydrogen bond donor, (b) electron withdrawing property and (c) hydrophobic / nonpolar properties. Their individual positive contribution is represented by blue cubes and their negative contribution is represented by red cubes. This study concerns three ligands: the most active (F45), the least active (F18) and the moderately active ligand (F49).

• Hydrogen Bond Donor Interaction

For the most active compound (F45), the red region is observed around carbon atoms bonded to nitrogen in the rhodanine ring and also around C24 carbon bonded to nitrogen in the piperazine ring. That of the moderately active compound (F49) is observed around the piperazine ring and also around the three carbon atoms located between the nitrogen of the rhodanine ring and that of the piperazine ring. This red zone is observed essentially around the cyclohexane of the F18 ligand. This indicates that substitutions at these positions by groups with more hydrogen bond donor property do not promote inhibitory activity of the 1t46 receptor. The blue regions of the F45 ligand are seen around rhodanine, the N52 atom close to the carboxylic function outside of the rhodanine ring and around the oxygen atom of the hydroxyl group of cyclohexane. For compound F45, these blue regions are located around the piperazine ring, the nitrogen close to the carboxylic group and the oxygen atom bonded to rhodanine. Compound F18 has blue areas around the nitrogen atoms located on the two rhodanine rings and around the carboxylic functions of the rhodanine rings. These areas indicate that substitution at these positions by groups with greater hydrogen bond donor properties promotes the inhibitory activity of 1t46 (Figure 4).

For compound F45 (more active), the red region is located around the sulfur atom included in the five-membered ring of rhodanine. Concerning the F49 molecule (moderately active), this zone is visible around the nitrogen close to cyclohexane but also around the carbon atoms and nitrogen located between the piperazine ring and the carboxylic function of the ligand. The F18 molecule (less active), the red cubes are observed around the oxygen atoms of the carbonyl groups of the rhodanines and also around cyclohexane. Substitution by groups having more hydrophobic properties at these positions does not promote the inhibitory activity of the 1t46 receptor.

Blue regions are seen around the nitrogen atoms of the rhodanine ring and the piperazine moiety of compound F45. For the F49 ligand, this region is observed carbon-oxygen (C = O) and carbon-sulfur (C = S) double bonds of rhodanine but also around the nitrogen of the piperazine ring. Those of the F18 molecule are located around the sulfur atoms. These areas indicate that the substitution at these positions by groups with strong hydrophobic properties promotes the inhibitory activity of the 1t46 receptor (Figure 5).

For the F45 molecule (more active), the red zone is located around the nitrogen of the rhodanine ring and also around the butyl group located between the piperazine ring and the rhodanine ring. For the moderately active compound F49, the rhodanine ring, nitrogen close to cyclohexane and the sulfur integrated into the rhodanine ring are the areas around which the red cubes are located. Concerning the F18 ligand, this zone is located around the rhodanine rings and also around the cyclohexane. These areas indicate that substitution with compounds with greater electron-withdrawing properties does not promote inhibitory activity of the 1t46 receptor. Blue areas are visible for compound F45 around nitrogen and oxygen atoms bonded to the carbonyl group including this group itself. For compound F49, blue cubes are observed around the carbon - oxygen (C = O) and carbon - sulfur (C = S) double bonds of rhodanine. The sulfur atoms of the rhodanine of compound F18 are located in a blue area. 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 1t46 receptor (Figure 6).

The growing numbers of genomic targets of therapeutic interest ³³ and macromolecules (proteins, nucleic acids) for which a three-dimensional (3D) structure is available ³⁴ make virtual screening techniques more and more attractive for projects to identify bioactive molecules ^{35, 36}. By virtual screening is meant any electronic search process in molecular databases allowing the selection of molecules. In this thesis, docking studies were performed using the Glide grid-based method ³⁷ of the Schrödinger suite. The virtual step-by-step screening was performed using the Glide HTVS, SP and XP anchoring methodologies described above. All mooring poses were noted with MM-GBSA approach, as implemented in the Prime program of the Schrödinger software suite. After the Prime MM-GBSA analysis we selected molecules with an energy of less than -50 kcal / mol. To take into account the flexibility of the protein, the set of molecules that were obtained were subjected to IFD protocol. Sequential virtual screening including HTVS, SP, XP prime MM-GBSA and IFD protocols allowed us to select a total of 21 hits. Table 10 presents some parameters of the 21 hits obtained, in particular the number of sites on the ligand which corresponds to the pharmacophoric hypothesis (Matched Ligand Sites) which made it possible to find them during the virtual screening.

All the molecules obtained have a fitness score greater than 1.5, which means that they are well aligned with the pharmacophoric hypotheses which made it possible to find them. In fact, the higher the value of the fitness score, the higher the quality of the structural alignment. The fitness scores for the leads are between 1.528 and 1.742. The compound Z1685826991 found by the pharmacophoric model AAAHR_6 with a fitness score of 1.742 is the highest classified. This is followed by the compound Z1683786797 found by the model AAAHR_7 with a fitness score of 1,738. The lead PV-001894667629 found by the pharmacophoric model AAAHR_5 shows the lowest fitness score with 1.528.

Compounds were subjected to a Glide-based three-level docking strategy in which all compounds were anchored by three steps of the docking protocol, High Throughput Virtual Screening (HTVS), Standard Precision (SP) and extra precision (XP). In the first step, Glide's 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 tracks from Glide SP were retained and were docked with Glide XP to refine the correct ligands. The glide energy, glide emodel and docking score parameters of HTVS, SP and XP for the selected leads are shown in Table 11.

The compounds obtained by XP screening are classified according to the docking score to be calculated. These values range from -12.031 kcal / mol to -10.959 kcal / mol. These values are all higher than those of imatinib, which is the reference molecule. This means that these 21 leads have an affinity comparable to that of the reference molecule, which confers better stability in the active site.

The 2D structures of the new leads obtained are shown in the following table:

To estimate the relative binding affinity of ligands and to compromise experimental binding affinities, free energies were calculated by the Prime / Molecular Mechanics / Generalized Born Surface Area (MM / GBSA) approach ³⁸. This method based on the complex obtained after Docking XP was used to calculate the ΔG_bind free enthalpy of ligands in the active site of the protein. The mean values of the free binding energy (ΔG_bind) obtained from the MM-GBSA calculations are presented in Table 13. The free enthalpy of binding ΔG_bind of the leads are between -82.896 and -48.197 kcal / mol. All the hits obtained have a binding energy greater than that of the reference ligand.

Van der Waals ΔG_vdW interactions (between -65.383 and -38.016 kcal / mol), electrostatic interactions or coulomb interactions ΔG_coulomb (between -30.244 and -3.18 kcal / mol) and lipophilic interactions ΔG_lipo (including between -34.648 and -19.472 kcal / mol) constitute the main energy factors favorable to the binding of ligands. With higher values, the contribution of the energy of the hydrogen bond ΔG_H-bond (between -2.37 and -0.566 kcal / mol) and that of the energy of packing ΔG_packing (between -10.973 and - 0.305 kcal / mol) is low in the free enthalpy of binding. The unfavorable energy contributions to ligand binding are the ΔG_covalent covalent interaction energies (between 2.566 and 18.012 kcal / mol) and ΔG_solvGB solvation (between 7.002 and 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.

In order to take into account, the flexibility of the molecules and the protein target, the 21 compounds taken out of the XP docking were subjected to an IFD docking. The ligands shown in Table 14 are based on the IFD score. This table summarizes the values of glide energy, glide emodel, docking score and IFD score and also presents the residues of the active cite which interfere in different types of interactions with each ligand. The glide energy, glide emodel, docking score and IFDscore respectively have values between -73.135 and -50.605 kcal / mol ; -123.445 and -71.983 kcal / mol ; -14.07 and -10.813 kcal / mol ; -645.847 and - 617.2 kcal / mol.

A more negative IFD score indicates better interaction of the inhibitor with the target protein. The IFD docking revealed that all compounds have a lower docking score than that of Imatinib (-7.747 kcal / mol) co-crystallized with the protein 1t46. On the other hand, only 2/3 of these compounds, ie 14 compounds, have an IFD score that is lower than that of Imatinib (- 621.974 kcal / mol). This suggests that these compounds have a higher binding affinity towards 1t46. The remaining seven (7) compounds have IFD scores between - 621.241 and - 617.2 kcal / mol and close to that of Imatinib, which implies that these compounds could have binding mode similar to that of reference molecule.

We found it interesting to see the positioning and to understand the mechanisms involved in the interaction of the best compounds with the target active site.

• Binding mode of the compound PV-000039856911

With an IFD score of -617,200 kcal / mol, the compound PV-000039856911 appears as the best inhibitor out of this screening. Visual analysis shows that this compound inhibits the 1t46 receptor by forming 4 hydrogen bonds with the residues Glu 671, Cys 673, Cys 809, Asp 810. The first hydrogen bond is established between the hydrogen atom of the nitrogen of the benzene group and the oxygen of the carboxylic group of residue Glu 671. The geometric parameters of this bond are d (H…O) = 2.08 Å ; NHO = 160.1°. The second hydrogen bond is established between the oxygen of the carboxylic group of the benzene group of PV-000039856911 and the hydrogen of the nitrogen near this ring according to the geometric parameters d (H…O) = 4.21 ; ∠NHO = 156.6°. The third and fourth hydrogen bond is established between the oxygen atom of the carboxylic group of PV-000039856911 the sulfur hydrogen atom of the residue Cys 809 and the hydrogen atom bonded to the nearest nitrogen of residue Asp 810. The respective geometrical parameters of these connections are d (H…O) = 1.86 ; ∠SHO = 160.6°; d (H…O) = 1.94 ; ∠NHO = 169.7°. The lead PV-000039856911 also establishes a π-π bond with the HIE 790 residue. The PV-000039856911 - 1t46 complex is also stabilized by numerous hydrophobic intersections precisely with the residues Leu 644, Ile 653, Val 654, Leu 783, Leu 799, Ile 808, Cys 809 and Phe 811 (Figure 7).

• Binding mode of compound Z1685831021

The compound Z1685831021 from the Enamine chemical library with the second lowest value of the IFD score (IFDscore = - 642,782 kcal/mol). This score can be explained by the presence of several bonds with the active site of the protein. It establishes three hydrogen bonds with the residues Glu 640 (d (H…O) = 2.79 Å; ∠NHO = 134°), Lys 623 (d (H…N) = 1.90 Å ; ∠NHN = 172.5°) et CYS 673 (d (H…N) = 1.87 Å ; ∠NHN = 166.5°). The established hydrogen bonds are all less than 3.1 Å and the angles greater than 120 °. This means that this compound can effectively inhibit human liver cancer cells. The compound Z1685831021 establishes a π-π bond with the residue Tyr 672 (Figure 9).

Numerous hydrophobic-type interactions also ensure the stability of the Z1685831021-1t46 complex. These interactions are observed between the inhibitor and the residues Leu 595, Val 603, Ala 621, Leu 637, Leu 644, Leu 647, Val 668, Tyr 674, Cys 673, Leu 799, Cys 809 and Phe 811 of the active site of the enzyme (Figure 10).

• Binding mode of compound Z1684496597

The compound Z1684496597 from the Enamine chemical library with the lowest value of the IFDscore (IFDscore = - 645.847 kcal/mol). This is lower than that of imatinib which is the reference compound. Its score results from the establishment of two hydrogen bridges between this compound and the active site of 1t46. The first is formed between the nitrogen atom of the first phenyl group of the inhibitor and the nitrogen atom of the residue Cys 673. The geometric parameters of this connection are d (H…N) = 1.85 Å ; ∠NHN = 176°. The second hydrogen bond is established between the nitrogen of 1,2,4-oxadiazole and the nitrogen atom of residue ASP 810 with a distance d (H…N) = 3.02 Å and an angle NHN = 174°. These strong interactions (< 3.1 Å) suggest that the compound Z1684496597 is capable of effectively inhibiting human liver cancer cells (Figure 11).

Inhibitor Z1684496597 also establishes two other bonds including a π-π bond between phenyl and the Phe 811 residue and a Pi-cation bond between the Lys 623 residue and 1,2,4-oxadiazole. It is also necessary to underline the intervention of residues Val 603, Ala 621, Leu 644, Val 654, Tyr 672, Cys 673, Leu 799, Cys 809 and Phe 811 in the stability of the complex 1t46 - Z1684496597 allowing the formation of numerous interactions hydrophobic type (Figure 12).

The 21 leads retained after molecular screening were aligned and then used with the 3D-QSAR model constructed above to predict their inhibitory activity (Table 15). The predicted activity of the twenty-one hits (21) varies between 0.2286 µM and 0.492 µM, which is higher than that obtained by Coulibaly et al.

Due to the importance of the bioavailability of a probable drug, we assessed absorption, distribution, metabolism, excretion, and toxicity (ADMET) characteristics via the Insilco ADME prediction tool, QikProp. For this purpose, therapeutically relevant descriptors were examined for the molecules tested among which QPlog Po / w, QPlog S, QPPCaco, QPlogBB, QPlogMDCK and PSA were selected to study their pharmacokinetic behavior.

The partition coefficient (QPlogPo / w) and the water solubility (QPlogS), which are critical for estimating the absorption and distribution of drugs in the body, vary between 1.956 and 4.964 and -7.602 to -3.697, respectively. Cellular permeability (QPPCaco), a key factor governing drug metabolism and its access to biological membranes, ranges from 95.338 to 3002.477.195. The predicted apparent permeability of MDCK cells, QPPMDCK, ranges from 38.999 to 8440.378 nm/s. The blood-brain partition coefficient (QPlogBB) ranges from -1.821 to 0.208. The percentage of human oral absorption of the compounds ranged from 73.821 to 100%. The polar surface (PSA) is another physicochemical characteristic which represents the zone accessible to the solvent of the compound which will interact with the solvent by dipole interaction or hydrogen bonding. It was observed that the hits obtained have a PSA between 56.319 and 119.13 Å, that is to say that they are within the acceptable limit (maximum 120 Ӓ). The pharmacokinetic properties of Enamine molecules are within the acceptable range defined for human use, indicating their potential as drug-like molecules.

4. Conclusion

The work presented in this study includes the development of pharmacophore models for receptor inhibition (PDB: 1t46) performed using PHASE. The best pharmacophores models retained after validation by enrichment studies have very satisfactory AUC, BEDROC, RIE and GH validation parameters. These models were subsequently used to screen a set of molecules from the public Enamine library.

An atom-based 3D-QSAR model was constructed from a set of 74 rhodanine derivatives subdivided into two subsets, 56 for the test set and 18 for the validation set. The model obtained presented a good statistical significance and an important predictive capacity (R² = 0.9606, Q² = 0.955, = 0.952). In addition, visualization of the 3D-QSAR models helped to correlate the chemical group / substituent at different molecular sites of interest with biological activities and provided guidance for structural modification for analog design more powerful. In addition, the pharmacophore model was used to carry out, from the Enamine database, a selection of new molecules not yet studied. The ligands from this high-scoring chemical database were then subjected to molecular docking studies (HTVS, SP, XP and IFD) providing an estimate of ligand-protein binding interactions essential for anti-cancer activities. Further study of 21 compounds from this screen, using ADME prediction, suggested that their pharmacokinetic properties are within an acceptable range. The present study should provide an effective guide for the systematic development of potent inhibitors of the protein kynase C-kit.

References