Virtual Screening of PEBP1 Inhibitors By Combining 2D/3D-QSAR Analysis, Hologram QSAR, Homology Modeling, Molecular Docking Analysis and Molecular Dynamic Simulations

doi:10.21203/rs.3.rs-997081/v1

Download PDF

Research Article

Virtual Screening of PEBP1 Inhibitors By Combining 2D/3D-QSAR Analysis, Hologram QSAR, Homology Modeling, Molecular Docking Analysis and Molecular Dynamic Simulations

https://doi.org/10.21203/rs.3.rs-997081/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Human Phosphatidylethanolamine Binding Protein 1 (hPEBP1) is a novel target affecting numerous cellular signaling pathways involving in the formation of metastases. It can be used in the treatment of many cases of cancer. For these reasons, the Pharmaceutical Companies use computational approaches, including Multi-QSAR (2D, 3D, and Hologram QSAR) analysis, Homology modeling and molecular docking analysis and molecular dynamic simulations to accelerate the drug discovery process. In this paper, QSAR modeling was conducted using two quantum chemistry optimization methods (AM1 and DFT levels). As per PLS results, we found that DFT (B3LYP(d, p)) method present high predictability according to 2D-QSAR, CoMFA, CoMSIA, and Hologram QSAR studies, with Q² of 0.81, 0.67, 0.79 and 0.67 respectively, and external power with R²_predof 0.78, 0.58, 0.66 and 0.56, respectively. This results has been validated by CoMFA/CoMSIA graphics, which suggests that electrostatic fields combined with hydrogen bond donor/acceptor fields are beneficial to the antiproliferative activity. While the Hologram QSAR models show the contributions of each fragment in improving the activity. The results from QSAR analyses revealed that Ursolic Acids with heterocyclic ring could improve the activities. However, Ramachandran plot validated the modeled PEBP1 protein. Whereas, molecular docking and MDs simulations revealed that the hydrophobic and hydrogen bond interactions are dominant in the PEBP1's pocket. These results, were used to predict the structures of three new compounds with high anticancer activity In Silico, Similarly, molecular docking stability studies and molecular dynamics simulations were conducted.

Spectroscopy

Materials Chemistry

Ursolic Acid heterocyclic

Multi-QSAR

DFT (B3LYP/6-31G) levels

PEBP1 target

Molecular dynamic simulations (MDs)

