Computational assessment of the reactivity and pharmaceutical potential of novel triazole derivatives: An approach combining DFT calculations, molecular dynamics simulations, and molecular docking

2.1 Data set and training/test separation

A database of Twenty-two mansonone E-linked triazole derivatives and their biological activity against topoisomerase II (expressed as IC₅₀ values) was brought from the literature(Huang et al., 2013). IC₅₀ values of dataset span the entire activity range from 1.49 µM to 98 µM were converted to pIC₅₀ (pIC₅₀ = -log (IC₅₀ [10^-6])) and were used as dependent variables in the QSAR analysis(El Masaoudy et al., 2023). The compound dataset was randomly divided into a training set (77 %, 17 compounds) which was used to build the QSAR models and a test set (23 %, 5 compounds) which was used to assess the predictive ability of the models. The structures and their biological activity are shown in Table 1.

Table 1 Molecular structure of mansonone E-linked triazole derivatives dataset and their corresponding biological activity values pIC₅₀.

Compounds			pIC50 (Exp)	HQSAR		CoMFA		CoMSIA		Affinity Kcal/mol
				pIC50 (calc)	Residus	pIC50 (calc)	Residus	pIC50 (calc)	Residus
	R1	R2
1	4-F-Ph	Cl	4.973	4.943	0.03	4.973	0.000	4.968	0.005	−10.4
2	4-Cl-Ph	Cl	4.914	4.943	−0.029	4.914	0.067	4.82	0.094	−10.5
3	Ph	Cl	4.735	4.573	0.162	4.735	−0.109	4.7	0.035	−9.9
4	4-CH3-Ph	Cl	4.703	4.85	−0.147	4.703	−0.013	4.87	−0.167	−10.6
5	4-CH2CH3-Ph	Cl	5.538	5.411	0.127	5.538	0.067	5.457	0.081	−10.2
6	−4-(CH2)2CH3-Ph	Cl	5.623	5.65	−0.027	5.623	0.011	5.589	0.034	−10.2
7	−4-(CH2)3CH3-Ph	Cl	5.15	5.66	−0.51	5.593	0.000	5.593	0.000	−9.9
8*	−4-(CH2)4CH3-Ph	Cl	5.593	5.694	−0.101	5.451	0.000	6.145	−0.694	−9.6
9*	4-OCH3-Ph	Cl	5.451	5.019	0.432	5.150	−0.137	5.289	−0.139	−9.6
10		Cl	4.301	4.314	−0.013	4.301	−0.016	4.374	−0.073	−10.1
11		Cl	5.569	5.617	−0.048	5.569	−0.018	5.625	−0.056	−10.7
12	−4-F-Ph	Br	4.91	4.85	0.06	4.910	−0.038	4.912	−0.002	−10.6
13	−4-Cl-Ph	Br	4.813	4.85	−0.037	4.813	−0.033	4.833	−0.020	−10.6
14	Ph	Br	4.533	4.573	−0.04	4.533	0.078	4.463	0.070	−10.1
15*	4-CH3-Ph	Br	4.533	4.85	−0.317	4.988	−0.279	5.037	−0.049	−10.1
16	4-CH2CH3-Ph	Br	5.306	5.411	−0.105	5.306	−0.007	5.226	0.080	−10.4
17	−4-(CH2)2CH3-Ph	Br	5.827	5.65	0.177	5.827	−0.064	5.896	−0.069	−9.6
18	−4-(CH2)3CH3-Ph	Br	5.593	5.66	−0.067	5.593	0.014	5.615	−0.022	−9.9
19*	−4-(CH2)4CH3-Ph	Br	5.630	5.694	−0.064	5.631	0.003	5.939	−0.308	−9.6
20	4-OCH3-Ph	Br	4.987	5.019	−0.032	4.987	0.002	5.026	−0.039	−9.8
21*		Br	4.643	4.314	0.329	4.644	0.045	4.432	0.212	−9.7
22		Br	5.676	5.617	0.059	5.676	0.000	5.629	0.047	−10.5

2.2 Structural optimisation and alignment of the dataset

All chemical structures of molecules in the dataset were constructed and drawn in SYBYL X2.0 software (Tripos, Inc., St. Louis, MO, USA). The Tripos force field was included to minimize the structural energy of molecules while accounting for Gasteiger-Huckel atomic partial charges. Then, a convergence threshold of the Powell gradient algorithm was set at 0.05 kcal/mol. To ensure the conformational stability of the molecules, the number of iterations has been increased to 10,000 iterations. On another side, the molecular alignment of the dataset is a crucial phase for QSAR modeling, and the best alignment is reflected on the statistical results of the generated 3D QSAR models (Lu et al., 2010; Pandey and Saxena, 2006). In this present work, the distill rigid method has been adapted thanks to its better results for the generation of the 3D-QSAR model(El Mchichi et al., 2022). The most active molecule N17 was taken as a template, the results are shown in Fig. 2.

Close

Fig. 2

(a) Molecular alignment by rigid distill of the dataset, (A) most active N17 as template, (B) Common core of the data set compounds and (c) dataset alignment employed by distill rigid for 3D-QSAR modeling.

Export to PPT

2.3 Regression analysis and validation of QSAR models

