3D-QSAR statistical results
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.
Table 2
Statistical parameters of CoMSIA models with different combinations of molecular fields.
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.
HQSAR statistical results
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.
Table 3
HQSAR analysis of various fragment distinctions use up fragment size (4 – 7)
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.
Table 4
Statistical parameters of HQSAR model using the fragments B+H with different fragment sizes.
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
|
Validation of QSAR models
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).
Table 5
Golbraikh and Tropsha test results
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.
Visualization contour carts CoMFA and CoMSIA
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.
Hologram QSAR Atomic Contribution Plot Interpretation
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.
Designed molecular docking results and analysis
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.
Lipinski rules and veber parametrs
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.
Table 7
lipinski and veber rules for newly designed compounds
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.
Molecular dynamics analysis
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.
Binding free energy calculation using prime molecularmechanics/generalized born surface area (MM-GBSA) approach
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.
Table 8
Post simulation of binding free energy for protein-ligand complexes estimated using MM-GBSA analysais
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.