Cancer continues its alarming progression worldwide with 18.1 million new cases and 9.6 million mortalities in 2018, according to the GLOBOCAN 2018 (Global cancer statistics) [1]. “Human phosphatidylethanolamine-binding protein 1 (hPEBP1)”, novel member of the PEBP family newly genes recently updated in the Universal Protein Resource UniProt (https://www.uniprot.org/), is implicated in many cellular processes, including “signal transduction, cell cycle, inflammation, adhesion, proliferation, differentiation, apoptosis, autophagy, circadian rhythm and mitotic spindle checkpoint” [2]. Furthermore, this protein encodes a member of PEBP1 family of proteins and regulate multiple signaling pathways, namely the MAPK, NF-kappa B [3], and “glycogen synthase kinase-3 ( GSK-3 signaling pathways) [4]. In addition, it can also inhibit the RAF1 kinase activity through inhibiting its activation and by dissociating the RAF1/MEK complex and acting as a competitive inhibitor of “MEK phosphorylation” [5]. Also, it is associated with many human cancers by acting as a metastasis suppressor gene. Likewise, with multiple query sequences in the genome, it is necessary for the commitment of scientists and organizations to explore novel approaches for discovering the specific role of PEBP1 protein, that would open us a new ways in cancer prevention and therapy.

Ligand-based Computational approaches like 2D-QSAR analysis are widely investigated to describing the relationship between biological activities and molecular descriptors. These molecular descriptors are physicochemical properties, various topological indices, Molecular shape, Structural fragments, and Quantum chemistry (QC) descriptors [6]. The selection of the accurate quantum chemical methods is of central importance in 2D-QSAR studies [7]. There are many quantum chemical methods used in the geometry's optimization of molecules, we can cite the fast semi-empirical method (AM1) in terms of calculation which takes into account only valence electrons with the adjustable parameters [8]. While “the density functional theory (DFT) method” including the electron-correlation to calculate energy of the system directly from the electron density [8]. Therefore, the choice of the quantum chemical method may produce errors in the calculation of the molecular descriptor and, ultimately, a lack of adjustment in 2D-QSAR model development. Whereas, In CoMFA approach, we calculate the steric and electrostatic potentials around the molecules and we relates fields variations calculated to the biological activity [9]. While, the CoMSIA strategy is introduced as extension of CoMFA analysis [10], with the five similarity indices which occupy the non-covalent interactions, namely hydrophobic effects, hydrogen-bond donor/acceptor potentials. However, as a result of the CoMFA and CoMSIA analyses, Contour areas are undertaken to locate the favored or unfavored regions surrounding molecules that would be needed to increase the biological activity. On the other hand, CoMFA and CoMSIA requires the best alignment of molecules to get better results, which poses a major problem in many cases of the 3D-QSAR studies. In order to overcome this problem, the newly Hologram QSAR (HQSAR) approach [11], which does not take into account molecular alignment in the QSAR model development is used. HQSAR analysis is based on the contribution of each fragment in the molecule to biological activity.

Homology modeling is commonly known approach as one of structure based methods that can predict with high accuracy the 3D structure of a protein from its amino acids sequence. Therefore, computational structural determination methods are needed to bridge this growing gap between the number of available sequences and the 3D structures of experimentally resolved proteins [12]. Furthermore, in order to explore the binding affinity of the comparative protein model, molecular docking simulation is widely used to select the best ligand-protein interactions [13]. On the other hand, molecular dynamics simulations are widely used to predict the stability of the protein-ligand interactions in the biological system [14].

Ursolic acid (3β-hydroxy-urs-12-en-28-oic acid) is a pentacyclic terpenoid, usually present in the stem bark, leaves or fruit peel [15]. Recently, it has been reported that ursolic acid is one of the naturally abundant pentacyclic triterpene acid has vast pharmacological activities including antidiabetic, antioxidative [16], antiHIV, antihepatodamage, antimalarial, antimicrobial, cardiovascular, immunomodulatory effects, anti-inflammatory, and antitumor activities [17],[18]. However, UA has the following characteristics: low toxicity, liver protection and potential anti-cancer metastasis [19, 20]. In addition, recent progress in the clinical trial proved that Ursolic acid exhibit numerous anticancer effects with acceptable toxicity [21–23]. However, it is imperative to continue the research of the mechanism of action and signaling pathway studies of ursolic acid as a new anticancer drug.

The aim of the present study is to compare the accuracy of 2D-QSAR study by different theory levels, namely AM1 and DFT/(B3LYP/6-31G) methods by using Partial Least Squares (PLS) analysis. Meanwhile, the 3D-QSAR analysis has been undertaken to select the structural requirement needed to enhance the antiproliferative activities of ursolic acids through CoMFA, CoMSIA and HQSAR analyses with these two quantum optimization methods. Whereas, the homology modeling was carried out on PEBP1 query sequence as new target protein and its binding pocket was investigated by molecular docking simulation to explore the potential affinity of ursolic acid heterocyclic derivatives against the modeled PEBP1 protein.

At the end of this study, we predicted the design of three new heterocyclic usolic acid derivatives with higher in silico activities than the most active compound (M30) of the series studied, and then evaluated their binding affinity and stability respectively through molecular docking and molecular dynamics simulations.

2.1 2D-QSAR modeling

2.1.1 Dataset treatment

A series of new synthesized Ursolic Acid (UA) derivatives (32 compounds of which only 30 with determined antiproliferative activities) has been selected from the literature [24] in order to develop, among others, 2D/3D and Hologram-QSAR models. In this study, the antiproliferative activities IC₅₀ values (in µM) were converted to -Log IC₅₀, as reported in Table 1, then, all compounds were sketched and initially optimized by MM2 force field in ChemBioOffice software (15.1.0, PerkinElmer, Waltham, MA, USA). After that, we selected randomly 25 compounds as training set for QSAR development and 5 compounds as a test set to check the predictability of the retained QSAR models.

2.1.2 Computational details

Ursolic Acids geometry optimization was performed with Gaussian 09 software [25] using two methods. The first method that has been used is semi-empirical, with the optimization that has been carried out by the "AM1 method", then the energy that has been calculated by the "ZINDO method". The second method adopted is empiric with the optimization that was determined by the DFT method [26] with (B3LYP/6-31G (d,p)) basis set [27, 28], and the “Time-Dependent Density Functional Theory (TD-DFT)” method [29, 30], which was used to calculate the absorption energy. On the basis of the calculation methods used, the descriptors determined are: the energy of lowest unoccupied molecular orbital (E_LUMO)”, “the energy of the highest occupied molecular orbital (E_HOMO)”, the frontier orbital energy (Egap), the dipole moment (DM), the energies E(_ZINDO) and E_(TD−DFT), the oscillation factor (f), the excitation energy (ΔE) and the maximum wavelength (λmax) (Suppllementary Data: Table S1). Supplementary topological descriptors were determined by using ChemBioOffice 15.1 program. Consequently, two separate matrices were obtained, each comprising 38 descriptors, and then statistical tools were undertaken to select the most relevant among them. Those data sets were applied to determine the influence of quantum chemical methods on the QSAR models performance.

The various calculations were performed using two computers: (Computer 1: Intel(R) core(TM) i5-4210M CPU @ 2.60GHz 2.60 GHz. Computer 2: Intel(R) core(TM) i5-4210M CPU @ 2.60GHz 2.60 GHz.

2.1.3 Statistical testing

2.1.3.1 Partial Least Squares method (PLS)

Statistical analyses were conducted through QSARINS [31] and XLSTAT software’s [25], then 2D-QSAR studies were investigated using PLS analysis at the 99% confidence level. The models quality were assessed with R-squared (R² coefficient), R² adj ( R² adjustment), the RMSE (Root Mean Square Error), Student-test and its probability (t-test,p-value), the Fisher-test (F value, p-value) [32, 33]. Furthermore, the collinearity between descriptors in the model was fixed by the correlation matrix and “the variance inflation factor (VIF)” [32, 34].

2.1.3.2 2D-QSAR model assessment

Before any employment of 2D-QSAR study in the prediction of new compounds, the relevant statistical metrics must check the predictive power of the retained 2D-QSAR model [35, 36]. The criteria commonly employed to check the internal reliability of the QSAR model is the “leave-one-out (LOO)” strategy [37]. The latter allows the Q² coefficient to be determined as shown in Equation 1. For an acceptable 2D-QSAR model, the Q² must be greater than 0.5 [37]. Furthermore, the “Y-randomization test” [38] was examined to avoid chance correlation between the retained descriptors and antiproliferative activity in the selected QSAR model. However, the external predictive capability was tested by the R²_pred metric (Equation 2).

\({\text{Q}}^{2}=1-\left[\frac{{\sum }_{\text{i}}{\left({\text{y}}_{\text{i}\text{t}\text{r}\text{a}\text{i}\text{n}}-{\text{y}}_{\text{i}\text{p}\text{r}\text{e}\text{d}\left(\text{L}\text{O}\text{O}\right)}\right)}^{2}}{{\sum }_{\text{i}}{\left({\text{y}}_{\text{i}\text{t}\text{r}\text{a}\text{i}\text{n}}-{\stackrel{-}{\text{y}}}_{\left(\text{L}\text{O}\text{O}\right)}\right)}^{2}}\right]\) Eq. (1) \({\text{R}}_{\text{P}\text{r}\text{e}\text{d}}^{2}=1-\left[\frac{{\sum }_{\text{i}}{\left({\text{y}}_{\text{i}\text{t}\text{e}\text{s}\text{t}}-{\text{y}}_{\text{i}\text{p}\text{r}\text{e}\text{d}\left(\text{t}\text{e}\text{s}\text{t}\right)}\right)}^{2}}{{\sum }_{\text{i}}{\left({\text{y}}_{\text{i}\text{t}\text{e}\text{s}\text{t}}-\text{y}\text{t}\text{r}\text{a}\text{i}\text{n}\right)}^{2}}\right]\) Eq. (2)

Here, y_i and y_ipred are respectively the experimental and predicted pIC₅₀ values of each molecule (Training set/Test set), \({\stackrel{-}{y}}_{Train}\) represent the average pIC₅₀ values of the molecules in the training set.

2.1.3.3 Applicability domain (AD)

Before any in silico prediction of the activities of new molecules, we should define the AD as reported by the Organisation for Economic Co-operation and Development (OECD) QSAR validation. The AD of the developed QSAR model is defined as the chemical space in which the biological activity of new molecules could predict with confidence [39]. The prediction is considered reliable only if the predicted compound is situated inside the domain of applicability. The data set was examined for the possible existence of outliers through William plot analysis. Herein, the plot shows the standardized residuals in relation to leverage. In this case, if the high lever (h) is greater than the warning lever (h*), the compounds have a significant influence on the predicted model [40, 41].

Table 1. Molecular structure of Ursolic Acid derivatives with their antiproliferative activities (pIC50)

Cpds ID	R 1	R 2	R 3	R 4	pIC50	Cpds ID	R 1	R 2	R 3	R 4	pIC50
1	OH	H	OH	H	4.89	16	=O	H	Imidazole	H	4.81
2	OAc	H	OH	H	4.14	17	=O	H	1-Methyl- Imidazole	H	4.70
3	OAc	=O	OH	H	4.47	18	=O	H	Triazole	H	4.71
4	=O	H	OH	H	4.61	19	=O	=O	Imidazole	H	4.93
5*	=O	=O	OH	H	4.26	20	=O	=O	1-Methyl- Imidazole	H	5.29
6	=O	H	OH	OH	4.87	21	=O	=O	Triazole	H	5.00
7*	OH	H	OMe	H	4.75	22	=O	H	Imidazole	Imidazole	5.29
8	=O	=O	OMe	OH	5.04	23*	=O	H	1-Methyl-Imidazole	1-Methyl-Imidazole	4.94
9	=O	H	OMe	OH	5.14	24	=O	H	Triazole	Imidazole	5.40
10	OAc	H	Imidazole	H	4.81	25	=O	=O	OMe	Imidazole	5.63
11	OAc	H	1-Methyl-Imidazole	H	4.69	26	=O	=O	OMe	1-Methyl-Imidazole	5.67
12	OAc	H	Triazole	H	4.63	27*	=O	=O	OMe	Triazole	5.09
13	OAc	=O	Imidazole	H	4.85	28	=O	H	OMe	Imidazole	5.23
14	OAc	=O	1-Methyl-Imidazole	H	4.81	29*	=O	H	OMe	1-Methyl-Imidazole	5.67
15	OAc	=O	Imidazole	H	4.99	30	=O	H	OMe	Triazole	5.72
*Test set

2.2. 3D-QSAR Methods

2.2.1. Molecular alignment

Since the quality of a 3D-QSAR model is based on the choice of the correct molecular alignment [42], for the comparison and the development of the optimal 3D-QSAR models, we have investigated alignment by two different optimization methods to reach the best quantum level that will allow to have the best alignment. Alignment I, as depicted in Fig (1.a). was studied by the semi-empirical (AM1) method and DFT(B3LYP) optimization method was investigated to get alignment II (Fig 1.b). From both alignments, the distill-rigid strategy as introduced in SYBYL-X 2.0 software [43] was carried out the select the common scaffold. The highly active Compound 30 (pIC50=5.72) was taken as molecule model for alignment and all remaining 29 molecules were superimposed on this.

2.2.2 CoMFA strategy

The CoMFA approach was used to develop the reliable 3D-QSAR model based on the aligned compounds. “The Steric (S) field was calculated by Lennard−Jones potentials and Coulombic (Electrostatic) potentials were calculated in 3D grid regularly spaced (2.0 Å between two consecutive points)” [40, 44]. A carbon-probe sp3 hybridized with + 1 charge was fixed to calculate the steric (S)/electrostatic (E) potentials effects [40]. To reduce noise in electrostatic and steric calculations, the default cut-off energy and Column filtering were fixed at 30 kcal/ mol and 1.0 kcal/ mol, respectively.

2.2.3 CoMSIA strategy

The CoMSIA strategy was done with the standard grid box as employed in CoMFA studies. Besides (S) and (E) potentials, the CoMSIA hydrophobic (H) and hydrogen-bond donor (HBD) and hydrogen-bond acceptor (HBA) descriptors were calculated to assess the similarity between compounds. The CoMSIA descriptors were calculated using a carbon-probe [40] with a +1 charge for H, HBD and HBA at each lattice points. We have set column filtering at 1.0 kcal/ mol for improving the performance of the CoMSIA model.

2.2.4. Partial least squares (PLS) analysis and 3D-QSAR models assessment

The PLS method was used to correlate the CoMFA and CoMSIA descriptors to biological activities for purposes of developing reliable 3D-QSAR models [40]. The Q² metric (equation 1) has been employed to set the optimum number of components (NOC) being to the high predictive power of the established 3D-QSAR models. As a result, the Q² metric (Q²> 0.5) allowed us to check the internal reliability of 3D-QSAR models [45, 46]. Furthermore, the external assessment of the developed 3D-QSAR models was assessed by the external set (Test set: 5 molecules). Therefore, the predictive R-squared coefficient (R²_pred) was calculated to assess the external reliability of the retained 3D-QSAR model.

2.3 Hologram QSAR (HQSAR)

3D-QSAR methods have demonstrated to be reliable tools in drug design. However, the complexity of the choice of optimization methods (stable geometries) and the structural alignment method makes serious problems in the development of 3D-QSAR model [11]. Furthermore, Hologram QSAR (HQSAR) is a new approach developed to solve the problems facing 3D QSAR analysis. HQSAR uses the molecular holograms (fragment fingerprints) as descriptors to get which fragments contribute to the activity [11]. In this study, HQSAR analysis was carried out by using the following parameters: hologram length (HL) values (53 to 401) and the fragment distinction; Atoms (A), Bonds (B), Connections (C), Hydrogen atoms (H), Chirality (Ch), and Donor/Acceptor (D/A). The setting atom count in fragments was set between 4 and 7 atoms in each fragment. The PLS technique was used to develop the HQSAR model and then was validated by Q², R² and RMSE values.

2.4 Homology modeling

2.4.1 Templates searching and modeled PEBP1 protein quality

The amino acids sequence of PEBP1 Homo sapiens (ID: P30086) has been uploaded from further homology modeling. The BLASTp program was enrolled to search the best templates protein with a high identity with the query sequence. After alignment search, we have chosen the crystalized structure (PDB ID: 6ens) as the template protein. This template was uploaded from PDB database (https://www.rcsb.org/). The homology model of PEBP1 was built with SWISS-MODEL web server [47], then the accuracy of the developed homology model was checked by Ramachandran plot [48] and the GMQE score (Global Model Quality Estimation) and Qmean score [49].

2.5 PEBP1 model -Ursolic Acid Docking

Once the homology model validation is completed, we have used CavityPlus web server to identify the binding sites necessary to get the best binding affinity. The cavity pocket is constructed with the following amino acid residues (Fig. 2): Asp-70, Asp-72, Ala-73, Pro-74, Ser-75, Lys-80, Tyr-81, Arg-82, Trp-84, His-85, His-86, Val-107, Pro-111, Thr-115, Gly-108, Gly-110, Ser-109, Pro-112, His-118, Tyr-120, Ser-142, Gly-143, Asp-144, Arg-146, Gly-147, Lys-148, Leu-180, Tyr-181, Leu-184, Ser-185. Subsequently, the highly active compounds 30 (pIC₅₀= 5.72), 29 (pIC₅₀= 5.67), 26 (pIC₅₀= 5.67) and less active compound 1 (Ursolic acid without heterocyclic ring) were taken as template ligands for docking studies. Based on the results obtained from 2D, 3D and hologram QSAR models, we are selected DFT(6-31G(d, p)) optimization method to elucidate the correct binding affinity of studied ligands. While, Autodock Tools was used to prepare the entry PDBQT files for Autodock vina docking program. The binding site of PEBP1 model has been enclosed by grid box with the following dimensions; the size of the points in each axis was set at (x = 60; y = 60; z = 60) then the Grid center was set at (Center x = 17.0, Center y =-6.928, Center z =20.733). The best poses were carried out based energy score (E in Kcal/mol). Discovery Studio software was used to visualize PEBP1-Ursolic Acid Interactions.

2.6 Molecular Dynamic simulation and system building

For better understanding of the stability of the PEBP1-Ligand complexes, molecular dynamic simulations (MDs) wer carried through the conformations obtained by molecular docking studies. The molecular dynamics (MD) simulations were carried out for the highest active compound M30 (pIC₅₀ = 5.72) and the novel predicted ligands (X1, X2, and X3) in the binding pocket of the modeled PEBP1 receptor. To do this, each PEBP1-ligand complex was solvated in a cubic water box following by the minimization of the solute-solvent system by MMFF94 force field with a maximum of 5000 steps with a gradient threshold of 1 kcal/mol/Å. After that, all simulations were run following a canonical ensemble (NTV) using the AMBER 7 F99 force field [50, 51]. Snapshots were taken at 250 fs intervals during 50 ps of simulation with a temperature of 300 K. All molecular dynamic simulation were performed using SYBYL-X 2.0. As a result, the potential and kinetic energies were plotted in order to evaluate the stability of the studied complexes during the 50 ps of simulations. On the other hand, the root-mean square deviation (RMSD) was investigated though calculation the distance between the initial minimized structure as the reference structure (t=0) and obtained conformation in each trajectory during 50ps of simulations. While, the compactness of the protein folding was elucidated by calculated the Radius of gyration (Rg).

3.1. 2D-QSAR analysis

We calculated a data set of 38 descriptors based on quantum mechanics methods (AM1 and DFT) and molecular mechanic methods (MM2) (for details, see Table A1, ESI†). The partial least squares (PLS) method allowed us to select a set of 14 descriptors from the 38 initial descriptors.

3.1.1 Selection of descriptors

Due to the multicollinearity problem inherent in the development of 2D-QSAR model, correlation analysis is an essential tool for refining the most correlated descriptors with the antiproliferative activity. From the analysis of the results in Fig. 3a, we can deduce that the two principal axes F1 and F2 obtained by DFT-based method represent 71.12% and 15.47%, respectively of the total variance, while those obtained by AM1-based method represent only 46.19 and 27.88 respectively (Fig. 3b). The correlation matrices for the both (AM1, DFT) methods are shown in Electronic Supplementary Material (ESI†: Tables A2 and A3).

Based on those results, we have selected four descriptors determined by AM1-based method and the topological index: Egap, C3, HBA and TV connectivity, whereas, four descriptors were selected by using the DFT and topological index methods: Egap, E_TD-DFT, f and HBA. The matrix dataset for the retained descriptors in both models is tabulated in Supplementary Material (Table A4).

3.1.2 PLS regression

Using the antiproliferative activities values (dependent variable) and the Quantum mechanics descriptors ( independent variables), the established 2D-QSARs models through PLS regression method are shown as the following below:

AM1-based method:

pIC ₅₀ = 7.98 – 0.41 * Egap + 1.108 * C3 + 0.155 * HBA + 1.13 * TVConn(x10⁶) Eq. (3)

N = 25 R ² = 0.77 R² adjusted = 0.73 RMSE = 0.19 F = 17.86 α = 1%

DFT-based method:

pIC ₅₀ = 8.50 – 0.418 * Egap + 0.00005 * E _TD−DFT + 26.25 * f + 0.19 * HBA Eq. (4)

N = 25 R ² = 0.87 R²adjusted = 0.85 RMSE = 0.14 F = 34.41 α = 1%

With N represents the number of molecules (in the training set), R2 is the coefficient of determination, R²adj represents the coefficient of determination adjusted, F is the global significance of the retained 2D-QSAR model and RMSE its Root Mean Square Error. The significance level α was set at 1%.

From the above results, we can observe that the statistical quality of the model in equation (4) is more significant than that obtained for the model in equation (3). The coefficients R² and R²adj in DFT-based model is remarkably superior than obtained in the AM1-Based model with a lower error RMSE too. However, The F value of each model indicates that the selected model in Eq. (4) is more significant than observed in Eq. (3). The values of the Student test and its resulting probability (p-value) for each descriptor included in the retained 2D-QSAR models are listed in Table 2. We can note that the p values of the t-test are all less than 0.01, therefore, the coefficients of both models are acceptable. Moreover, the VIF values in the eq. (3) and (4) indicates that there is no collinearity between descriptors in both 2D-QSAR models (Table 2). The observed versus predicted activities from both retained 2D-QSAR models are shown in Fig. 4. The predicted antiproliferative activities of the compounds with their residues from Eq. (3) and (4) are given in Supplementary Material (Table A5).

Table 2

Significance of all descriptors by using Student test and the VIF values.
Descriptors	Coeff	Standard Error	t-test	p-value	VIF
AM1- Based model
Egap	- 0.41	0.12	- 3.22	0.004	2.81
HBA	0.15	0.04	3.07	0.006	1.91
TVConn	1135377	396908	2.86	0.01	2.13
C3	1.11	0.39	2.78	12	1.21
DFT – Based model
Egap	- 0.41	0.07	- 5.89	< 0.001	2.27
E_TD−DFT	0.000051	0.000015	3.43	0.001	2.45
f	26.25	7.07	3.71	0.001	1.13
HBA	0.19	0.05	3.70	0.003	3.88

3.1.3 Predictive powers of retained 2D-QSAR models

The internal predictive capacity of the two selected 2D-QSAR models was tested by cross validation (LOO) procedure. Q² value (0.55) of the Eq (3) is relatively lower than that observed in Eq (4) (0.81). These results show that AM1 based model has a lower predictive power than the DFT-based model. The correlations between observed vs predcited activities values through internal and external validations are illustrated in Fig. 4.

The randomization test was used to avoid the chance correlation between activities and descriptors [52]. As we can see in Table 3, the initial models show high stability compared to randomized models, which allowed us to conclude that there is a low probability to find a chance correlation in the retained 2D-QSAR models. However, the selected test set was used in this study to prove the external predictive power of the both selected 2D-QSAR models. The model obtained by DFT level has better external predictive power (R²_pred=0.78) than that obtained by the AM1 level (R²_pred=0.68).

Table 3

Values of Q² and R² metrics after multiple Y-Randomization Test
Iteration	Randomization test for both obtained 2D-QSAR models
	DFT based model		AM1 based model
	R²	Q²	R²	Q²
Initial model	0.87	0.81	0.77	0.57
Random 1	0.29	-0.07	0.08	-0.32
Random 2	0.29	-0.09	0.14	-0.59
Random 3	0.19	-0.16	0.16	-0.20
Random 4	0.06	-0.14	0.08	-1.02
Random 5	0.09	-0.29	0.19	-1.89
Random 6	0.06	-0.29	0.04	-1.33
Random 7	0.02	-0.17	0.16	-2.54
Random 8	0.14	-0.69	0.14	-0.49
Random 9	0.34	-0.01	0.12	-0.24
Random 10	0.17	-0.78	0.06	-0.46

The theory level of descriptor calculation is a critical step in 2D-QSAR model development. Thus, from the above presented results, we can see the effect of the chosen level of theory on performing 2D-QSAR models. The best statistical quality was obtained by DFT–based method (Eq. 4) with R-squared (R²)value of 0.87 and the best predictive internal/external power with Q² = 0.81 and R²pred=0.78, respectively. The second 2D-QSAR model is based on the AM1-based method, with the significant coefficients R²= 0.77, Q² is 0.55, and R²pred=0.68. The significant difference in the results obtained for both developed 2D-QSAR models could probably be explained by the fact that the semi-empirical AM1 method does not include the effects of instantaneous electronic correlation [53], compared to DFT method (B3LYP/6-31G). Our results are in agreement with those of the other comparative studies [54],[55], which have proved that the DFT (B3LYP/6-31G) method is the most appropriate method to determine the various properties of molecules.

From both 2D-QSAR models, the energy Egap is negatively correlated with antiproliferative activity, indicating that a decrease in Egap of Ursolic Acid derivatives is favorable to increase the antiproliferative activity. In order to ensure the significance of the energy Egap = E_LUMO – E_HOMO in the DFT based model, an analysis of correlation was carried out. As shown in Fig 5, the E_HOMO (AM1) energies are weakly correlated with E_HOMO (DFT) energies (R² = 0.40), while the relation between E_LUMO (AM1) and E_LUMO (DFT) energies is highly correlated (R² = 0.95).

The main objective of this analysis is to prove that the E_LUMO energies calculated by the two methods (AM1 and DFT), have no effect on the variation of Egap energies, consequently, the observed variations are due to the variation of E_HOMO energies. As we see in Figure 6, for the highest active compound 30 (pIC₅₀ = 5.72) as the template molecule to elucidate the molecular orbitals (LUMO and HOMO), the orbital LUMO of both quantum level is located around heterocylic group (acceptor) at R4 position, while the orbital HOMO (Donor) is shown around the Ursolic acid scaffold near to R3 position. From this analysis, we can deduce that the significance of Egap in the DFT-Based method is more accurate than that calculated in the AM1-Based method. To assess the importance of the heterocyclic group in the calculation of the Egap energy, another comparative study based DFT method was carried out with or without heterocyclic group. Ursolic acid (compound 1) without heterocyclic group has a high Egap value of 6.28 compared to compound 30 (Egap = 4.22), which indicates that the presence of heterocyclic groups like imidazole is favorable for reducing Egap energy, thus an increasing in the antiproliferative activity. In addition, the E_TD−DFT descriptor presents a positive contribution to the antiproliferative activity, which means that the more the energy of the molecule is lower more the molecule is stable. The active conformer 30 could validate this analysis, which presents the lowest E_TD−DFT= - 44620 eV compared to the less active compound M1 (E_TD−DFT=-38035 eV). The observed E_TD−DFT energies of the two compounds could be due to heterocyclic rings grafted to Ursolic acid (compound 1), which contribute favorably to improve its antiproliferative activities.

Other topological descriptor retained in the model (Eq. (3)) is the Mullikan charge (C3) calculated by AM1 method, negatively correlated to pIC₅₀, indicating that the presence of the N-alkyltriazole in R4 conjugated to an unsaturated ketone in R1 affords a positive charge of the compound 30, which is favorable to improve the antiproliferative activity. Whereas TVconn, and Oscillator Strength (f) correlates positively to pIC₅₀. Another descriptor involved in AM1-based model is Hydrogen Bond Acceptor (HBA), which is a topological descriptor derived from the acceptor of nitrogen and oxygen atoms (in this study the heterocyclic group).

3.1.4 Applicability domain (AD)

The Graphs of the domain of applicability of developed 2D-QSAR model based DFT method is depicted in Fig. 7. The analysis of this Figure shows that the compound 1 has a leverage above the threshold h* (h*= 0.6) and therefore these compounds are considered outside of the applicability domain.

3.2 3D-QSAR analysis

3.2.1 CoMFA results

The statistical metrics of obtained CoMFA models for both studied alignment methods are illustrated in Table 4. The CoMFA model generated by Align II (DFT- based optimization method) shows better performance than that obtained in Align I (AM1-based optimization method) with values of Q², R², and R²pred respectively equal to 0.67, 0.92, and 0.58, and a low value of the RMSEtrain of 0.1. In both alignment methods, we find a slight difference between the steric and electrostatic field's contributions. The observed and predicted activities with their residual values are presented in Supplementary Material (Table A6) and the relationship between the observed versus predicted pIC₅₀ values of training / test set is shown in Fig. 8.

3.2.2. CoMSIA results

In the case of CoMSIA study, the statistical results are presented in Table 4. From the five similarity indices of CoMSIA, we have identified that E, HBD, and HBA descriptors are favorable to predict the antiproliferative activity of Ursolic Acid derivatives. As we can see in Table 4, the Align II showed better results than Align I with Q², NOC, R², F, RMSEtrain and R²pred of 0.79, 2, 0.88, 86.46, 0.12 and 0.66, respectively, and a slightly difference observed between E, HBD, and HBA fields contributions. The observed and predicted activities with their residual values are presented in supplementary material (Table A6) and the relationship between the observed versus predicted pIC₅₀ values of training and test sets are shown in Fig. 8.

3.2.3 HQSAR results

According to HQSAR analysis (Table 4), the identification of the fragment that best predicts the antiproliferative activity indicates that the model using best fragment distinction (atoms (A), bonds (B), connections (C) and Donor/Acceptor (DA)) gives significant results with a hologram length equal to 83. The selected hologram shows high-validated coefficients with values of Q², NOC, R², RMSE and R²pred equal to 0.67, 4, 0.87, 0.15 and 0.56 respectively. As per Table 4, we can see that no significant change in the observed results of the two alignment methods. From this analysis, we can conclude that HQSAR model is independent of molecular alignment, which in accord with literature [11, 56]. The plot of the observed versus predicted activities of training/test sets is shown in Fig. 9.

Table 4

Results of CoMFA, CoMSIA and HQSAR studies
Statistical parameters	Align I (AM1-Based)			Align II (DFT-Based)
Statistical parameters	CoMFA	CoMSIA	HQSAR	CoMFA	CoMSIA	HQSAR
Q²	0.62	0.70	0.67	0.67	0.79	0.67
R²	0.88	0.84	0.86	0.92	0.88	0.87
NOC	5	2	5	3	2	5
F	53.84	62.7	83 (HL)	58.7	86.49	83 (HL)
RMSEtrain	0.12	0.15	0.15	0.1	0.12	0.15
R²pred	0.54	0.62	0.56	0.58	0.66	0.56
Probability of R²cv	0	0		0	0	-
Fields contributions
Steric (S)	0.48	-	-	0.53		-
Electrostatic (E)	0.51	0.4	-	0.46	0.45	-
H-bond donor (HBD)	-	0.32	-	-	0.3	-
H-bond Acceptor (HBA)	-	0.27	-	-	0.23	-

3.2.4 Contour maps of CoMFA, CoMSIA, and HQSAR

The CoMFA and CoMSIA contour plots were shaped in order to show the contributions of each selected field in enhancing the antiproliferative activity of Ursolic acid derivatives. The maps obtained represents regions with standard coefficients over 80% (favorable) or under 20% (unfavorable). The active conformation 30 (pIC₅₀ =5.72) is taken as template molecule to elucidate the contour maps.

3.2.4.1 CoMFA contour plots

As per CoMFA contour areas, the green colors show the regions near to R4 position, indicating that this region must be occupied by bulky steric groups (Fig 10.a) which would enhance the antiproliferative activity. Another green region is observed around R3 position. These results could be validated, on the one hand, by the active conformation 29 (pIC₅₀=5.67) that contains the 1-methylimidazole group at R4 and the methyl group at R3, which explains its potential activity. On the other hand, the lack of steric groups at the R3 and R4 positions considerably reduces the antiproliferative activity. Furthermore, the least active compounds 1, 2, 3, 4, and 5 contain no steric group at these positions, reflecting the weaker antiproliferative activities. Therefore, the yellow region, showed at R1 position, indicating area where an increase in steric volume would decrease the antiproliferative activity. In the case of CoMFA electrostatic contour maps, the blue color areas represent the favorable regions and the red colors areas represent the unfavorable regions for the antiproliferative activity. As you see in Fig 10.b, the big contour blue is around the triterpenoid scaffold, indicates that hydrogen bonds would be the reason for high electropositive interactions around the triterpenoid. The red region near to R3 and R1 positions in the active compound 30 denotes that the electronegative group (Ketone group) has likely contributed to enhanced activity.

3.2.4.2 CoMSIA Contour maps

In Figure 11.a, we observe that the CoMSIA/E contour map is located at R3 position, which indicates that electropositive groups like methyl is favorable (blue color) to the increase of the activity. Similarly, the electronegative (red color) region observed near R4 position provide that the presence of nitrogen atoms (N atom) provide by the heterocyclic rings (imidazole, Triazole, 1-methyl imidazole) are favorable to enhance the is conducive to activities enhancement. Moreover, the favored regions of CoMSIA/HBD, are displayed by cyan color and the unfavorable regions are shown by purple color. Figure 11.b, shows a big cyan region around R4 position suggests that heterocyclic rings like imidazole group may contribute favorably to the HBD interactions, consequently improving the antiproliferative activity. Furthermore, the CoMSIA/HBA presented in Fig 11.c indicates that the favored areas are visualized by magenta color and the red color represents the unfavorable regions to the activity. From Fig 11.c, we can deduce that the ketone group at R1 and R3 positions is favorable to enhance the activity, while the nitrogen atom acceptor is detrimental for the activity.

3.2.4.3 HQSAR contribution maps

In addition to the CoMFA and CoMSIA contour maps, HQSAR provides complete contribution maps showing where molecule fragments could contribute to the antiproliferative activity of ursolic acid derivatives. Figure 12 shows the colored contributions for the highly active compounds 30 (pIC₅₀=5.72) and 26 (pIC₅₀=5.67). The contribution map is codifying by different colors, the extreme green color represents positive contributions for the antiproliferative activity, the extreme red color represents the negative contributions, while the observed white color represents the neutral contributions for the activity.

3.3 Homology-based modeling

The sequence identity between PEBP1 target and the protein (Template PDB ID: 6ENS) has a value of 86.49%, which fulfills the condition proposed by Rost.B et al (identity>25%) [57]. Therefore, the crystal structure of the template is more suitable for homology modeling of PEBP1 sequence. The initial alignment of PEBP1 sequence with the template protein was examined using T-Coffee (http://tcoffee.crg.cat/) (Fig. 13.a). The Quality of PEBP1 model was examined by the Ramachandran plot (Fig. 13.b) and Table A7 (Electronic Supplementary Material, ESI†), indicating that 100% of residues are positioned in the favored regions. In addition, the GMQE score of 0.97 close to 1 expected the accuracy of alignment to build PEBP1 model. In addition, the QMEAN value of -0.10 indicate the best accuracy of the obtained PEBP1 model with the experimental data. Those results indicate that PEBP1 model is stable and reliable to explain the further binding modes of Ursolic Acids using molecular docking.

3.4 Molecular docking studies

Docking simulation is the most approach used to screen the binding modes of the retained Homology PEBP1 model (Fig. 13.c). Hence, it was employed in this study to elucidate the predicted binding pocket using a set a set of compounds of ursolic acid derivatives. To investigate the binding affinity of the ursolic acid hetrocyclic derivatives, we are selected the active compounds 30, 29, 26 (with heterocylic ring) and the less active compounds 1 (without hetrocylic ring). As shown in Table 5, the active compounds 30 (pIC₅₀ = 5.72), 29 (pIC₅₀ = 5.67), 26 (pIC₅₀ = 5.67), and the less active compound 1 (pIC₅₀ = 4.89) showed higher binding energies, which are equal to -10, -9.1, -9.2, and -8.4 Kcal/mol, respectively. The molecule 30 has the best binding energy (-10 Kcal/mol) with the highest antiproliferative activity compared to others compounds. To better understand the interactions in the selected liaison pocket, we have represented graphically the highest active compound 30 to elucidate the different interactions relevant to the molecular docking.

Figure 14.a shows the complex PEBP1-M30 with the dominance of hydrophobic and hydrogen bond interactions. As we see, the hydrophobic interactions have been revealed with the key residues PRO112, TYR181, TYR81, PRO74, LEU184, TRP84 and ALA 73. Meanwhile, the residues TRP173, LYS113, TRP84, GLY141, GLY147, and TRP84 contribute to the hydrogen bond interactions. Overall, residues PRO112, TYR181, TYR81, PRO74, TRP84 were conserved throughout docking simulations. The binding mode of the PEBP1-M30 is very strong and the bond distances are: TRP173 (2.76 Å), LYS113 (2.06 Å), TRP84 (1.90 Å), PRO112 (4.51 Å, 4.18 Å) and 4.03 Å), TYR181 (5.10 Å and 5.19 Å), TYR81 (4.45 Å and 4.93 Å), and PRO74 (4.51 Å, 5.42 Å and 5.37 Å) (see Table 5).

Structurally, residue TRP84 interacts with carbonyl group at R3 position of compound M30, while residues TRP173 and LYS113 interact with the imidazole group at R3 position via hydrogen bond interactions. A strong hydrophobic network is observed around the triterpenoid scaffold by Pi-Pi T-shaped, alkyl and Pi-alkyl interactions. On the other hand, the less active compound (1) shows (Fig. 14.b) no hydrogen bond interactions, which could explain its low antiproliferative activity compared to the active compound 30. These results are in accordance with the 2D-QSAR and 3D-QSAR analysis, which clearly suggest that the heterocycle groups like imidazole contribute favorably to improving the antiproliferative activity. The 2D and 3D represnetations of the docked molecules M26 and M29 are displayed in the ESI† (Fig A1, and A2).

Table 5

Molecular Docking results of the homology PEBP1 model.
Cpds ID	Energies score (Kcal/mol)	Antiproliferative Activities (pIC₅₀)	Amino acid Residues H-Bond interactions	Amino acid Residues Hydrophobic interactions
30	-10	5.72	TRP173 (2.76 Å), LYS113 (2.06 Å), TRP84 (1.90 Å)	PRO112 (4.51 Å, 4.18 Å) and 4.03 Å), TYR181 (5.10 Å and 5.19 Å), TYR81 (4.45 Å and 4.93 Å), PRO74 (4.51 Å, 5.42 Å and 5.37 Å)
29	-9.2	5.67	GLY141 (5.44 Å)	TYR181 (5.13 Å), LEU184 (5.83 Å), TRP84 (4.95 Å, 4.13 Å and 5.81 Å), LEU184 (5.83 Å), ALA 73 (5.06 Å and 5.40 Å), TYR81 (5.10 Å)
26	-9.1	5.67	GLY147 (3.65 Å), TRP84 (5.22)	TYR181 (4.81 Å, 4.28 Å and 5.88 Å), PRO112 (4.07 Å and 5.35 Å), PRO74 (5.28 Å), TYR81 (4.76 Å and 6.69 Å), TRP84 (6.38 Å)
UA (1)	-8.4	4.89	-	TYR81 (4.09 Å, 5.87 Å), TRP84 (5.43 Å, 5.23 Å), PRO (4.08 Å, 4.48 Å)

3.5. Molecular electrostatic potential (MEP) analysis

The importance of MEP lies in its ability to predict binding affinity of new candidate drugs [58]. This approach indicates the intermolecular specificity of molecules, necessary to predict their chemical reactivity and also some types of ligand-protein interactions [59]. Electrostatic molecular potential maps are formed by moving a charged probe across the surface of molecule. At each point, the potential electrostatic energy of the probe is computed and the surface of the molecule is coded accordingly by the specific colors. The blue color indicates the extreme electropositive region susceptible to react with a nucleophilic group and the red color represents the extreme electronegative region likely to react with an electrophile group. The observed medium color regions represent the middle of the charge between both extremes’ cases. In this study, the molecular potential electrostatic maps (MEP) were generated for M1: Ursolic acid, less compound (pIC₅₀ = 4.89) and for the active compound in the series (30: pIC₅₀ = 5.72) using DFT (B3LYP/6-31G(d,p)) level. Analysis of the results shows, in the one hand, that the ESP surfaces of compound 30 shows that extremely negative regions exist around N38, N39, O29, and O26 atoms, while the extremely positive surfaces are near H85, H83, and H42 atoms (see Figure 15.a). Moreover, the observed electronegative maps on atoms N38, N39, and O26 are connected with residues TRP173, LYS113, and TRP84 respectively of the homology model PEBP1 and the electropositive surfaces observed on atoms H85, H83, and H42 grafted with the amino acids of PRO112 and TYR181, which is in agreement with the molecular docking study. On the other hand, The MEP maps of compound 1 shows that extreme electropositive surfaces are around atoms H62 and H69 (Figure 15.b), while the electronegative surfaces are observed around atoms O22 and O27. From this analysis, we can see that heterocyclic Ursolic acids (highest active compound (30)) keep better stability in the active pocket compared to the less active compound 1 (Ursolic acid without heterocycle group), these results are in accordance with CoMFA and CoMSIA analysis.

3.6 Newly predicted compound with potent anticancer activity

Based in the results obtained by 2D-QSAR, 3D-QSAR, Hologram QSAR (HQSAR) and molecular docking simulations, a set of new molecules was predicted based in DFT (B3LYP/6-31G(d,p)) optimization method as depicted in Table 6. Also, the predicted compounds (X1, X2, X3) were subjected to molecular docking simulations to test their stability against modeled PEBP1 protein. The results of these simulations are presented in Table 6. From the analysis of this table, we can notice that all predicted compounds (X1, X2, X3) present a higher binding affinity than that observed for the compound with the highest activity in the series of molecules studied (M30). While, the stability of those predicted compounds was also performed by ulterior Molecular dynamics simulations. The 3D and 2D interactions of the predicted compounds with PEBP1 protein were shown in Fig. A3, A4 and A5 (Supplementary materials).

Table 6: New predicted compounds with their anticancer activities

3.7 Molecualr dynamic simulations (MDs)

3.7.1 Kinetic and potential energies analysis

To elucidate the dynamic interactions of PEBP1 complex with active compound (M30) and the predicted compounds (X1, X2, X3), we performed 50 ps molecular dynamic simulation studies. The kinetic and potential energies plots of studied compounds (M30, X1, X2, X3) with PEBP1 protein (see Fig. 16.a and b) suggested a well-equilibrated complex system through the simulation. It is well understood that in the middle of the simulation process the system has a tendency to be reaches equilibrium for all studied complexes. From this analysis, we noticed that the all studied complexes were stable with no significant difference. Moreover, these stabilities can also prove the higher negative binding energy of the predicted compounds (X1, X2, X3) obtained through molecular docking simulations (Table 6). Due to the fact that the higher the value of binding energy, more stable will be the stability of the ligands on the binding pocket of the protein, we can be seen (Table 6), that the binding energy of predicted compound X2 (E= -11.9 Kcal/mol has better stability on the PEBP1 binding pocket, what is in agreement with the potential energy analysis (Fig. 16b).

3.7.2 Root Mean Square Deviation and Radius of gyration analyses

To evaluate the stability of the studied complexes and characterize the molecular interactions between ligands and PEBP1, we calculated the root-mean-square deviation (RMSD) through 50 ps of simulations. RMSD calculation was based upon the complex backbone atoms over 50-ps dynamic simulation. RMSD graph of PEBP1, M30, X1, X2, clearly shows a slight deviation in conformation by 1.1 to 1.6 nm after 30 ps of simulation (Figure 16c). While the complex PEBP1-X3 show the lower RMSD value around 0.6 nm during the 50 ps of MD simulations. Comparing the RMSD graphs of the four binding complexes reveals that, complexes with predicted ligands (X1, X2, X3) tends to be well stabilized and equilibrated complexes than active compound (M30). In addition, results of the present study do not show any larger conformational changes in the structure of the PEBP1 protein when subjected to Molecular dynamic simulations.

Radius of gyration (Rg) was determined to understand the level of compaction in the structure of PEBP1 in the presence of three hits (X1, X2, X3) and active compound in the series (M30). The Rg is defined as the mass weighted root mean square distance of a collection of atoms from their common center of mass [60]. As depicted in Fig. 16d, the Rg values of PEBP1-X1, PEBP1-X2, PEBP1-X3 and PEBP1-M30 varies between 1.75 nm and 2 nm (Figure 16d). It shows that PEBP1-X3 has marginally lower Rg values, however the difference cannot be considered significant.

In this work, we studied the anticancer activity of a series of 30 molecules of Ursolic Acid derivatives using the semi-empirical method AM1 and the empirical method DFT. These methods have allowed us to build predictive 2D, 3D, and Hologram QSAR models, that link biological activities to molecular descriptors, using statistical analysis. The retained 2D-QSAR results proved that DFT (B3LYP(d, p)) method is more accurate than AM1 method with a high Q2 value of 0.81 and a maximum external validity R2pred value of 0.78. These results suggest that the DFT (B3LYP(d, p)) method is the most appropriate approach to predict the antiproliferative activities of Ursolic acid derivatives. Thus, the same quantum level were used to select the appropriate molecular alignment, giving the best 3D and Hologram QSAR models. The results found that molecular optimization by DFT (B3LYP(d, p)) is more appropriate to build the 3D-QSAR model, however, the contour maps of the CoMFA and CoMSIA analysis proved that steric, electrostatic, and hydrogen bond donor/acceptor fields are beneficial to the antiproliferative activity. Moreover, we found no significant change in the performance of the HQSAR model for the two theoretical levels. Homology modeling was investigated to build the 3D structure of target protein PEBP1 and then validated it by Ramachandran plots and Qmean score for further molecular docking studies. As a result, the best binding affinities between the active Ursolic acid derivatives against the Modeled PEBP1 protein were obtained, then validated by the molecular potential maps to select each atomic contribution responsible to the observed binding affinity. Conseduently, the above studies, allowed us to predict three new triterpenoids derivatives with potent anticancer activity. In addition, the molecular docking and molecular dynamics silmulations were suggested the high conformational stability of the predicted compounds against the modeled PEBP1 protein.

hPEBP1: Human Phosphatidylethanolamine Binding Protein 1

QSARs: Quantitative structure activity relatioships

2D: two dimentional

3D: tree dimentional

DFT: Density Functional Theory

AM1: Austin Model 1

MDs: Molecular Dynamics

CoMFA: Comparative molecular field analysis

CoMSIA: Comparative molecular similarity indices analysis

HQSAR: Hologram Quantitative structure activity relatioships

GMQE score: Global Model Quality Estimation

OECD: Organisation for Economic Co-operation and Development.

Acknowledgements:

Realized with the support of the National Center for Scientific and Technical Research (CNRST - Morocco) as part of the Research Excellence Awards Program (No. 34USMBA2017).

Funding: This work was realized with the support of Sidi Mohamed Ben Abdellah University of Fez - Morocco

Conflicts of interest/Competing interests: The authors declare no competing interests.

Availability of data and material: The Data is intergated in the Manuscript.

Code availability: Not applicable

Authors' contributions:

MS did most of the practical work as part of a PhD thesis supervised by HT and prepared the manuscript. HT designed and coordinated the study, participated in article preparation, corrected the manuscript and edited the final version and submitted it for publication. TA contributed to data analysis. FL participated in study designed, helped to improve the manuscript and critically revised the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate Compliance with ethical standards.

Bray F, Ferlay J, Soerjomataram I et al (2018) Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer. Journal for Clinicians 68:394–424. https://doi.org/10.3322/caac.21492
Schoentgen F, Jonic S (2018) PEBP1/RKIP: from multiple functions to a common role in cellular processes. arXiv preprint arXiv:1802.02378
Beshir AB, Ren G, Magpusao AN et al (2010) Raf kinase inhibitor protein suppresses nuclear factor-κB-dependent cancer cell invasion through negative regulation of matrix metalloproteinase expression. Cancer Lett 299:137–149. https://doi.org/10.1016/j.canlet.2010.08.012
Al-Mulla F, Bitar MS, Al-Maghrebi M et al (2011) Raf Kinase Inhibitor Protein RKIP Enhances Signaling by Glycogen Synthase Kinase-3β. Cancer Res 71:1334–1343. https://doi.org/10.1158/0008-5472.CAN-10-3102
Rath O, Park S, Tang H et al (2008) The RKIP (Raf-1 Kinase Inhibitor Protein) conserved pocket binds to the phosphorylated N-region of Raf-1 and inhibits the Raf-1-mediated activated phosphorylation of MEK. Cell Signal 20:935–941. https://doi.org/10.1016/j.cellsig.2008.01.012
Yan X-F, Xiao H-M, Gong X-D, Ju X-H (2006) A comparison of semiempirical and first principle methods for establishing toxicological QSARs of nitroaromatics. J Mol Struct (Thoechem) 764:141–148. https://doi.org/10.1016/j.theochem.2006.02.018
Puzyn T, Leszczynski J, Cronin MTD (2010) Recent advances in QSAR studies: methods and applications. Challenges and advances in computational chemistry and physics. Springer, Dordrecht; New York
Cronin MTD, Livingstone D (2004) Predicting chemical toxicity and fate. CRC Press
Kubinyi H (1998) Comparative Molecular Field Analysis (CoMFA). The encyclopedia of computational chemistry, 1, 448-460
Ghemtio L, Zhang Y, Xhaard H (2012) CoMFA/CoMSIA and pharmacophore modelling as a powerful tools for efficient virtual screening: Application to anti-leishmanial betulin derivatives. In: Virtual Screening, 55–82
Lowis DR (1997) HQSAR: A New, Highly Predictive QSAR Technique. Tripos technical notes 1:17
Muhammed MT, Aki-Yalcin E (2019) Homology modeling in drug discovery: Overview, current applications, and future perspectives. Chem Biol Drug Des 93:12–20. https://doi.org/10.1111/cbdd.13388
Safarizadeh H, Garkani-Nejad Z (2019) Molecular docking, molecular dynamics simulations and QSAR studies on some of 2-arylethenylquinoline derivatives for inhibition of Alzheimer’s amyloid-beta aggregation: Insight into mechanism of interactions and parameters for design of new inhibitors. J Mol Graph Model 87:129–143. https://doi.org/10.1016/j.jmgm.2018.11.019
Muneeswaran G, Pandiaraj M, Kartheeswaran S et al (2018) Molecular dynamics simulation approach to explore atomistic molecular mechanism of peroxidase activity of apoptotic cytochrome c mutants. Informatics in Medicine Unlocked 11:51–60. https://doi.org/10.1016/j.imu.2018.04.003
Woźniak Ł, Skąpska S, Marszałek K (2015) Ursolic Acid—A Pentacyclic Triterpenoid with a Wide Spectrum of Pharmacological Activities. Molecules 20:20614–20641. https://doi.org/10.3390/molecules201119721
Mourya A, Akhtar A, Ahuja S et al (2018) Synergistic action of ursolic acid and metformin in experimental model of insulin resistance and related behavioral alterations. Eur J Pharmacol 835:31–40. https://doi.org/10.1016/j.ejphar.2018.07.056
López-Hortas L, Pérez-Larrán P, González-Muñoz MJ et al (2018) Recent developments on the extraction and application of ursolic acid. A review. Food Res Int 103:130–149. https://doi.org/10.1016/j.foodres.2017.10.028
Hua S-X, Huang R-Z, Ye M-Y et al (2015) Design, synthesis and in vitro evaluation of novel ursolic acid derivatives as potential anticancer agents. Eur J Med Chem 95:435–452. https://doi.org/10.1016/j.ejmech.2015.03.051
Tang Q, Liu Y, Li T et al (2016) A novel co-drug of aspirin and ursolic acid interrupts adhesion, invasion and migration of cancer cells to vascular endothelium via regulating EMT and EGFR-mediated signaling pathways: multiple targets for cancer metastasis prevention and treatment. Oncotarget 7:. https://doi.org/10.18632/oncotarget.12232
Mlala S, Oyedeji AO, Gondwe M, Oyedeji OO (2019) Ursolic Acid and Its Derivatives as Bioactive Agents. Molecules 24:2751. https://doi.org/10.3390/molecules24152751
Qian Z, Wang X, Song Z et al (2015) A Phase I Trial to Evaluate the Multiple-Dose Safety and Antitumor Activity of Ursolic Acid Liposomes in Subjects with Advanced Solid Tumors. Biomed Res Int 2015:1–7. https://doi.org/10.1155/2015/809714
Sultana N (2011) Clinically useful anticancer, antitumor, and antiwrinkle agent, ursolic acid and related derivatives as medicinally important natural product. J Enzyme Inhib Med Chem 26:616–642. https://doi.org/10.3109/14756366.2010.546793
Iqbal J, Abbasi BA, Ahmad R et al (2018) Ursolic acid a promising candidate in the therapeutics of breast cancer: Current status and future implications. Biomed Pharmacother 108:752–756. https://doi.org/10.1016/j.biopha.2018.09.096
Leal AS, Wang R, Salvador JAR, Jing Y (2012) Synthesis of novel ursolic acid heterocyclic derivatives with improved abilities of antiproliferation and induction of p53, p21waf1 and NOXA in pancreatic cancer cells. Bioorg Med Chem 20:5774–5786. https://doi.org/10.1016/j.bmc.2012.08.010
Gaussian 09 (2009) R.A.: 1, mj frisch, gw trucks, hb schlegel, ge scuseria, ma robb, jr cheeseman, g. Scalmani, v. Barone, b. Mennucci, ga petersson et al., gaussian. Inc Wallingford CT. 121, 150
Kohn W, Sham LJ (1965) Self-Consistent Equations Including Exchange and Correlation Effects. Phys Rev 140:A1133–A1138. https://doi.org/10.1103/PhysRev.140.A1133
Lazrak M, Toufik H, Bouzzine SM et al (2018) The computational study of bridge effect in D-π-A photosensitive dyes, based on triphenylamine. IOP Conf Ser: Earth Environ Sci 161:012021. https://doi.org/10.1088/1755-1315/161/1/012021
Lazrak M, Toufik H, Bouzzine SM, Lamchouri F (2020) Bridge effect on the charge transfer and optoelectronic properties of triphenylamine-based organic dye sensitized solar cells: theoretical approach. Res Chem Intermed 46:3961–3978. https://doi.org/10.1007/s11164-020-04184-x
Casida ME (1995) Time-Dependent Density Functional Response Theory for Molecules. In: Recent Advances in Computational Chemistry. WORLD SCIENTIFIC 155–192. https://doi.org/10.1142/9789812830586_0005
Ennehary S, Toufik H, Bouzzine SM, Lamchouri F (2020) Effect of the alkyl chain length on the optoelectronic properties of organic dyes: theoretical approach. J Comput Electron 19:840–848. https://doi.org/10.1007/s10825-020-01486-6
Gramatica P, Chirico N, Papa E et al (2013) QSARINS: A new software for the development, analysis, and validation of QSAR MLR models. J Comput Chem 34:2121–2132. https://doi.org/10.1002/jcc.23361
Stitou M, Toufik H, Akabli T et al (2018) 2D-QSAR method of lupane-type saponins the treatement of cancer cell line. RHAZES: Green Applied Chemistry 2:33–45
Akabli T, Toufik H, Yasri A et al (2018) Combining ligand-based and structure-based drug design approaches to study the structure-activity relationships of a β-carboline derivative series. Struct Chem 29:1637–1645. https://doi.org/10.1007/s11224-018-1141-1
O’brien RM (2007) A Caution Regarding Rules of Thumb for Variance Inflation Factors. Qual Quant 41:673–690. https://doi.org/10.1007/s11135-006-9018-6
Golbraikh A, Tropsha A (2000) Predictive QSAR modeling based on diversity sampling of experimental datasets for the training and test set selection. Molecular diversity. 231–243. https://doi.org/10.1023/A:1021372108686
ReenuVikas (2016) Evaluating the role of electron-correlation in the external prediction of the toxicity of nitrobenzenes towards Tetrahymena pyriformis. New J Chem 40:2343–2353. https://doi.org/10.1039/C5NJ02552D
Stitou M, Toufik H, Bouachrine M et al (2019) Machine learning algorithms used in Quantitative structure-activity relationships studies as new approaches in drug discovery. In: 2019 International Conference on Intelligent Systems and Advanced Computing Sciences (ISACS), IEEE, Taza, Morocco, 1–8. https://doi.org/10.1109/ISACS48493.2019.9068917
Rücker C, Rücker G, Meringer M (2007) y-Randomization and Its Variants in QSPR/QSAR. J Chem Inf Model 47:2345–2357. https://doi.org/10.1021/ci700157b
Roy K, Kar S, Ambure P (2015) On a simple approach for determining applicability domain of QSAR models. Chemometr Intell Lab Syst 145:22–29. https://doi.org/10.1016/j.chemolab.2015.04.013
Pourbasheer E, Aalizadeh R, Shokouhi Tabar S et al (2014) 2D and 3D Quantitative Structure–Activity Relationship Study of Hepatitis C Virus NS5B Polymerase Inhibitors by Comparative Molecular Field Analysis and Comparative Molecular Similarity Indices Analysis Methods. J Chem Inf Model 54:2902–2914. https://doi.org/10.1021/ci500216c
Hadni H, Elhallaoui M (2020) 2D and 3D-QSAR, molecular docking and ADMET properties in silico studies of azaaurones as antimalarial agents. New J Chem 44:6553–6565. https://doi.org/10.1039/C9NJ05767F
Halim SA, Zaheer-ul-Haq (2015) Structure based 3D-QSAR studies of Interleukin-2 inhibitors: Comparing the quality and predictivity of 3D-QSAR models obtained from different alignment methods and charge calculations. Chem Biol Interact 238:9–24. https://doi.org/10.1016/j.cbi.2015.05.018
SYBYL/QSAR. Molecular Modelling Software, Tripos Inc., 1699 S. Hanley Road, St. Louis, MO 63944, USA
Stitou M, Toufik H, Bouachrine M, Lamchouri F (2020) Quantitative structure–activity relationships analysis, homology modeling, docking and molecular dynamics studies of triterpenoid saponins as Kirsten rat sarcoma inhibitors. Journal of Biomolecular Structure Dynamics 39:152–170. https://doi.org/10.1080/07391102.2019.1707122
Matysiak J, Niewiadomy A (2017) QSAR models of antiproliferative activity of imidazo[2,1-b][1,3,4]thiadiazoles in various cancer cell lines. Mol Diversity 21:211–218. https://doi.org/10.1007/s11030-016-9705-8
Wang X, Yan J, Wang M et al (2018) Synthesis and three-dimensional quantitative structure-activity relationship study of quinazoline derivatives containing a 1,3,4-oxadiazole moiety as efficient inhibitors against Xanthomonas axonopodis pv. citri. Mol Diversity 22:791–802. https://doi.org/10.1007/s11030-018-9837-0
Arnold K, Bordoli L, Kopp J, Schwede T (2006) The SWISS-MODEL workspace: a web-based environment for protein structure homology modelling. Bioinformatics 22:195–201. https://doi.org/10.1093/bioinformatics/bti770
Ye W-L, Zhang L-X, Guan Y-D et al (2019) Virtual screening and experimental validation of eEF2K inhibitors by combining homology modeling, QSAR and molecular docking from FDA approved drugs. New J Chem 43:19097–19106. https://doi.org/10.1039/C9NJ02600B
Benkert P, Tosatto SCE, Schomburg D (2008) QMEAN: A comprehensive scoring function for model quality assessment. Proteins: Struct Funct Bioinf 71:261–277. https://doi.org/10.1002/prot.21715
Paulsen JL, Anderson AC (2009) Scoring ensembles of docked protein: ligand interactions for virtual lead optimization. J Chem Inf Model 49:2813–2819. https://doi.org/10.1021/ci9003078
Chaube U, Bhatt H (2017) 3D-QSAR, molecular dynamics simulations, and molecular docking studies on pyridoaminotropanes and tetrahydroquinazoline as mTOR inhibitors. Mol Diversity 21:741–759. https://doi.org/10.1007/s11030-017-9752-9
Walton IM, Cox JM, Benson CA et al (2016) The role of atropisomers on the photo-reactivity and fatigue of diarylethene-based metal–organic frameworks. New J Chem 40:101–106. https://doi.org/10.1039/C5NJ01718A
Frimand K, Jalkanen KJ (2002) SCC-TB, DFT/B3LYP, MP2, AM1, PM3 and RHF study of ethylene oxide and propylene oxide structures, VA and VCD spectra. Chem Phys 279:161–178. https://doi.org/10.1016/S0301-0104(02)00457-3
Pasha FA, Srivastava HK, Singh PP (2005) Comparative QSAR study of phenol derivatives with the help of density functional theory. Bioorg Med Chem 13:6823–6829. https://doi.org/10.1016/j.bmc.2005.07.064
Deeb O, Clare BW (2008) Comparison of AM1 and B3LYP-DFT for Inhibition of MAO-A by Phenylisopropylamines: A QSAR Study. Chem Biol Drug Design 71:352–362. https://doi.org/10.1111/j.1747-0285.2008.00643.x
Myint KZ, Xie X-Q (2010) Recent Advances in Fragment-Based QSAR and Multi-Dimensional QSAR Methods. IJMS 11:3846–3866. https://doi.org/10.3390/ijms11103846
Rost B, Sander C (1996) Bridging the protein sequence-structure gap by structure predictions. Annu Rev BioPhys BioMol Struct 25:113–136
Weiner PK, Langridge R, Blaney JM et al (1982) Electrostatic potential molecular surfaces. Proceedings of the National Academy of Sciences 79:3754–3758. https://doi.org/10.1073/pnas.79.12.3754
Rathi PC, Ludlow RF, Verdonk ML (2019) Practical High-Quality Electrostatic Potential Surfaces for Drug Discovery Using a Graph-Convolutional Deep Neural Network. Journal of medicinal chemistry 63:8778–8790. https://doi.org/10.1021/acs.jmedchem.9b01129
Monajjemil M, Oliaey AR (2009) Gyration Radius and Energy Study at Different Temperatures for Acetylcholine Receptor Protein in Gas Phase by Monte Carlo, Molecular and Langevin Dynamics Simulations. Journal of Physical & Theoretical Chemistry of Islamic Azad University of Iran. 5:195–201

Download PDF

Editorial decision: Major revisions needed
27 Nov, 2021
Reviews received at journal
22 Oct, 2021
Reviewers invited by journal
21 Oct, 2021
Editor invited by journal
21 Oct, 2021
Editor assigned by journal
21 Oct, 2021
First submitted to journal
19 Oct, 2021

You are reading this latest preprint version

Virtual Screening of PEBP1 Inhibitors By Combining 2D/3D-QSAR Analysis, Hologram QSAR, Homology Modeling, Molecular Docking Analysis and Molecular Dynamic Simulations

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials & Methods

2.1 2D-QSAR modeling

2.1.1 Dataset treatment

2.1.2 Computational details

2.1.3 Statistical testing

2.1.3.1 Partial Least Squares method (PLS)

2.1.3.2 2D-QSAR model assessment

2.1.3.3 Applicability domain (AD)

2.2. 3D-QSAR Methods

2.2.1. Molecular alignment

2.2.2 CoMFA strategy

2.2.3 CoMSIA strategy

2.2.4. Partial least squares (PLS) analysis and 3D-QSAR models assessment

2.3 Hologram QSAR (HQSAR)

2.4 Homology modeling

2.4.1 Templates searching and modeled PEBP1 protein quality

2.5 PEBP1 model -Ursolic Acid Docking

2.6 Molecular Dynamic simulation and system building

3. Results And Discussions

3.1. 2D-QSAR analysis

3.1.1 Selection of descriptors

3.1.2 PLS regression

3.1.3 Predictive powers of retained 2D-QSAR models

3.1.4 Applicability domain (AD)

3.2 3D-QSAR analysis

3.2.1 CoMFA results

3.2.2. CoMSIA results

3.2.3 HQSAR results

3.2.4 Contour maps of CoMFA, CoMSIA, and HQSAR

3.2.4.1 CoMFA contour plots

3.2.4.2 CoMSIA Contour maps

3.2.4.3 HQSAR contribution maps

3.3 Homology-based modeling

3.4 Molecular docking studies

3.5. Molecular electrostatic potential (MEP) analysis

3.6 Newly predicted compound with potent anticancer activity

3.7 Molecualr dynamic simulations (MDs)

3.7.1 Kinetic and potential energies analysis

3.7.2 Root Mean Square Deviation and Radius of gyration analyses

4 Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1