The PLS regression, executed using SYBYLY-X2.0 software, was employed to quantify the linear relationship between biological activity (pIC50) and molecular descriptors. To determine the optimal number of components and the cross-validation coefficient Q², as well as to statistically evaluate the established models (Vong et al., 1988), cross-validation was performed using the Leave-One-Out (LOO) method. However, to assess the statistical confidence of the generated models, a non-validation was carried out by calculating the coefficient of determination R2 and the standard error of estimation (SEE). It is important to note that a reliable evaluation of 3D-QSAR and HQSAR models through internal cross-validation alone is insufficient. Therefore, external validation is required, such as employing the Golbraikh-Tropsha criteria and Roy's metrics(Tabti et al., 2023d). Certain conditions need to be satisfied for the successful execution of this external validation, including(Tabti et al., 2023a):

Q² > 0.5, R² > 0.6, (R² – R²₀)/R² < 1, 0.85 ≤ K ≤ 1.15 (Golbraikh-Tropsha criteria);

r² _m < 0.2, r² _m > 0.5 (Roy's metrics).

On another side, The evaluation of the robustness of a QSAR model is based on carrying out a Y randomization test(Khan et al., 2019). In this test, a new model is generated by randomly shuffling the vector of the dependent variable. This procedure is repeated several times. The resulting QSAR model is only considered valid if the new QSPR models generated by randomization have low R² and Q² values.

2.4 The applicability domain (AD)

To choose an appropriate QSAR model for new compounds, it is essential to delineate its applicability domain of application AD((Tropsha et al., 2003). Model predictions should only be considered reliable for compounds included in this AD. Two simple approaches to define this scope are range extrapolation, based on the calculation of Euclidean distances and leverage (hi) for each chemical compound. Model predictions should only be considered reliable for compounds included in this scope. Among the simple approaches to define this field of application is the calculation of Euclidean distances and leverage (hi) for each chemical compound according to the following formula: $$h_i=x_i{\left(X^{T}X\right)}^{-1}{x}_i^T$$ The vector xi represents the k model parameters of the compound in question, while the matrix X is an n × k matrix containing the model parameters for each of the n training compounds. A leverage value exceeding 3(k + 1)/n is considered significant. This indicates that the model prediction, which can potentially compromise its reliability.

2.5 Limitations of the models developed

The models have notable limitations that require full recognition for informed interpretation of their results. First, these models are confronted with small data sets. Furthermore, the compounds examined only show moderate activity, which could potentially influence predictions in favor of this type of molecules. Another constraint lies in the absence of double-blind (nested) cross-validation, which could lead to overestimation of model performance(Baumann and Baumann, 2014). Additionally, using simple cross-validation, particularly Leave-One-Out (LOO), carries a risk of overfitting the training data, thereby jeopardizing the models' ability to generalize reliably(Chtita et al., 2021). To overcome these limitations, the combination of HQSAR, CoMFA and CoMSIA models presents significant advantages in terms of molecular modeling. This global approach exploits the complementarity of the information provided by each model, thus improving the reliability of the predictions. Additionally, it strengthens cross-validation, increasing confidence in the results. This synergy facilitates the precise identification of molecular descriptors relevant for drug design and improves the consideration of chemical diversity, thereby broadening the activity prediction capacity.

2.6 Molecular docking studies and re-docking validation

A molecular docking analysis using the program Autodock Tools vina 1.5.7 was carried out to evaluate the interaction and mode of binding between our studied compounds and the enzyme Topo IIα(Hajji et al., 2021). The 3D X-ray crystal structure of human DNA topoisomerase II α (PDB ID: 1ZXM, resolution = 1.87 Å) was downloaded from the Databank (https://www.rcsb.org). 1ZXM structure is a good target in the ATP catalytic pocket of topoisomerase IIα (Oyedele et al., 2020). The protein structure was prepared by the removal of 280 water molecules and the co-crystallized ligand (phosphoaminophosphonic acidadenylate ester (ANP)). Then the calculated Gasteiger charges and the polar hydrogen atoms were added and the non-polar hydrogens were merged using AutoDock Tools(Tabti et al., 2022c). The protein ligand interactions were framed in a grid defined over 2,0 Å in all three directions using the AutoGrid algorithm(Guo et al., 2005). In addition, the size of the grid space has been taken as the default of 0.375 Å while the central grid box is approximately (35.286 Å, −5.030 Å and 38.362 Å). The choice of best conformation was according to their good affinity towards the protein(Tabti et al., 2023c).

Contact Atlas (http://pca.mbgroup.bio/) as indicated Fig. 3 identifies the amino acids that form strong contacts, the amino acids found (first shell) are Glu-87, Asn-91, Ala-92, Asn-95, Arg-98, Asn-120, Ile-125, Ile-141, Phe-142, Ser-148, Ser-149, Asn-150, Gly-161, Arg-162, Asn-163, Gly-164, Tyr-165, Gly-166, Ala-167, Lys-168 Lys-215 Gln-376 and Lys-378 largely contributed to the affinity towards the selected protein. This vital information allows us to assess the contact patterns of protein ligands.

Close

Fig. 3

The human DNA topoisomerase IIα (PDB 1ZXM) co-crystal structure asteroid plots, The inner shell residues were those that had direct contact with the ligand molecule, whereas the outer shell residues were those that had indirect interactions with the ligand. The number of contacts residues formed with the ligand was represented by the size of the circular nodes.

Export to PPT

2.7 Dynamic molecular analysis and post dynamic binding free energy analysis

The GROMACS 2019 software package with the Charmm force field was used to run a 100 ns MD simulation to validate the stability of the best complex molecular docking results under physiological conditions(Tabti et al., 2022d). The molecular topology file and force field parameters for the ligand were generated with the CHARMM general force field (CGenFF). The system of docked complexes were immersed in a cubic box of TIP3 water molecules with a margin distance of 10 Å. System neutralization was achieved by adding Na + counterions. The entire system was taken through 50,000 energy minimization steps using the steepest descent integrator to release the internal strain energies of the entire system. Then, NVT and NPT balancing were performed for 0.1 ns at 300 K and 1 ns at 1 bar pressure, respectively. System temperature and pressure were controlled by the Berendsen thermostat and the Parrinello/Rahman barostat respectively(Berendsen et al., 1984) (Parrinello and Rahman, 1981). The balanced system was used to perform 100 ns MD simulations at 300 K and 1 bar pressure. All bond lengths were constrained using the LINCS algorithm (Hess, 2008). The Particle mesh Ewald (PME) method was used to handle long-range Coulomb interactions(Essmann et al., 1995). The estimation of MM-PBSA binding free energies was performed by the g_mmpbsa tool on all MD simulation trajectories (Issar et al., 2015; Kumari et al., 2014). The contribution of each residue to the total binding free energy was also estimated using the MMPBSA.py tool(Miller III et al., 2012).

2.8 Pharmacokinetics and Toxicity studies

The study of ADME pharmacokinetics (absorption, distribution, metabolism, excretion) and toxicity is a key factor in the development of drugs(Van De Waterbeemd and Gifford, 2003). Early evaluation of ADME and Toxicity can significantly improve drug success rate, save drug development time costs, reduce drug side effect and toxicity, and further guide rational clinical use drugs(Nel, 1987). In this part of research, SwissADME and pkCSM online servers were used to predict the ADME and Toxicity parameters, drug-likeness and difficulty of synthesis of new compounds designed by the best established model, in order to confirm their possibility of become a drug candidate (Daina et al., 2017; Pires et al., 2015).

2.9 Computational reactivity details

Three different approaches were utilized to investigate the regioselectivity of 1,3DC between Mansonone's azide and 3,3,3-Trifluoroprop-1-yne (Fig. 17). These methods included the Domingo approach, the Gazquez and Mendez approach based on the HSAB principle, and activation energy calculations(Domingo et al., 2013; Mendez and Gazquez, 1994). By employing these approaches, we were able to predict the favored stereoisomer with success. The computations were conducted by the DFT method at the B3LYP/6-31G* level of theory with Gaussian G09 (Becke, 1988; Lee et al., 1988). Additionally, frequency calculations were performed to validate the stationary points and ensure that minima and transition states have zero and one imaginary frequency (Nacereddine et al., 2010). The global reactivity indices such as electronic chemical potential μ, chemical hardness η, global electrophilicity ω, and global nucleophilicity N were estimated using equations suggested by Parr (Parr et al., 1999; Parr and Yang, 1989) and Domingo (Chamorro et al., 2013; Domingo et al., 2013). Specifically, the one-electron energies of the frontier molecular orbitals were employed to evaluate the electronic chemical potentials and chemical hardnesses of the reactants examined in this study.

After determining the values of μ and η, ω was calculated using the following formula:

Nucleophilicity index N has been recently introduced on the basis of the HOMO energies(Domingo and Pérez, 2013):

Because of its low HOMO energy, tetracyanoethylene (TCE) is frequently used as a reference molecule in a wide range of molecules (Jaramillo et al., 2008; Kohn and Sham, 1965).

The chemical softness S was introduced as the inverse of the chemical hardness (Domingo et al., 2002):

In case an amount equivalent to one electron is transferred, the nucleophile becomes a radical cation, while the electrophile becomes a radical anion. Interestingly, analysis of the atomic spin density (ASD) at the radical cation and the radical anion gives a picture of the distribution of the electron density in the electrophile and the nucleophile when they approach each other along the reaction progress.

Based on these observations, in 2014, Domingo proposed the Parr functions P(r) (Chamorro et al., 2013; Domingo et al., 2013), which are given by the following equations:

for electrophilic attack.

for nucleophilic attacks.

Where ρ_s^rc (r) is the ASD at the r atom of the radical cation of a considered molecule and ρ_s^ra (r) is the ASD at the r atom of the radical anion. Each ASD gathered at different atoms of the radical cation and the radical anion of a molecule provides the local nucleophilic P_k^- and electrophilic P_k⁺ parr functions of the neutral molecule.

With these electrophilic and nucleophilic Parr functions at hand, we can redefine the local electrophilicity ω_k, and the local nucleophilicity N_k indices as follows:

The meaning of electron flow in the reaction can be estimated by calculating the dual descriptors γ₁ and γ₂. The closer the y value gets to zero, the greater the preference for the interaction. The equations to calculate γ₁ and γ₂ are as follows:

3.1 CoMFA and CoMFA statistical analysis

The relationship between steric and electrostatic fields and the biological activity of the data set was established using the CoMFA method. The statistical results of the CoMFA model are shown in Table 2. The CoMFA model had a Q² cross-validation of 0.702 for five Optimal Components. PLS analysis gave the conventional R² coefficient of determination value of 0.988, and an SEE standard error estimate value of 0.063. The steric field contributed in the CoMFA model was 67.6 % while the electrostatic field was 32.4 %.

Bold line indicates best significant optimized model.

For the CoMSIA modeling study, the possible models with different possible combinations of fields based on the five different field descriptors were analyzed.

The most reliable models, which met the criterion of (Q² > 0.5), along with their respective parameters and contribution fractions of different fields, are summarized in Table 2. The CoMSIA/SEA model exhibited the highest Q² value of 0.82, utilizing five optimal components. It also achieved a higher R² value of 0.979 and a lower SEE value of 0.086. The contribution of the electrostatic, steric, and acceptor H-bond fields in the CoMSIA model were 50 %, 25.2 %, and 24.7 %, respectively. Additionally, the CoMFA model demonstrated significant static results with a dominant steric field contribution of 67.6 %.

These results strongly support the crucial role of steric fields in the activity of triazole derivatives linked to mansonone E. Overall, the constructed models yielded outcomes within an appropriate range, affirming the robustness and reliability of the 3D-QSAR models.

3.1.1

3.1.1 The domain of applicability

The domain of applicability is defined as the chemical space encompassing the compounds in the model's training set. In this study, the analysis of the area of applicability is based on the “Leverage” method, which is based on the variation of the standardized residuals of the dependent variable in relation to the leverage effect. All observations have standardized residuals in the interval [-2; 2]. It is notable that Fig. 4 does not include any outliers, since the “Leverages” obtained remain below the threshold of 0.52 for the CoMFA model and below the threshold of 0.51 for the CoMSIA model. However, in the case of the CoMFA model, compound 22 exceeds the limit of 2 in terms of standardized residues. Overall, both models tested indicate the absence of outliers.

Close

Fig. 4

The Williams plots of the model CoMFA and CoMSIA.

Export to PPT

3.2 HQSAR statistical results

To generate HQSAR models, the different combinations of fragment distinctions such as atomic number (A), bond type (B), atomic connection (C), hydrogen bond donor–acceptor (DA), hydrogen (H) and chirality (Ch) with a fragment size of 4 to 7 atoms were carried out, then a size change was carried out from 1 to 4 atoms up to 12 to 15 atoms the results obtained are listed in Fig. 5.The top 24 statistically, the best model was model M32 which suggests that parameters such as bond type (B) and hydrogen bond donor–acceptor (DA) were important for generating a model.

Close

Fig. 5

Statistical parameters comparison (a) for the different HQSAR model distinctions and (b) for different size fragment of best Model (M32).

Export to PPT

For further model refinement (M32), the selected model fragment size was explored varying in size from 1 to 4 atoms up to 12–15 atoms. A significant difference was marked in the statistical parameters of the models, as shown in Fig. 5(B). The best model was with a fragment size of 5 atom minimum to 8 atom maximum where Q² and R² were found to be 0.919 and 0.962, respectively, with five componentes optimales and the hologram length was 97. These results verify the impact of fragment size on the robustness of the model.

3.3 QSAR models validation

The validation of the model was continued by carrying out a randomization test Y. We carried out ten random mixtures of the independent vector Y, leading to low values of R² and Q². This observation suggests that the results obtained from the origin model are not influenced by structural dependence or random correlation in the dataset. According to Table 3, the average values of R² were 0.0965, 0.1453 and 0.0649, while those of Q² were 0.0666, 0.0687 and 0.0329 for the ten randomized models CoMFA, CoMSIA and HQSAR, respectively.

The predictive reliability of the established models was assessed by an external validation protocol using the reserved test set. The statistical results of the external validation according to the Tropsha and Globraikh criteria and the Roy metrics have been collected in Table 4.

The R2test values ​​obtained by the CoMFA, CoMSIA and HQSAR models were 0.895, 0.765 and 0.970 respectively. The optimal models generated for the test set gave values ​​of R2m, R'2m, ̅R2m, ΔR2m and ΔR20 of 0.766, 0.607, 0.687, −0.159 and −0.083 (CoMFA); of 0.481, 0.431, 0.456, −0.049 and −0.052 (CoMSIA), 0.818, 0.865, 0.841 and 0.047 and 0.013, (HQSAR) for antiproliferative activity, respectively. It can be emphasized that the values ​​of R'm2, Rm2, and ̅R2m, do not satisfy the criterion greater than 0.5 for the CoMSIA model, while all of them are respected by the CoMFA and HQSAR models. Furthermore, the values ​​of the R2o and R'2o parameters were 0.915 and 0.998; 0.903 and 0.955, and 0.994 and 0.982 for the CoMFA, CoMSIA and HQSAR models, respectively. These parameters are close to 1 for the three generated models. The R2o and R'2o parameters were used to calculate the relations (R2/R2o)/R2 and (R2/R'2o)/R2, the values ​​were obtained from −0.023 and −0.115, −0.181 and −0.248, and −0.025 and −0.012 for the CoMFA, CoMSIA, HQSAR models, respectively.

Generally, we can confirm the credible advantage of the CoMFA and HQSAR model for the external validation phase, while the CoMSIA model showed weaknesses in some external validation static parameters. These results showed that the QSAR models are considered predictive.

The correlation and the residual between the calculated and experimental pIC₅₀ activity values for the CoMFA, CoMSIA and HQSAR models are presented in Fig. 6, respectively.

Close

Fig. 6

Plot of correlation and residual of experimental pIC₅₀ values vs those calculated by (a) HQSAR, (b) CoMFA, and (c) CoMSIA models.

Export to PPT

3.4 Visualization contour carts CoMFA and CoMSIA

In this part, the information derived by the CoMFA and CoMSIA models have been visualized in the form of contour maps using the StDev*Coefficient field type. The most active compound N17 was used as reference structure as shown in Fig. 7 and Fig. 8.

Close

Fig. 7

Contour map of CoMFA analysis (a) Steric, (b) electrostatic with a grid spacing of 2.0 Å in combination with compound N17.

Export to PPT

Close

Fig. 8

Contour map of CoMSIA analysis Steric, electrostatic and hydrophobic with a grid spacing of 2.0 Å in combination with compound N17.

Export to PPT

Contribution levels for Advantaged and disadvantaged regions default to 80 % and 20 %, respectively. In the CoMFA steric contour map, the favorable steric field is represented by green (80 % contribution) and yellow (20 % contribution) represents the unfavorable steric field, as shown in Fig. 7a, areas of yellow color surrounding the phenyl ring of compound N17 proves to be sterically unfavorable. Two areas of green color near the propane ring of the phenyl ring suggest that the addition of bulky groups could improve the activity. We can notice that the absence of yellow and green cards near the mansonone group shows that this group has no steric impact on the biological activity studied. In the electrostatic contour map, the favorable electron-donating group is denoted by blue (80 % contribution) and red represents the favorable electron-withdrawing groups for the activity (20 % contribution). As shown in Fig. 7b, Five blue areas covers around oxygen atom of mansonone, bromine atom, core nitrogen atom of triazole and near propene group bonded to triazole, these areas show that donor groups electrons increase activity. However there is a two red colored areas near the ortho and para position of the phenyl ring which indicates the presence of an electron-withdrawing group at this position could enhance activity.

Contour maps for the CoMSIA model, as shown in Fig. 8a-8c, represent how steric, electrostatic, and hydrogen bond acceptor interactions affect the inhibitory activity of the compounds under study. Among these interactions, the steric interaction has the major contribution to the activity (50 %). The CoMSIA steric contour map (is shown in Fig. 8a) a single green area around the phenyl ring and another green near propane which reflects the addition of larger and smaller groups in the position of phenyl and propane, respectively, are of interest. The CoMSIA electrostatic contour map (Fig. 8b) is less extensive than the corresponding CoMFA one where it shows single electron poor in this area. A small red area near the propane group substituting the triazole nucleus shows the importance of adding electron-rich groups in this region. In the hydrogen bond acceptor contour map (Fig. 8c), a large red contour around the naphthoquinone rings, and two other red p in the vicinity of triazole nucleus indicating that the insertion of hydrogen bond acceptor group in this area could increase activity. A small magenta outline covering the ortho position of the phenyl ring shows negative impact of hydrogen bond acceptor group on biological activity.

3.5 Hologram QSAR atomic contribution plot interpretation

The HQSAR model contribution map shows the degree of contribution of each atom, the more active N17 and less active N10 compounds were selected as reference structures (Fig. 9). The HQSAR hologram helps us to better understand and identify the impact of each atom fragment on the inhibitory activity of antiproliferative. As shown in Fig. 9, the common backbone appears in cyan color, and the effect of each atom on biological activity was indicated by a color. The green or yellow color on the contribution map reflects positive impacts, while the orange and red presets represent negative influences. Fragments of atoms colored white indicate intermediate contributions in biological activity.

Close

Fig. 9

Hologram quantitative structure activity relationship model contribution map of the most compound N17 and the less compound N10.

Export to PPT

It can be seen in Fig. 9 Both structures contain a mansonone E linked to triazole moiety as a common core. In addition, some of the hydrogen atoms on the naphthoquinones cycle and the pyranic cycle appear yellow and green, which have a significant contribution to enhance inhibitor activity. Moreover, the bromine atom appears in green, which confirms the important role of this atom in improving the activity, this result confirms the results obtained by Hong and his collaborators. The atoms of the studied compounds do not appear in red indicating that there is no negative contribution. In addition, the atoms and the substituent appear in white to indicate a neutral contribution in the substituent R1.

3.6 Designed molecular

Analysis of the contour map of the 3D QSAR models and the HQSAR hologram contribution map has derived interesting information on the structural features that affect and govern the antiproliferative activity of triazole mansonone E derivative. The results obtained from the QSAR study carried out have been schematized in Fig. 10.

Close

Fig. 10

Summary of structure–activity relationship derived from QSAR study.

Export to PPT

Investing the information generated by the models, and according to the results of external validation, the CoMFA model was selected to design seven new compounds taking into account compound N17 as a reference. Predicted pIC₅₀ values were calculated using three models CoMFA, CoMSIA and HQSAR (Table 5). The predicted activity of all designed molecules was high (pIC₅₀ > 5.624) and found to be quite similar based on CoMFA and CoMFA/ESA models and slightly different from CoMSIA model. For this, it was taken as an optimal model to predict the pIC₅₀ activity of the new molecules, for example the activity predicted by CoMSIA for the T1 molecule is pIC₅₀ = 8.837 were found to higher than the activity predicted by CoMSIA of compound N17 the highest active compounds from the database pIC₅₀ = 8.103.

Table 5 Structure, predicted pIC_50, and synthetic accessibility of newly designed compounds.

Compds	Structure	Activity predicted pIC50			Synthetic accessibility	Bioavailability Score
		CoMSIA	CoMFA	HQSAR
T1		4.570	6.631	4.326	3.81	0.85
T2		6.871	6.271	5.624	4.60	0.56
T3		6.436	6.179	6.465	4.65	0.85
T4		6.463	6.172	6.621	4.66	0.85
T5		6.348	6.006	6.198	4.45	0.56
T6		8.103	5.998	8.980		0.85
T7		6.667	5.972	7.398	4.88	0.85
N17		5.794	5.907	6.581	4.40	0.85

3.7 In silico pharmacokinetic ADME, and toxicity results

The results of the ADME parameters and the toxicity profile of the candidate compounds and the N17 molecule are summarized in Tables 6 and 7 respectively. The estimated Caco/2 permeability value was very high for all the compounds tested (greater than 0.64), the intestinal absorption was greater than 90 %, indicating the good absorption of these compounds(Tabti et al., 2023b). The human VDss was low (less than 0.71 log (L/kg)(Hosen et al., 2023). At the level of penetration, the compounds screened revealed low penetration of BBB and CNS. In addition, these compounds examined, proved to be no-cytochrome P450 (CYP/2D6) inhibitors and P450 (CYP3D4) inhibitors when they may be a CYP3A4 substrate and not a CYP/2D6 substrate estimated total clearance results showed acceptable results.

The boiled Egg plot between TPSA and Log P to predict cerebral penetration and gastrointestinal absorption of selected molecules is shown in Fig. 11. It can be seen from the plot that compound N17 exhibits the BBB but it was located in the GI absorption region. All newly designed molecules exerted a high HIA. Prediction of substrates (PGP) and non/substrates (PGP/) of the permeability glycoprotein (PGP) were also estimated by the boiled egg model.

Close

Fig. 11

Boilled /Egg plot for all designed compounds (BBB – Blood/brain barrier; PG+// – P/glycoprotein; HIA – human intestinal absorption).

Export to PPT

At the level of the estimation of the toxicity profile, the results obtained are collected in Table 7. The Ames test revealed that all the candidates could be considered as non/mutagenic. The maximum tolerated doses of these designed compounds were lower than those of Compound 17 (0.7 log mg/kg/day). All candidates showed oral LD₅₀ values in rats of between 2.614 and 2.291, which are higher than that of compound N17 (2.239). For the chronic LOAEL model in rats, the compounds tested showed LOAEL values between 1.21 mg/kg-bw/day and 0.393 mg/kg-bw/day. All candidates were evaluated as non-skin sensitization.

3.8 Docking results

In order to gain structural insight on the binding mechanism of all newly predicted compounds into the protein binding site, molecular docking analyzes were performed using autodock tools software. The crystal structure of the ATPase domain of human DNA topoisomerase IIα (PDB:1ZXM) was chosen. In addition we also docked Compound N17 and Doxorubicin, to compare our results. Table 8 summarizes all the results obtained. Generally the analysis of the docking results of the newly designed molecules revealed very high docking scores between −9.6 and −10.6 kcal/mol in interval), which shows the good positioning of these compounds in the active receptor site. The best engineered molecule showed −10.3 kcal/mol higher affinity than Compound N17 (-9.6 kcal/mol) and Doxorubicin (-10.1 kcal/mol).

The results obtained show that the residues Asn91, Ser-148, Asn-150, Arg-162, Ile-125, Ser-149, Gly-164 and Lys-168 are the most frequently interacting with the topoisomerase II residues of compound T1 and ANP. The 3D visualization of the interactions of the molecules T1, 17 and Doxorubicin has been displayed in Fig. 12. Compound T1 forms a higher number of conventional hydrogen bonds (eight bonds) with the following residues Thr-165, Asn-163, Gln-376, Gly-166, Gly-164, Arg-162, Ser-148, and Asn-150, and two non-classical hydrogen bonds with Asn-91, Asp-94, while the hydrophics interaction are absent. The highest compound in the data set created three conventional hydrogen bonds Asp-94, Asn-150 and Lys-168. more hydrophobic interactions are formed with Phe-1, Arg-98, Ile-141, Phe-142, and Thr-159. Concerning the control compound Doxorubicin produced four hydrogen lisions with Asn-91, Asn-95, Arg-98, and Thr-215.

Close

Fig. 12

3D visualization of the main interactions between the active site of Human DNA topoisomerase II α (PDB ID: 1ZXM) with the compounds T1, N17, ANP, and Doxorubicin.

Export to PPT

Close

Fig. 12

3D visualization of the main interactions between the active site of Human DNA topoisomerase II α (PDB ID: 1ZXM) with the compounds T1, N17, ANP, and Doxorubicin.

Export to PPT

To validate the docking method to use, were docked the ligand ANP in its active site and by calculating the RMSD parameter, the results obtained were displayed on the structure of Fig. 13 showing a low RMSD value of 1.09 Å, this which indicates The results indicate the reliability of the molecular docking used to predict binding affinity for ligands.

Close

Fig. 13

3D visualization of docking validation by the redocking method of the co-crystallized ligand in the topoisomerase IIα receptor (PDB: 1ZXM).

Export to PPT

3.9 MD simulation results

In order to further analyze the atomic behavior of molecular interactions dynamically contributing to their stable complexed conformation and to further explore the binding mode of the complexes resulting from molecular docking. The docking results of the complex bound to compound T1, 17 and doxorubicin was simulated along a 100 ns trajectory. Judgment of the stability of the screened complexes was made based on parameters including Root-Mean-Square Deviation (RMSD), the root mean square fluctuation (RMSF).

The backbone RMSD plot of three simulated systems versus time has been shown in Fig. 14 (A) the RMSD value of protein bound to compound 17 and doxorubicin showed an almost similar evolution where a gradual increase in the range of 0 ns-18 ns, then stabilization at a relatively constant level for the remaining time around 2.1 Å, these indicating that the conformation of the complexes relative to compound N17 and Doxorubicin were stable throughout the MD simulation trajectory. Concerning the complex linked to compound 17, after 3 ns the RSMD values increased in the first 3 ns, then showed fluctuations in a range between 1.5 Å and 2.5 Å, indicating that the system had reached a stable state. But it should be noted that in the last 15 ns the value of RMSD increases sharply and stabilizes around 2.9 Å. These results showed more stability of the newly designed compound T1 compared to the most active compound N17 of the data set.

Close

Fig. 14

(A) RMSD plot. (B) RMSF plot of the complexes screened during 100 ns MD simulation and (C) Radius of gyration plot of the complexes screened during 100 ns MD simulation.

Export to PPT

The average residue position fluctuation is calculated by root mean square fluctuation values (RMSF) to assess the flexibility of protein chain residues during of the simulation MD. The RMSF plot of protein bound to the screened ligands is shown in Fig. 14(B). It can also be observed that the evolution of RMSD for the three systems was almost similar and also showed the same behavior.

Generally, three ranges of residues were recorded with a fluctuation greater than 3 Å, are from Tyr-151 to Glu-155, from Glu-281 to Glu-283 and from Lys-342 to Asn-345, while the rest of the protein backbone residues have weak and insignificant fluctuations.

The radius of gyration (Rg) is a parameter that expresses the compactness and shape of the complex protein structure, and is therefore an indicator of the stability of the ligand. As shown in Fig. 14(C), the protein ligated with the compound designed by our model showed almost stable Rg values around 3.65 Å during the last 80 ns of simulation. We can observe that the control compound, doxorubicin, also showed stability throughout the simulation but at a value of around 3.9 Å. On the other hand, the Rg of the protein bound to compound 17 initially stabilized around 3.65 Å, but progressively decreased after 75 ns and stabilized around 3.4 Å.

3.10 PCA and FEL analysis

The statistical technique known as Principal Component Analysis (PCA) serves to simplify intricate datasets by reducing their dimensionality. It has utility in investigating changes in protein conformation and identifying significant motions in molecular dynamics simulations. Meanwhile, Free Energy Landscape (FEL) analysis provides a valuable approach for delving into the thermodynamics and kinetics of protein–ligand interactions. This analysis yields insights into the binding process and can even predict the affinity of a ligand for a protein.

The outcomes of energy landscape analysis, as revealed in Figs. 15 and 16, unveil a diverse array of characteristics. Specifically, the presence of two energy wells of varying sizes in the first compound (T1), as illustrated in Fig. 15, might signify conformational stability and reduced flexibility. Likewise, the existence of a solitary energy well in the second compound (N17) (Fig. 15 (b)) suggests a more stable and less flexible conformation. Finally, in Fig. 16, two medium-sized energy basins are evident. In light of the FEL and PCA analyses, complex studies may exhibit increased stability and reduced flexibility. In particular, the T1 complex, characterized by either fewer energy wells or a larger energy well, may be indicative of a more stable and less flexible conformation.

Close

Fig. 15

Plots pca and fel for (a) compound t1 and (b) compound 17.

Export to PPT

Close

Fig. 16

Plots pca and fel for doxorubicin.

Export to PPT

3.11 Free binding energy

The free energy of binding is a dominant factor in estimating the activity of drug candidates. As much as, the lower the value of G_bind, the more stable the complex. The validation of the binding affinity of the compounds T1, 17 and doxorubicin with the protein binding site was carried out by calculating the free energy of binding by the MM-GBSA method. The results obtained have been gathered in Table 9 for the beginning (0 ns) and the end of simulation (100 ns).

The binding free energy of the compounds Pred01, 17, Doxorubicin were at the beginning of the simulation −63.0195 kJ/mol, −87.6201 kJ/mol, and −62.7450 kJ/mol, and at the end of the simulation trajectory −75.2439, −65.1043, and −85.7266. It is interesting to note that the free energy and decreases during the simulation for the compound T1 and Doxorubicin. Which reflects the stability of interactions has been improved in a dynamic environment. Van der Van der Waals energy (E_vdW) and Energy of lipophilicity (E _lipo) and Couloubian energy (ΔG_Coulomb) contributed much more to the binding free energy than the other energies, indicating that the hydrophobic interaction and electrostatic interaction Coulomb type played a key role in the complex system. Although the interaction of Hydrogen bond energy (ΔG _H-bond) and packing energy (ΔG_{Bind packing}) were relatively less important. However, the positive values of the polar solvation energy (ΔG_SolvGB), indicating that it was not conducive to ΔG_bind, the reason was that an excessive binding pocket could cause exposure of the ligand to the solvent.

3.12 Reactivity cycloaddition 1,3 of compound T1

This section focuses on the examination and evaluation of the mechanism behind the 1,3-dipolar cycloaddition reaction of T1, a newly designed compound from mansonone azide and 3,3,3-trifluoroprop-1-yne (illustrated in Fig. 17). For the study of the reaction, two of the regioselective attacks were considered between the reactants and their theoretical parameters were calculated using density functional theory (DFT) calculations performed at the B3LYP/6-31G (d, p) theoretical level.

Close

Fig. 17

Cycloaddition 1.3-dipolar reaction between 3,3,3-Trifluoroprop-1-yne and Mansonone azide.

Export to PPT

3.12.3

3.12.3 Theoretical elucidation of the regioselectivity of the reaction by

In order to highlight the preferential mode of cyclization (isomer 1 or isomer 2) and consequently the major product of the cycloaddition reaction of mansonone triazole methyl trichloride, we located the transition states and calculated the activation barriers for the two possible reaction paths.

The transition states TS1 and TS2, corresponding to the two regioisomers Isomer 1,4 and Isomer 1,5 respectively, have been located at the computational level B3LYP/6-31G*(d,p). The two transition states were confirmed by the presence of one and only one eigenvalue negative in the matrix of force constants, i.e. there is a single imaginary frequency in the Hessian matrix and it corresponds to the mode of vibration corresponding to the formation of the two new bonds. The calculations yielded, as expected, analogous energy barriers for regiochemistry 1,4 and 1,5 (Fig. 18), resulting in 17.3842 and 15.9698 kcal mol-1, respectively. This energetic deference explains that isomer 1 is the favorable isomer for formation. The optimized transition state structures are shown in Fig. 18. Additionally, an analysis of the geometries at the TS is given in Table 12 and Fig. 19.

Close

Fig. 18

The activation energy barriers for cycloaddition of a mansonone azide with 3,3,3-Trifluoroprop-1-yn.

Export to PPT

Table 12 Optimized geometry parameters for the TS 1,4 and TS 1,5 in concerted mechanism for reaction.

Optimized Parameters
TS 1,4			TS 1,5
Distance bond	Value(Angstroms)		Distance bond		Value(Angstroms)
R(16,17)	1.262		R(1,2)		1.234
R(16,26)	2.7504		R(1,3)		1.4674
R(16,27)	2.2106		R(1,23)		2.2012
R(17,18)	1.1654		R(1,24)		2.6266
R(17,26)	2.5855		R(1,25)		2.7329
R(17,27)	2.6404		R(2,7)		1.0692
R(18,26)	2.1815		R(2,23)		2.7526
R(18,27)	2.7779		R(2,24)		2.5831
R(26,27)	1.235		R(2,25)		2.1418
Angles	Value(Degrees)		R(20,23)		1.4639
A(13,16,17)	119.1223		R(23,24)		1.2595
A(13,16,27)	131.6954		R(24,25)		1.1698
A(17,16,32)	119.8582		Angles		Value(Degrees)
A(26,16,32)	50.8633		A(2,1,3)		157.5263
A(16,17,18)	140.3347		A(1,2,7)		160.0067
A(18,26,28)	100.9443		A(20,23,24)		121.2944
A(27,26,28)	153.631		A(23,24,25)		139.0246
A(16,27,18)	53.0167
A(17,27,32)	123.1696
A(18,27,32)	147.8052
A(26,27,32)	162.9151

Close

Fig. 19

Structures of transition states TS1,4 and TS1,5 of cyclization, the distances are given in Å.

Export to PPT

The 1,3-dipolar cycloaddition of 3,3,3-Trifluoroprop-1-yne was investigated using DFT at the B3LYP/6-31G*(d,p) level. The energy barriers for the coupling of mansonone triazole and 3,3,3-Trifluoroprop-1-yne have been calculated. As shown in Table 13, the activation barriers for the reaction lead to both the 1,4 and 1,5 regioisomers of 1,2,3-triazole, which are very close, respectively with a preferred regioisomer 1,4.

Table 13 Total energies (a.u), energy activation, and frequency of the stationary points involved in the 1, 3-DC reaction.

Reactants	Isomer 1,4	Isomer 1,5
E Dipole (ua)	−3422.4892	−3422.4892
E Dipolarophile (ua)	−414.3898	−414.3898
Product	−3836.93044	−3836.879558
E(TS) ua	−3836.8534	−3836.8513
Energy activation Kcal/mol	15.9698	17.3842
Frequency (cm-1)	−392.18	−387.12
Δd	0.3	0.6

For TS1 and TS2, the binding orders are quite different (2.18 Å −2.21 Å) and (2.14 Å −2.20 Å). This means that the two processes are asynchronous (formation of a bond is advanced compared to the second). However, the process passing through TS1,4 is slightly more asynchronous (Δd = 0.3) than that passing through TS1,5 (Δd = 0.6). Moreover, the very low energy difference between cycloaddition isomer 1,4 and isomer 1,5 shows that these two regioisomers are isoenergetic.