DOI: https://doi.org/10.21203/rs.3.rs-1233496/v1
The p53 is a tumor suppressor protein that adjusts cell cycle and growth arrest as well as genes that restore DNA damage and apoptosis. Murine double minute 2 (MDM2) is a main p53 antagonist. We created a novel QSAR model using a series of highly active spiro [pyrrolidin-3,2-oxindoles] that consisted of 29 compounds that were experimentally validated to inhibit the MDM2-p53 interaction. Three optimal models have been developed CoMFA/E+S, CoMSIA/S+H+A and HQSAR have revealed good statistical results, but the CoMSIA mode only which validates all the external validation tests applied successfully. Based on the CoMSIA/S+H+A model was carefully chosen to design four compounds with values of inhibitory activity greater than the highly active compound in the data set. The Newly designed compounds were docked in the target receptor binding site (ID: 4LWU). The newly designed compound Pred 01 showed the highest affinity with a value of -9.4kcal / mol, while compound N°04 which represents the data set and control compound (Nutlin-3) showed binding energies of the order of -8.8 kcal / mol and -8.2kcal / mol, respectively. In addition, the roles of lipinski and veber were estimated, the results obtained demonstrate that the proposed molecules involve good oral bioavailability and an ability to diffuse through different biological barriers. For in-depth study, The Pred01 / receptor, N°04 / receptor and Nutlin-3 / receptor complexes were selected via dynamic simulation analyzes with a simulation time of 100 ns and, also their free binding energy was examined operating the MM-GBSA approach. The molecular docking results obtained accentuate the crucial residues responsible for the ligand / protein interaction, providing insight into the mode of interaction. The MD simulation analysis confirms the conformational stability of the selected complexes during the MD trajectory, and the fluctuations recorded are insignificant. The results of MM-GBSA reveal that the new compound Pred 01 exhibits the lowest free energy, which confirms the result of molecular docking.
The term ferroptosis was first coined in 2012 to characterize a form of programmed cell death induced by the accumulation of lethal iron-dependent lipid peroxidation [1]. Ferroptosis can be affected by nutritional factors, pharmacological agents, and physiological functions that directy affect these factors. Its may have potential uses in the treatment of cancer, induction of ferroptosis in vivo may have therapeutic value. Multiple pathways can initiate ferroptosis, most of which act directly through the mechanism of using covalent inhibitors to bind glutathione peroxidase 4 (GPX4). It is a uniquely designed selenoprotein that is able to reduce lipid hydroperoxides (L-OOH) in fatty alcohols (L-OH) and thus neutralize their toxicity[2, 3].
In recent years, some studies have shown that the p53 protein may be involved in the control of ferroptosis through the accumulation of lipid hydroperoxides or toxic iron [4–6]. The p53 transcription factor is one of the most usually mutated in cancer, and approximately 50% of human cancers are linked to its mutations or deletions [7, 8]. In 2015, the first report of p53 response to ferroptosis was published, the authors discovered that p53 promotes ferroptosis in fibroblasts and some cancer cells that are similar to human breast cancer MCF7[9].
On the other hand, the activity of p53 can be regulated by overexpression of double minute mice 2 (MDM2) and may induce a tendency for tumor growth [10]. Hence, the use of a small molecule inhibitor to reactivate p53 function by disrupting MDM2-p53 interactions is now considered a promising strategy for cancer treatment. BDP (benzodiazepinediones) and nutlines have been identified as the most potent p53-MDM2 interaction abrogators so far Fig. 1[11].
In recent years, the combination of spirooxindole and p pyrrolidine has garnered substantial interest. McConnell et al. developed a series of spirooxindoles fused with 3,2-pyrrolidine, which promotes the formation of functional groups to form p-p interactions and hydrogen bond types and provided great properties to these compounds [12]. The mentioned properties of the spirooxindole scaffold fused to 3,2-pyrrolidine reveal the possibility of introducing it synchronously to construct a novel dual GPX4 / MDM2 inhibitor. New strategies for designing novel molecules against protein targets are currently in vogue for the development of anticancer treatments.
The QSAR protocol is an process to quantifie a relationship between activity of inhibitors and the structure using mathematical and computer tools for predict the activity biologic of new compounds [13, 14]. In the present research, we performed 2D / 3D QSAR analysis on series of nitroisoxazole derivatives as MDM2-p53 antagonists. In the case of 3D / QSAR modeling, the CoMFA (comparative molecular field analysis) methods and CoMSIA (comparative molecular similarity index analysis) methods were applied to establish the 3D / QSAR model. In addition, a 2D model has been established using the Halogram QSAR method, which is a recent method which. On the contrary 3D-QSAR, avoids the calculation of the molecular descriptors of the compound and it takes a long time therefore this method does not take into account the complexations of the 3D structure and the alignment of molecules to predict their activity. The best model built was based to discover new molecules to disrupt p53-HDM2 interactions. Moreover these proposed molecules, the most active compound and the control compound Nutlin-3 were docked at the MDM2 protein binding site to understand the mode of ligand / protein interaction as well as to identify the amino acids that govern activity. An MD simulation analysis were performed to access the stability of the complexes studied and validate the molecular docking results. The ADME / Tox properties of the modified compounds were also predicted to confirm the biodispobinility of the populated molecules as a drug candidate.
In this research, we selected a database of 29 compounds nitroisoxazole containing spiro [pyrrolidine-oxindole] derivatives from the literature[15]. The basic structure used is replaced by an isoxazole in a small molecule inhibiting GPX4 for the development of new anti-tumor drugs. A fluorine substitution has also been found in the basic structure studied, where it reveals an essential role in improving the efficacy of drugs which has become widespread in recent decades[17, 18]. In particular at the level of the optimization of the molecular properties of a drug and of influencing the selectivity thereof. target by modulating its conformation, permeability, hydrophobicity, metabolism and lipophilicity[19, 29]. In addition, in 2019, eleven fluoride-containing drugs were authorized by the FDA, more than half of which contain a trifluoromethyl group[20]. This suggested to us that the strategic integration of a trifluoromethyl group can enhance the inhibitory activity of a drug to some level.
The database structures were designed and optimized and again saved in mol2 format utilizing SYBYL X2.0 software[21]. Biological activity data of nitroisoxazole derivatives in these literatures are expressed as Ki (10-6 M), and we have converted them to logarithmic pki (-Log ki). The structures and activity (pKi) used in this work for were collected in Table 1. The 29 compounds were haphazardly divided into 23 compounds for the training set and 6 compounds for the test set in order to facilitate the 3D-QSAR analysis[22].
Compound N° |
R1 |
R2 |
R3 |
R4 |
pKi exp |
CoMFA |
CoMSIA |
HQSAR |
||||
(10−6M) |
pKi pred |
Residu |
pKi pred |
Residu |
pKi pred |
Residu |
||||||
1 |
Train |
H |
H |
Ph |
CF3 |
5.893 |
5.890 |
0.003 |
5.919 |
-0.026 |
5.882 |
0.011 |
2 |
Train |
5-F |
H |
Ph |
CF3 |
6.309 |
6.458 |
-0.149 |
6.371 |
-0.062 |
5.998 |
0.311 |
3 |
Train |
5-Cl |
H |
Ph |
CF3 |
6.167 |
6.138 |
0.029 |
6.124 |
0.043 |
5.998 |
0.169 |
4 |
Train |
6-Cl |
H |
Ph |
CF3 |
6.620 |
6.484 |
0.136 |
6.549 |
0.071 |
5.921 |
0.699 |
5 |
Train |
5-Br |
H |
Ph |
CF3 |
5.991 |
5.990 |
0.001 |
5.981 |
0.01 |
5.998 |
-0.007 |
6 |
Test |
6-Br |
H |
Ph |
CF3 |
6.585 |
5.570 |
1.015 |
6.41 |
0.175 |
5.921 |
0.664 |
7 |
Test |
5-NO2 |
H |
Ph |
CF3 |
5.434 |
5.892 |
-0.458 |
5.46 |
-0.026 |
5.841 |
-0.407 |
8 |
Train |
5-Me |
H |
Ph |
CF3 |
5.848 |
5.886 |
-0.038 |
5.901 |
-0.053 |
5.999 |
-0.151 |
9 |
Train |
H |
H |
3-F-C6H4 |
CF3 |
5.943 |
5.942 |
0.001 |
5.941 |
0.002 |
5.895 |
0.048 |
10 |
Test |
H |
H |
4-F-C6H4 |
CF3 |
5.987 |
5.597 |
0.39 |
5.824 |
0.163 |
5.79 |
0.197 |
11 |
Train |
H |
H |
4-Cl-C6H4 |
CF3 |
5.558 |
5.543 |
0.015 |
5.534 |
0.024 |
5.79 |
-0.232 |
12 |
Train |
H |
H |
2-Br-C6H4 |
CF3 |
5.000 |
5.025 |
-0.025 |
4.998 |
0.002 |
4.969 |
0.031 |
13 |
Test |
H |
H |
3-Br-C6H4 |
CF3 |
5.556 |
5.863 |
-0.307 |
5.459 |
0.097 |
4.969 |
0.587 |
14 |
Train |
H |
H |
4-Br-C6H4 |
CF3 |
5.528 |
5.502 |
0.026 |
5.523 |
0.005 |
5.79 |
-0.262 |
15 |
Train |
H |
H |
4-NO2-C6H4 |
CF3 |
5.171 |
5.165 |
0.006 |
5.213 |
-0.042 |
5.606 |
-0.435 |
16 |
Train |
H |
H |
2-Me-C6H4 |
CF3 |
5.000 |
5.09 |
-0.09 |
5.047 |
-0.047 |
4.969 |
0.031 |
17 |
Test |
H |
H |
4-Me-C6H4 |
CF3 |
5.827 |
5.59 |
0.237 |
6.069 |
-0.242 |
5.790 |
0.037 |
18 |
Train |
H |
H |
2-OMe-C6H4 |
CF3 |
5.000 |
5.02 |
-0.02 |
5.131 |
-0.131 |
4.732 |
0.268 |
19 |
Train |
H |
H |
4-iPr-C6H4 |
CF3 |
5.646 |
5.663 |
-0.017 |
5.694 |
-0.048 |
5.69 |
-0.044 |
20 |
Train |
H |
H |
2,5-Cl2-C6H3 |
CF3 |
5.000 |
5.017 |
-0.017 |
4.967 |
0.033 |
4.996 |
0.004 |
21 |
Test |
H |
H |
3,4-Cl2-C6H3 |
CF3 |
5.627 |
6.03 |
-0.403 |
5.294 |
0.333 |
5.717 |
-0.09 |
22 |
Train |
H |
H |
3,4-(OMe)2-C6H3 |
CF3 |
5.434 |
5.423 |
0.011 |
5.321 |
0.113 |
5.477 |
-0.043 |
23 |
Train |
H |
H |
2-furyl |
CF3 |
5.924 |
5.908 |
0.016 |
5.873 |
0.051 |
5.978 |
-0.054 |
24 |
Train |
H |
H |
2-thienyl |
CF3 |
5.866 |
5.857 |
0.009 |
5.908 |
-0.042 |
5.978 |
-0.112 |
25 |
Train |
H |
H |
α-naphthyl |
CF3 |
5.078 |
4.97 |
0.108 |
4.984 |
0.094 |
5.009 |
0.069 |
26 |
Train |
H |
Me |
Ph |
CF3 |
5.000 |
4.996 |
0.004 |
4.988 |
0.012 |
5.038 |
-0.038 |
27 |
Train |
H |
allyl |
Ph |
CF3 |
5.000 |
4.979 |
0.021 |
4.994 |
0.006 |
5.281 |
-0.281 |
28 |
Train |
H |
Ac |
Ph |
CF3 |
5.000 |
4.979 |
0.021 |
4.985 |
0.015 |
4.987 |
0.013 |
39 |
Train |
H |
H |
Ph |
Ph |
5.199 |
5.252 |
-0.053 |
5.23 |
-0.031 |
5.980 |
-0.781 |
The structures of the selected compounds were minimized using the conjugate gradient under the standard Tripos Field force field implanted in SYBYL-X2.0 [23] with partial atomic charges of Gasteiger-Hückel by the method of Powell with a convergence criterion of 0.05 Kcal / mol. The maximum iterations number is 10.000. Molecular alignment is an essential stage in 3D-QSAR for associated the orientation of molecules in 3D space. The highest activity compound N°04 was used as a template for alignment. The result of the structural alignment of all inhibitor molecules and the template compound structure is presented in Fig. 2.
3D-QSAR study included the usage of CoMFA and CoMSIA, They are based on the idea that a compound of similar structure interacting with the receptor site in the same mode implies that the molecular fields involving them must also be similar, which shows that the biological activity of the compounds is intimately connected to the molecular field that neighborhoods them. Using the Lenard-Jones function to compute the steric field (S), as well as the Coulomb function to determine the electrostatic field (E), the CoMFA model focuses primarily on van der Waals interactions[24].
Since most of interactions are very complex, and sometimes the steric field and the electrostatic field are insufficient to build the model. For this reason, we had to integrate the CoMSIA method[25]. The essential principle of the CoMSIA model is similar to CoMFA, with the addition of a hydrophobic field (H), a hydrogen bond donor field (D) and a hydrogen bond acceptor field (A). The grid 3D used the spacing grid 2.0Å with the three Cartesian directions. The Gaussian function is used to calculate the energy of the field molecular. The attenuation factor was set at 0.3.
The relation between the activity and the descriptors CoMFA, CoMSIA and HQSAR was quantified by the method of the partial least squares (PLS) integrated in SYBYL X2.0 [26]. We carried out a Leave-one-out (LOO) to select the best value of validation coefficient cross Q² and the best principal component number (NOC). The value of the correlation coefficient (R²), the value F, the standard deviation of estimate (SEE) and the percentage of contribution were selected by non-validation. The 3D-QSAR model is predictive and reliable if Q²> 0.5, r2> 0.6, the SEE value should be as little as possible and the F value should be as large as possible [27].
The HQSAR method is a novel two-dimensional QSAR method that provides an unequaled advantage in molecular design. HQSAR can certainly display where optimization is required, as well as which groups of atoms contribute positively or negatively to molecular activity and find the appropriate structural substitutions in molecule design. The fragment parameters and fragment size were modified during the optimization technique of the HQSAR model. The fragmentation parameter defines how the hologram's topology information is mapped. The atomic number (A) indicates the type of the atom; The type of bond (B) indicates the chemical bond formed between atoms; Atomic connection (C) distinguishes the chemical bond in the fragment; hydrogen atom (H) represents hydrogen atom of the fragment; Chirality (Ch) represents stereochemical information about atomic chirality; hydrogen bond donor / acceptor (DA) signifies the hydrogen bond donor or acceptor of the fragment. In this method we chosen 12 holographic lengths (53, 59, 61, 71, 83, 97, 151, 199, 257, 307, 353 and 401) to build the different HQSAR models[28].
The selected 3D-QSAR and HQSAR models must also be validated by external validation by using the test set reserved, The external validation was carried out according to Golbraikh-Tropsha criteria and using Roy metrics [30–32]. The conditions that the 3D-QSAR models and the HQSAR model must meet certain conditions for a more reliable external predictive capacityare as follows:
- Q2> 0.5, R2> 0.6, (r2 - r20) / r2 <1, (r² - r’²0) /r² <0.1, 0.85 ≤ k’ ≤1.15, 0.85 ≤k ≤ 1.15.
Where:
k and k′ are respectively the slope of the predicted values to those observed and the observed values to those predicted for the test set with zero intercepts/
- Δrm² <0. 2, Δr’m² <0. 2, ȓ²m> 0.5, in Roy's criteria.
Where
- r2 is the square correlation coefficient of the predicted and observed activity values of the test phase.
- r20 and r’20 are respectively the coefficients for determining predicted values and those observed and the observed values and those predicted for the test set with zero interception.
In this work, molecular docking was used to understand the mechanism of interactions between nitroisoxazole derivatives and the chosen receptor code was obtained from the RCSB protein database (https : //www.rcsb.org/) PDB: 4LWU. The AutoDock program performed all molecular docking steps[32]. Within a week of removing all water molecules. The polar hydrogen atoms and Gasteiger charges were added to the receptor protein chain. A grid box spacing of 0.708 Å was taken, with a number point in x, y, and z-dimensions of -11.639, 21.472, and -4.528, respectively. Molecular docking generated nine conformations based on docking affinity, the conformation with the lowest binding energy was then selected.
Furthermore, re-docking of X-ray bound ligand (20U) into their own active site was executed to validate the accuracy of the molecular docking protocol. A molecular docking protocol is reliable if the root mean square deviation (RMSD) is less than 2Å[33].
The Desmond program, a 2005 explicit solvent MD package with the OPLS fixed force field, was used to run MD simulations on the best docked ligand-protein complexes[34]. To simulate the system, an orthorhombic box with periodic boundary conditions was utilized as the solvent and a preset SPC (Simple Point Charge) water model as the solvent. An electrically neutral simulation system was constructed of 0.15 M NaCl was developed. Using the Nose-Hoover thermostat algorithm and the Martyna-Tobias-Klein Barostat method, we stabilized the temperature at 300 K and the pressure at 1 atm throughout the simulation[36–38]. For 100 ns, the simulation was run under TNP conditions, and the trajectory information was acquired using the Berendsen thermostat and barostat methods throughout the remaining 2 ps[38].
In order to calculate binding free energies of ligand-protein complexes, molecular mechanics and generalized Born surface area (MM-GBSA) were operated. A simulation with the VSGB solvation model and the OPLS3e force field was run with a 100-step sample size on frames 0-501 to estimate the Prime MM-GBSA binding free energy. For Prime MM-GBSA's free binding energy (kcal/mol), the additivity principle was employed in which each component's separate energy modules were combined to estimate the binding free energy of the protein and ligand.
The molecular properties of a compound influence pharmacodynamics and its pharmacokinetics. Hence optimize their absorption, distribution, metabolism and excretion (ADME) in the human body as a drug. Lipinski's rule validation demonstrate good oral bioavailability and a balance between the aqueous solubility of a compound and its ability to passively diffuse through different biological barriers.
furthermore, Veber's rule indicates that molecular flexibility and polar surface area (PSA) are crucial factors of oral bioavailability[39]. This molecular flexibility governed by the number of rotating bonds. Reduced molecular flexibility and low PSA are important predictors of good oral bioavailability that molecular flexibility determined by the number of rotating bonds. The reduction in molecular flexibility and the low PSA are important predictors of good oral bioavailability[40].
A CoMFA analysis was generated with static and electrostatic fields and, 31 possible combinations of five descriptors for the CoMSIA analysis, as the dependence on the five CoMSIA descriptors results in a decrease in the significance of the model[41]. We selected the models with values Q² and R² were respectively greater than 0.5 and 0.6 Table 2.
Model |
No validation |
Leave-One-Out (LOO) |
Field distribution (%) |
||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
R2 |
SEE |
F |
R²pred |
Q2 |
NOC |
S |
E |
H |
D |
A |
|||
CoMFA |
E+S |
0.987 |
0.066 |
208.357 |
0.553 |
0.553 |
6 |
44.2 |
55.8 |
- |
- |
- |
|
CoMSIA |
H + A |
0.988 |
0.064 |
223.422 |
0,559 |
0.533 |
5 |
- |
- |
39.6 |
- |
60.4 |
|
S + H + A |
0.988 |
0.065 |
217.030 |
0,794 |
0.563 |
5 |
15.8 |
- |
31.2 |
53 |
|||
H + D + A |
0.982 |
0.079 |
143.318 |
0,945 |
0.521 |
5 |
- |
- |
28.5 |
25.9 |
53.0 |
||
S + E + H + A |
0.990 |
0.06 |
256.781 |
0,511 |
0.511 |
5 |
11.5 |
29.6 |
21.4 |
- |
37.5 |
||
S + E + H + D |
0.981 |
0.081 |
136.993 |
0,994 |
0.502 |
6 |
12.8 |
37.9 |
24 |
25.4 |
- |
||
S + H + D + A |
0.982 |
0.078 |
148.929 |
0,898 |
0.541 |
5 |
11.8 |
- |
23.4 |
23.4 |
41.4 |
||
S + E + H + D + A |
0.986 |
0.070 |
183.314 |
0,855 |
0.517 |
5 |
9.4 |
23.3 |
17.3 |
19 |
31 |
In the case of the CoMFA model, the values for the optimal number of components (ONC), Q², R2, F and SEE were obtained, at 6, 0.553, 0.987, 208.357 and 0.066, respectively. The contributions of steric and electrostatic fields to the CoMFA model were 44.2 and 55.8%, respectively. Among the seven CoMSIA models, the best model is the CoMSIA / S+H+A model with steric, hydrophobic and acceptor hydrogen bond fields, the Q², R², F and SEE values were determined at 0.563, 0.988, 217.030 and 0.065, respectively. The contributions of steric and hydrophobic fields and acceptor hydrogen bond to the CoMFA model were 15.8, 31.2 and 53%, respectively, signaling that the acceptor hydrogen bond field created the higher contribution to the binding affinity.
The CoMFA and CoMSIA models have been validated by the reserved test set. The R²pred values of CoMFA and CoMSIA were calculated to be 0.553 and 0.794, respectively. Which shows the strong predictability of the CoMSIA model compared to the CoMFA model.
The construction of the HQSAR model was adopted by three factors, containing the fragment size, the fragment distinction and hologram length. The statistical parameters of the established HQSAR models are presented in Table 3-4, the parameters R², Q² and SEE were used to display the validity of the models.
Model |
Fragment |
Q² |
SEECV |
R² |
SEE |
Best HL |
NOC |
---|---|---|---|---|---|---|---|
13 |
B+H |
0.660 |
0.323 |
0.796 |
0.251 |
257 |
4 |
14 |
B + Ch |
0.650 |
0.331 |
0.709 |
0.211 |
257 |
4 |
15 |
B + DA |
0.591 |
0.354 |
0.835 |
0.225 |
257 |
4 |
35 |
B + H + Ch |
0.660 |
0.323 |
0.796 |
0.251 |
257 |
4 |
36 |
B + H + DA |
0.591 |
0.354 |
0.835 |
0.225 |
257 |
4 |
37 |
B + Ch + DA |
0.591 |
0.354 |
0.835 |
0.225 |
257 |
4 |
45 |
A + B + H + Ch |
0.494 |
0.418 |
0.925 |
0.161 |
199 |
6 |
48 |
A + C + H + Ch |
0.508 |
0.388 |
0.877 |
0.195 |
353 |
4 |
54 |
B + H + Ch -+DA |
0.591 |
0.354 |
0.835 |
0.225 |
257 |
4 |
57 |
A + B + C + H + Ch |
0.602 |
0.371 |
0.962 |
0.115 |
97 |
6 |
The most excellent hologram model was generated using a histogram length of 257 having four optimal components. The descriptors used to create the model were B (bonds) and H (hydrogen atom), the addition of descriptor Ch (chirality) on this model does not add any advantage in terms of the statistical parameters. The best model generated had a Q² cross-validation coefficient value of 0.66 and a non-validation R² value of 0.796 with a standard error of 0.323. and the fragment size of 4 – 7.
Model |
Fragment |
Q2 |
SEE cv |
R2 |
SEE |
Best LH |
NOC |
---|---|---|---|---|---|---|---|
1 - 4 |
B + H |
0.284 |
0.497 |
0.722 |
0.310 |
353 |
6 |
2 - 5 |
B + H |
0.336 |
0.465 |
0.803 |
0.253 |
199 |
5 |
3 - 6 |
B + H |
0.618 |
0.352 |
0.848 |
0.222 |
257 |
5 |
4 - 7 |
B + H |
0.660 |
0.323 |
0.796 |
0.251 |
257 |
4 |
5 - 8 |
B + H |
0.633 |
0.327 |
0.759 |
0.265 |
199 |
3 |
6 - 9 |
B +H |
0.600 |
0.351 |
0.820 |
0.235 |
199 |
4 |
7 - 10 |
B + H |
0.523 |
0.372 |
0.762 |
0.263 |
307 |
3 |
8 - 11 |
B + H |
0.585 |
0.347 |
0.783 |
0.251 |
97 |
3 |
9 - 12 |
B + H |
0.580 |
0.349 |
0.782 |
0.252 |
257 |
3 |
10 - 13 |
B + H |
0.587 |
0.347 |
0.795 |
0.244 |
97 |
3 |
11 - 14 |
B + H |
0.543 |
0.397 |
0.969 |
0.103 |
257 |
6 |
12 - 15 |
B + H |
0.606 |
0.348 |
0.925 |
0.151 |
267 |
4 |
13 - 16 |
B + H |
0.648 |
0.349 |
0.973 |
0.097 |
199 |
6 |
The predictive capacities of the QSAR models were validated externally using the six molecules reserved for the set of tests [29]. As it was mentioned after that the parameters Q2 and R2, obtained from internal validation, are necessary but insufficient to confirm the stability and the predictive capacity of the models.
The results of the external validation are listed in Table 5. The best QSAR models for the test set gave the R2pred, r²m, r’²m, \(\stackrel{-}{r}\)²m and Δr2m values of 0.555, 0.554, 0.215, 0.385 and -0.337 (CoMFA), 0.794, 0.795, 0.574,, 0.634 and -0.320 (CoMSIA), 0.559, 0.286, 0.381, 0.470 and -0.178 (HQSAR) and the slope regression slope values k and k’ of 1.011 and 0.981 (CoMFA), 1.014 and 0.985 (CoMSIA), and 1.052 and 0.949 (HQSAR) for MDM2–p53 activity inhibitory, respectively. The values of the parameters r2o and r’2o are 0.975 and 0.924 (CoMFA); 0.962 and 0.959 (CoMSIA), and 0.500 and 0.458 (HQSAR) for inhibitory activity, respectively, were used to calculate the relationships (r2-r2o) / r2 and (r2-r’2o) / r2values were obtained from -0,765 and -0.673 (CoMFA), -0.211 and -0.204 (CoMSIA), and 0.106 and 0.182 (HQSAR).
Statistical Parameters |
3D-QSAR |
2D-QSAR |
||
CoMFA |
CoMSIA SHA |
HQSAR |
||
Tropsha and Golbraikh’s criteria |
r0² |
0.975 |
0,962 |
0.500 |
r'0² |
0.924 |
0,956 |
0.458 |
|
R²pred |
0.555 |
0,794 |
0.559 |
|
k |
1,011 |
1,014 |
1,052 |
|
(R² + R²0) /R² |
-0,765 |
-0,211 |
0,106 |
|
k’ |
0,981 |
0,985 |
0,949 |
|
(R² + R’²0) /R² |
-0,673 |
-0,204 |
0,181 |
|
Roy’s criteria |
rm2 |
0.554 |
0,795 |
0.286 |
r'm2 |
0.215 |
0,574 |
0,381 |
|
\(\stackrel{-}{r}²\)m |
0.385 |
0,634 |
0,470 |
|
Δrm² |
-0,337 |
-0,320 |
-0,178 |
|
Δr0² |
0,051 |
0.006 |
0,042 |
Considering the values of the performance parameters provided by the external validation (Table 5), it is obvious that among the three models generated, only the CoMSIA model respected all the criteria of Tropsha and Glorbaikh, and the metrics of Roy. these results validated the predictive power of the CoMSIA / S+H+A model to estimate the biological activities of new compounds.
The correlation plots between the experimental pki and predicted pki and the residuals of the CoMFA and CoMSIA and HQSAR models are displayed in Fig. 3a, 3b and 3c, respectively.
The ability to visualize the results as three-dimensional contour plots is an intriguing feature of CoMFA and CoMSIA modeling. The contour maps were created by multiplying the coefficients and standard deviations associated with each CoMFA or CoMSIA column by a scalar product. The resulting maps depict regions with scaled coefficients greater than or equal to 80% (favored/unfavored). They will give us with sufficient information about the regions that are favorable or unfavorable for the compounds' activity, allowing us to tweak and strengthen the inhibitors previously reported.
The steric and electrostatic fields for CoMFA analysis are illustrated by contour plots on the most active compound N°04 as seen in Fig. 4. In the case the steric field CoMFA, green contours indicate regions where bulky groups improve biological activity, while yellow contour indicate regions where bulky groups decrease activity. The presence of two significant zones of green contours around the chlorine atom of R1 substitution and the CF3 group of R3 substitution indicates that the substituents at these locations have favorable steric interactions. A small green contour between the para and ortho phenyl ring positions of the R4 substitution suggests that the bulky substituents will rise proportionally biological activities. Three yellow polyhedrons concentrated around the meta position of the R3 substitution phenyl ring suggest that less bulky substituents in the position will improve biological activity (Fig. 4a).
The Fig. 4b was the contour map of the electrostatic field, the electronegative and electropositive regions were represented by blue and red contours, respectively. Red and blue contours indicate regions where bulky groups improve biological activity, while yellow contour indicate regions where bulky groups decrease activity. Three red contours around the CF3 group, as well as the NO2 and methyl of the isoxazole group, signal that these positions are appropriate for electronegative groups to considerably increase biological activity (Fig. 4b). Two small blue contour regions between pyrrole and isoxazole and near the core of R3 substituent phenyl indicate positively charged substituents may improve activity.
In the case of the CoMSIA model, the steric, hydrophobic, hydrogen bond acceptor contours are presented in Fig. 5. In the steric field, green contours indicate regions where bulky groups improve biological activity, while yellow contour indicate regions where bulky groups decrease activity. In the map of steric contours of CoMSIA, we noticed two large contours one near the CF3 group of R4 substitution, and the other covers the ortho and meta position of R3 substitution phenyl which indicates that a bulky substitution is required in this region. Furthermore, we could also see two large yellow contours covering the two faces of the isoxazole ring showing that the substitution in this position is unfavorable to the increase in activity Fig. 5a. In the hydrophobic field, the yellow contour indicates the regions where the hydrophobic areas are favorable, while the gray contour indicate the regions where the hydrophile regions are favorable. the hydrophobic grey contours covering the R4 substitution CF3 group, the phenyl para position of the R3 substitution and also NO2 group, suggest that the substitution of hydrophilic groups in these positions could enhance biological activity. while yellow outlines localize near the chlorine atom of group R1, and around oxindole show that the substitution of the hydrophobic groups could enhance the activity Fig. 5b.
The hydrogen bonds acceptor field (A) plays a major role in the bioactivity compared to the steric and hydrophobic fields (percentage contribution 53%). The cyan colored contour and the purple colored contour denote the hydrogen bond acceptor favorable and unfavorable region. As shown in Fig. 5c, a large cyan contour near the NO2 group indicates that the hydrogen bond acceptor group at this position suitable for activity. The unfavorable hydrogen bond acceptor red contour was noticed near the ring of the phenyl ring, indicating that the hydrogen bond acceptor group is not favorable at this position.
On the most active compound N°04, as well as the least active compound N°20, an atomic contribution graph generated by the HQSAR model is shown Fig. 6. The different colors indicate the degree which each atom contributes to the total biological activity. The negative contributions are illustrated by the red, red-orange, and orange patterns, while the positive contributions are displayed by the yellow, green, and blue-green patterns. Grey colors have been used to indicate intermediate contributions. The most common substructure has been indicated in cyan. The appearance of red or orange phenyl hydrogen indicates that it contributed negatively to the activity. The hydrogen atoms of isoxazole's methyl group were whitish, indicating an intermediate contribution.
The results of established 3D-QSAR models were used for the design of new potent inhibitors, as shown in the Table 6. Compound N°04 exhibited the highest activity and was taken as a template to design new compounds. In the phenyl ring at the R3substitution, there was the favorable steric contour between the para and meta position, suggesting that bulky substitutions increased biological activity. As a result, cyclohyptane was introduced at this position, and the Pred01 compounds were designed.
Table 6 Structure, predicted pKi and affinity of newly designed compounds
The four designed compounds, the most potent inhibitor (N°04), as well as the control were docked in the binding site of protein, for the purpose of detecting the key amino acids in the biological activity studied, as well as to acquire in-depth insight into crucial interaction modes. We selected the conformation with the highest affinity score value among nine output conformations. According to the docking result, the newly designed compound Pred01 had the highest binding affinity -9.4 kcal / mol, while compound N°04 and Nutlin showed an affinity of -8.4 kcal / mol and -9.2 respectively Table 6. The amino acid His92 is in a suitable position to be involved in the formation of hydrogen bond of compound NO2 (2.48 Å) and compound Pred01 (2.83Å) with the oxygen atom of NO2 group and the nitrogen atom of isoxazole respectively. These observations are precisely in agreement with the contour maps of the H-binding acceptor in the CoMSIA odels. In addition, the compound Pred03 form a hydrogen bond with the same His92 residue on distance 2.56 Å. These results indicating this amino acid plays an essential role in controlling activity (Fig. 7). Thus, these residues may be critical indications for rational design to find out more potent MDM2–p53 inhibitors.
The docking methodology was validated by the redocking technique by redecking the co-crystallized ligand, then calculating the RMSD between the reference co-crystalline and the re-docked co-crystallized as shown in Fig. 8. The result of the docking showed that the docking of the co-crystallized ligand shows an RMSD value equivalent to 0.87 Å, indicates that the methodology is accurate in predicting the binding affinity for ligands.
The lipinski rules are widely used to assess drug-like compounds. According to these rules, compounds that do not meet at least two of the following criteria are very likely to have problems with absorption or oral permeability[42]. The Veber and lipinski properties of the newly designed compounds are collected in Table 7.
Comp. |
Lipinski rule |
Veber rules |
Synthetic accessibility |
|||||
---|---|---|---|---|---|---|---|---|
MW |
NBA |
NBD |
Log P |
Molar Refractivity |
Nrot Bond |
PSA |
||
Pred01 |
532.90 |
9 |
2 |
6.01 |
136.10 |
4 |
112.98 |
5.33 |
Pred 02 |
534.49 |
9 |
2 |
6.53 |
143.90 |
5 |
112.98 |
5.44 |
Pred 03 |
496.94 |
7 |
2 |
4.47 |
138.83 |
5 |
122.21 |
5.44 |
Pred 04 |
468.85 |
8 |
2 |
2.82 |
124.96 |
5 |
139.28 |
5.01 |
N°04 |
492.84 |
9 |
2 |
5.52 |
123.47 |
4 |
112.98 |
5.04 |
According to Table 7, the most soluble compound in aquatic and lipid solutions is Pred 01, followed by Pred03, N°04, Pred01, and Pred02. Giving to Lipinski's rules, these compounds should have no difficulty with oral bioavailability. We note that all compounds have PSA and Nrot Bond values of less than 140A2 and 10 respectively, indicating that they are well predicted for oral bioavailability and membrane transport. We can convince these compounds to comply with Veber's rule. As shown in Table 7, the synthetic accessibility values of the new compounds are located within a range closest to the representative molecule of the dataset N°04 and can be considered encouraging.
In this part, we have selected three complex N°04/receptor, Nutlin-3/receptor and Pred 01/receptor. For a study molecular dynamics simulations were performed in 100 ns to assess conformational stability very much closer to physiological environmental conditions .The RMSD trajectories of free protein and ligand-protein complexes were used to assess the stability of the selected systems as were illustrated in Fig. 10.
RMSD values of free protein alpha carbons showed changes during the first simulation of 70 ns, then became stable after 70 ns for all three systems. For the ligand N°04, Nutlin-3 and Pred 01, there were fluctuations at the start, then all the complexes gradually reached equilibrium at 56ns, 58ns and 75ns of simulation respectively. During the last 35 ns of simulation, the RMSD values of compound Pred 01, Nutlin-3 and compound N°04 varied around 1.6Å, 2.4Å and 4.2Å respectively, indicating that compound Pred 01 could form a more stable complex with the target protein binding site Fig. 9.
The evolution of RMSF for the three complexes, indicating the local changes in the protein chains, generally the three complexes fluctuate weakly such that does not exceed 3Å, but the complex relating to compound Pred 01 fluctuates more weakly than the others and does not exceed 2Å Fig. 10.
The interactions of protein/ligand and the graphs of contacts were monitored throughout the simulation are also illustrated in Fig. 11. During the MD simulation equilibrium steps the protein-ligand interactions were evaluated to identify the crucial residues involved in the interaction for the three ligands. For compound N°04, the most stable interactions were hydrophobic conctates with His92, Leu50, Ile95 and Ile 57 accounting for 30, 28, 26, 17% of the whole Simulation Trajectory, respectively. Other hydrogen bond interactions formed with His92 and Gln55 were also observed, His92 formed a hydrogen bond with the nitrogen dioxide group, with percentages of 12 and 5% respectively, indicating the hydrophobic interaction formed with these residues played an important role in the inhibitory activity interacting with the protein.
In addition, the most significant water bridges were observed with Gln55, Tyr51 and Gly54 residues accounting for 23, 21 and 13% respectively of all MD simulation time, indicating the important role of these bindings in the receptor-ligand interaction.
In cas of the Nutlin-3 exhibited hydrogen bonds with Gln55 and Tyr63 with percentages of 30% and 3% respectively. In general, the most required interactions are hydrophobic types formed with Ile95, Gln55, Ile57 Val89 and Leu 50 for 57, 40, 34, 29 and 20%n respectively. Nevertheless as for the compound Pred01 established hydrogen bond, hydrophobic and water bridges interactions with Val89 and His92 and His 69 for over 30% and 45% 11% of the simulation time, respectively. Moreover the residue Tyr63 created the most important interaction for 58%.
Figure 11 depicts a temporal representation of interactions and contacts. The top right corner panel shows the total number of interactions established by the protein with the ligand during MD, whereas the top left panel illustrates the residues that interact with the ligand in each route frame. Tyr63, Val89, and His92 are the residues that establish more specific interaction with the ligand Pred01.
The binding free energy was estimated by Coulomb energy (Coulomb), covalent bond (Covalent), hydrogen bond (Hbond), lipophilic bond (Lipo), π-π Packing interaction, the solvent the generalized binding and the binding of Vander Waals (VDW) the values of the contribution of each element on the free energy of total binding between the compounds N°04, Nutlin-3, Pred 01 and the binding pocket of the target protein were has been shown in Table 8.
Comp. |
MMGBSA (kcal/mol) |
|||||||
---|---|---|---|---|---|---|---|---|
ΔGBind |
ΔGCoulomb |
ΔGCovalent |
ΔGLipo |
ΔGSolv GB |
ΔGvdW |
ΔG H bond |
ΔG Bind Packing |
|
Compound N°04 |
-43.461 |
2.831 |
0.983 |
-14.568 |
3.082 |
-34.878 |
-0.127 |
-0.782 |
Nutlin-3 |
-71.384 |
-9.172 |
3.182 |
-27.396 |
19.002 |
-55.852 |
-0.289 |
-0.859 |
Pred 01 |
-60.700 |
0.131 |
2.460 |
-20.501 |
5.471 |
-46.089 |
-0.197 |
-1.970 |
The value of average ΔGBind was -43.46 kcal / mol for the complex relating to N°04, -71.38 kcal / mol for the complex relating to Nutlin-3 and -60.70 kcal / mol for the Pred 01 complex. In general, the newly designed Pred 01 compound and the control compound showed the lowest total binding free energy. Among all interactions, the contributions of Vander Waals energy and lipofil binding were greater than those of other elements with values of -55.85 kcal / mol and -27.39 kcal / mol for Nutlin-3 and -46.09 kcal / mol and -20,50 kcal / mol for Pred01 respectively. The Hydrogen bond interaction values for Pred 01 and the most active compound N°04 were -0.20 and -0.13 kcal / mol, respectively, respectively. Which clearly indicated that the newly designed compound formed more stable hydrogen bonds with binding site residues than Compound N°04.
A database of 29 compounds was used to modulate inhibitory activity against GPX4/MDM2 - p53. The established QSAR model was determined by examining the structure-activity relationship of nitroisoxazole containing spiro [pyrrolidine-oxindole]. The three best models obtained CoMSIA / S+H+A, CoMFA / E+S and HQSAR / B + H with statistically significant indicators giving R² = 0.988; Q² = 0.560, R² = 0.987; Q² = 0.553 and R² = 0.796; Q² = 0.660 respectively. Based on, the CoMSIA model passed all external validation tests and met the criteria of Tropsha and Globraikh. It was selected to design four new candidates with more potent inhibitory activity than the most active molecule in data set N°04.
In addition, the newly predicted pred01 compounds, N°04 and Nutlin-3 as the control compound were selected for the molecular docking study and the 100 ns MD simulations to analyze the mode of ligand interaction in the binding site of the receptor and assess the stability of the docked complexes respectively. To obtain more information, the binding free energy was examined using the MM-GBSA method. The value of the free energy confirms the stability of the Pred 01 ligand showing the lowest energy. All of these interesting results can be invested in and used to design a nitroisoxazole derivative drug containing spiro [pyrrolidine-oxindole] as more potent inhibitors of GPX4 / MDM2 - p53.
We are grateful to the “Association Marocaine des Chimistes Théoriciens” (AMCT) and “Moroccan Centre of Scientific and Technique research” (CNRST) for their pertinent help concerning the programs.
Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript
Competing Interests The authors have no relevant financial or non-financial interests to disclose.
The authors declare that they have no conflict of interest.
Ethical approval This chapter does not contain any studies with human participants or animals performed by any of the authors.
Author information
Affiliations
School of Sciences, Moulay Ismail University, Meknes, Morocco
Kamal Tabti, Larbi Elmchichi, Abdelouahid Sbai, Hamid Maghat & Tahar Lakhlifi
High School of Technology, Moulay Ismail University, Meknes, Morocco
Mohammed Bouachrine
Contributions
Kamal Tabti: Data curation, Writing – original draft; Larbi Elmchichi: Visualization, Investigation; Abdelouahid Sbai: Conceptualization, Methodology, Software; Hamid Maghat: Supervision. Mohammed Bouachrine Software, Validation; Tahar Lakhlifi: Writing – review & editing. All authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Corresponding author
Abdelouahid Sbai