Design of new α-glucosidase inhibitors through a combination of 3D-QSAR, ADMET screening, molecular docking, molecular dynamics simulations and quantum studies

2.1 Selection of dataset

Based on previous research, thirty-one molecules of 4-amino-1,2,4-triazole have a potent and selective α-glucosidase inhibitory activity selected to investigate the bioactivity against DM diseases by theoretical methods (Yeye et al., 2020). The target derivatives were subjected to a 3D-QSAR analysis. They were divided into two groups: twenty-five molecules were randomly selected with their activity as a training set to create a quantitative model. At the same time, six molecules were retained as a test set to check the efficiency of the recommended models. The chemical structures of the selected molecules and their corresponding pIC₅₀ values (pIC₅₀ =  − log IC₅₀), where IC₅₀ shows how much medication is needed to inhibit a biological process by half, are shown in Fig. 3 and Table 1.

Close

Fig. 3

4-amino-1,2,4-triazole analogues structure (1–31).

Export to PPT

Table 1 Chemical structures of 4-amino-1,2,4-triazole derivatives and their pIC₅₀ values.

N°	Structure	pIC50	N°	Structure	pIC50
1a		5.582	17		5.198
2		5.264	18		5.276
3		5.264	19		5.485
4		5.270	20		5.599
5		5.241	21a		5.625
6		5.369	22a		5.258
7		5.368	23a		5.461
8		5.263	24a		5.633
9		5.271	25a		5.644
10		5.262	26		5.360
11		5.370	27		5.636
12		5.196	28		5.264
13		5.369	29		5.184
14		5.604	30		5.267
15		5.270	31		5.201
16		5.607	aTest set compounds

2.2 Geometry minimization and molecular alignment

This paper carefully sketched the derivatives of 31 4-amino-1,2,4-triazole molecules by using SYBYL-X 2.0 package. Next, they have been optimized by utilizing the Tripos force field (Clark et al., 1989) with gradient convergence criteria of 0.01 kcal/mol (Purcell and Singer, 1967). The Gasteiger-Hückel atomic partial charges were calculated and added to the 31 targeted molecules using Powell's approach (Tsai et al., 2010). Molecular alignment is a pivotal step employed to implement the 3D-QSAR models. It was done using a common core and compound C25, the database's more active molecule, as a reference structure. The ALIGN DATABASE algorithms (Bouamrane et al., 2022a) reachable in SYBYL-X 2.0 package was used for the alignment of the 31 molecules under study.

2.3 CoMFA and CoMSIA studies

The CoMFA and CoMSIA approaches were conducted primarily to create a 3D QSAR model that might be employed to estimate the α-glucosidase activity of new proposed inhibitors. Both techniques are built on the concept of fields surrounding aligned molecules (Khaldan et al., 2022a-d). First, the CoMFA method focuses on calculating steric (S) and electrostatic (E) energies by using a carbon atom hybridized sp³ and a net charge of + 1.0, with an energy cutoff computation default value of 30 kcal/mol (Cramer et al., 1988). Second, the CoMSIA technique was performed to complement the CoMFA technique by calculating more physicochemical fields as hydrophobic (H), hydrogen bond donor (D), and acceptor (A), as well as electrostatic and steric fields (Klebe el al., 1994) and keeping the same parameters as those used in CoMFA. During the computation of CoMFA and CoMSIA techniques, the column filtering was fixed at 2.0 kcal/mol, and the correction factor was set at 0.3, the default value (Ståhle and Wold, 1988).

2.4 3D-QSAR models production

The 3D-QSAR models were set up using Partial Least Squares (PLS) (Wold, 1991) technique. Numerous statistical parameters that characterize the models were identified using leave-one-out (LOO) cross-validation and non-cross-validation techniques available in PLS approach. The LOO cross-validation technique was conducted to yield the cross-validation correlation (Q²) and the optimal number of components (N) parameters. On the other hand, the non-cross-validation technique was done to produce the correlation coefficient (R²), F-test value (F), standard error of estimate (SEE), as well as the five fields. The values of Q² (>0.5) and R² (>0.6), as well as the lower value of SEE, helped to identify the optimal CoMFA and CoMSIA models. Moreover, the external validation of selected models was executed to assess their proficiency and predictive power using six 4-amino-1,2,4-triazole molecules as a test set (Bouamrane et al., 2022b; Khaldan et al., 2021a-d; Bhatia et al., 2017).

2.5 Validation of 3D-QSAR models

One of the main purposes of using a 3D-QSAR model is to forecast the α-glucosidase activity of non-synthetic molecules evaluated by their robustness and predictive ability. For this reason, the recommended models have been submitted to additional validation. Six compounds have been employed in the external validation approach as test set molecules. Similarly to the training set, these molecules were sketched and optimized by the same method. The α-glucosidase activity of the six studied molecules used in the validation process was calculated using the recommended models. The prediction capacity of the elaborated models was investigated by determining the external validation coefficient (r²pred) that computed using the equation (1) (Golbraikh and Tropsha, 2002):

Where PRESS is the quadratic difference between the actual and estimated activity values of the compounds in the test set, and SD is the total of the squared deviations between the activity values in the test set and the mean activity values of the training set. Furthermore, Golbraikh and Tropsha (Golbraikh and Tropsha, 2002) calculated additional statistical parameters such as coefficients ${\mathrm{r}}_0^2$ and ${{\mathrm{r}}_0^{\prime }}^2$ , slope k and k' to ensure the strength of the established models using the equations (2), (3), (4) and (5), respectively:

Where ${\mathrm{r}}_0^2$ and ${{\mathrm{r}}_0^{\prime }}^2$ are squared correlation coefficient values of estimated vs actual, and actual versus estimated activity for the test set compounds at zero intercept, respectively.

k and k' are the slopes of the estimated against observed and observed vs estimated activity plots at zero intercepts for the test set, respectively.

A study by Roy (Roy, 2007) has recently demonstrated that it is necessary to compute the parameters ${\mathrm{r}}_{\mathrm{m}}^2$ and ${{\mathrm{r}}_{\mathrm{m}}^{\prime }}^2$ , which are the difference between $r^2$ is and ${\mathrm{r}}_0^2$ values, and $r^2$ and ${{\mathrm{r}}_0^{\prime }}^2$ values, respectively. The model satisfied Roy's criteria can be considered as a robust model and might be applied to estimate the activity of the new suggested inhibitors. The parameters ${\mathrm{r}}_{\mathrm{m}}^2$ and ${{\mathrm{r}}_{\mathrm{m}}^{\prime }}^2$ are calculated by applying equations (6) and (7), respectively:

2.6 Y-randomization test

The Y-randomization test (Rücker et al., 2007) is another validation approach used to evaluate the capabilities of the recommended models and to obviate the problem of overfitting in the training set. Indeed, the approach randomly permuted the pIC₅₀ values in the training set. Thus, a new model was elaborated after each permutation. After the Y-randomization test, the smaller Q² and R² values show the power of the original models, while the higher values show the fragility of the models.

2.7 Design of novel 4-amino-1,2,4-triazole inhibitors and prediction of their α-glucosidase activity

3D-QSAR contours created by the recommended models were studied to determine the regions of the molecule responsible for bioactivity as α-glucosidase inhibitors. Thus, the original molecule can be fully or partially modified to obtain a better activity. This process was done by applying the SYBYL-X 2.0 software. Then, the new molecules were subjected to the same methodology in executing the 31 proposed molecules, such as sketching, minimization and alignment. Finally, α-glucosidase activity of the new targeted inhibitors was predicted using the recommended 3D-QSAR models.

2.8 In silico pharmacokinetics ADME/Tox and drug likeness prediction

The proposed 4-amino-1,2,4-triazole molecules were subjected to further in-depth study using ADME-Tox analysis and drug-likeness prediction, to assess the bioavailability and drug potency. Several pharmacokinetics parameters such as absorption, skin permeation, blood–brain barrier penetration, cytochromes P450 (CYP), toxicity, and synthetic accessibility, as well as a drug-likeness prediction as Lipinski, Veber, and Ghose rules were determined using SwissADME (Daina et al., 2017) and pkCSM (Pires et al., 2015) online servers.

2.9 Molecular docking simulation

Molecular docking methodology is conducted to obtain a greater comprehension of the binding modes and make educated guesses about how a ligand and receptor interact (Shi et al., 2022; Piekus ́-Słomka et al., 2022; Khaldan et al., 2020). In this investigation, molecular docking was carried out by means of two programs; Autodock Vina (Trott and Olson, 2010) and Autodock tools 1.5.6 (Hunter et al., 2001). The crystal structure of isomaltase from Saccharomyces cerevisiae was downloaded from the protein data bank (PDB ID: 3a4a) (Khaldan et al., 2022a; Rahim et al., 2015) with a resolution of 1.60 Å. Using the software Discovery Studio 2016 (Dassault Systèmes BIOVIA, 2016), the retrieved receptor was created by eliminating all of the molecules of water and introducing the polar hydrogen atoms. The active site has been identified, and its coordinates are as follows: x = 21.595, y = -7.436 and z = 24.042. Within the 3a4a receptor's pocket, a grid box with the coordinates x  = 30, y = 30, and z = 30 was made, with a 1 Å space between grid points. An extended PDB format, named PDBQT, is applied to coordinate files that include atom types and atomic partial charges using the Autodock 1.5.6 tools. The non-bonded rotation and flexible ligands are assigned using torsion angles. The DFT approach was used to optimize the C25 and M1 derivatives, which are the subject of a molecular docking study, with the B3LYP/6–311++G (d,p) level and employing the Gaussian G09 software (Frisch, 2009). Then, they are docked in the active site of the targeted receptor. Finally, the derived findings were evaluated using the tools PyMol (DeLano, 2002) and Discovery Studio 2016 (Dassault Systèmes BIOVIA, 2016).

2.10 Molecular dynamics (MD) simulations

In this approach, the web-based CHARMM-GUI was executed to build the system (Jo et al., 2008; Brooks et al., 2009; Lee et al., 2016) interface with the CHARMM36 force field (Best et al., 2012). The NAMD 2.13 (Phillips et al., 2005) package was executed to achieve all MD simulations. The periodic boundary and the TIP3P explicit solvation model were used (Jorgensen et al., 1983). It was built on the basis of web-based technology. With a dimension of the dimensions, conditions were fixed at 118.76 Å, 147.65 Å, and 133.50 Å in x, y, and z, respectively. The CHARMM general force field was executed to generate the parameters for the best docking findings (Yu et al., 2012). Then, the system used was neutralized with 9Na⁺ ions. The minimization, equilibration, and production strategies were part of the molecular dynamic strategies. The present MD simulations were run with an integration time step of 2 fs. The isothermal-isobaric (NPT) ensemble was employed for the generation, whereas the canonical (NVT) ensemble was applied for the equilibration. Throughout the entirety of the MD simulation's 100 ns, the Nose'-Hoover Langevin piston barostat (Nosé and Klein, 1983; Nosé, 1984) was applied to set the pressure to 1 atm, with a period of 0.1 ps and a Langevin piston decay of 0.05 ps. The Langevin thermostat (Grest and Kremer, 1986) set temperature at 298.15 K. With pair list misbehaviors 16 and a distance cutoff of 12.0, nonbonded interactions across short distances were smoothly truncated, while Lennard Jones interactions were easily shortened at 8.0. The particle-mesh Ewald (PME) approach (Darden et al., 1993; Essmann et al., 1995) was utilized to address long-range electrostatic interactions. Grid spacing of 1.0 was applied to arrange all simulation cells. The SHAKE technique was implemented to limit all covalent bonds that contained hydrogen atoms (Ryckaert et al., 1977). Lastly, the same process was adopted for all MD simulations to maintain uniformity.

2.11 Binding energy computations

The relative binding energy computations were executed using the one-average molecular mechanics generalized Born surface area (MM/GBSA) (Genheden et al., 2012; Wang et al., 2019) technique implemented in the MOLAICAL code (Bai et al., 2021), where the ligand (L) links to the receptor (R) to create the complex (RL) according to the equation (8):

This can be expressed by several interactions' contributions,

Where the following gas-phase molecular mechanics changes $\left({\mathrm{\Delta }\mathrm{E}}_{\mathrm{M}\mathrm{M}}\right)$ , solvation Gibbs energy $\left({\mathrm{\Delta }\mathrm{G}}_{\mathrm{S}\mathrm{o}\mathrm{l}}\right)$ , and conformational entropy $\left(-T\mathrm{\Delta }S\right)$ are calculated according to the equation (9): ${\mathrm{\Delta }\mathrm{E}}_{\mathrm{M}\mathrm{M}}$ is the total alteration of the electrostatic energies ${\mathrm{\Delta }\mathrm{E}}_{\mathrm{e}\mathrm{l}\mathrm{e}}$ , the internal energies ${\mathrm{\Delta }\mathrm{E}}_{\mathrm{i}\mathrm{n}\mathrm{t}}$ (bonded interactions), and the van der Waals energies ${\mathrm{\Delta }\mathrm{E}}_{\mathrm{v}\mathrm{d}\mathrm{W}}$ ; ${\mathrm{\Delta }\mathrm{G}}_{\mathrm{S}\mathrm{o}\mathrm{l}}$ is the combined of both the polar solvation (determined by the generalized Born model) and the nonpolar solvation (determined by the solvent-accessible surface area) and $-T\mathrm{\Delta }S$ is derived by the standard mode investigation, however, since we are only concerned with relative binding energies, this part was left out to reduce the computational cost. During the MM/GBSA computations, the surface tension constant of 0.03012 kJ mol^-1Ų was used.

2.12 Global indices of the 4-amino-1,2,4-triazole molecules using DFT approach

The minimum-energy geometries of the investigated compounds were achieved using the Density Functional Theory (DFT) technique and the Gaussian G09 package (Frisch, 2009) by using B3LYP/6–311++G (d,p) basis set. The chemical potential (μ), chemical softness (S), global electrophilicity (ω), chemical hardness (η), and global nucleophilicity (N), called global reactivity indices, were calculated for the new reported inhibitors. They were used to figure out the most nucleophilic and electrophilic inhibitors for the compounds under study. Practically, μ index was derived according to the frontier molecular orbital LUMO and HOMO by using the equation (9). The chemical softness (S) and chemical hardness (η) were computed using the expressions (10) and (11) (Pearson and Songstad, 1967), respectively. The global electrophilicity (ω) (Parr et al., 1999) and global nucleophilicity (N) (Domingo et al., 2008) were determined by means the equations (12) and (13), respectively.

Since it has a lower HOMO energy value, tetracyanoethylene (TCE) was chosen as a computation reference(Domingo et al., 2008; Domingo and Pérez, 2011).

2.13 Conceptual density Functional (CDFT) approach

Conceptual Density Functional Theory (CDFT) and backing of the DFT are frequently used to investigate compounds' local reactivity (Berk et al., 2022). The Parr functions P(r) which contain the nucleophilic ${\mathrm{P}}_{\mathrm{k}}^{-}$ and electrophilic ${\mathrm{P}}_{\mathrm{k}}^{+}$ are valuable tools for determining local reactivity of organic compounds. They were calculated using the equations (14) and (15), respectively:

Where ${\mathrm{\rho }\mathrm{s}}^{\mathit{rc}}\left(r\right)$ is the radical cation's atomic spin density (ASD) at the r atom of a particular chemical, and ${\mathrm{\rho }\mathrm{s}}^{\mathit{ra}}(r)$ is the radical anion's ASD at the r atom.

The local nucleophilic ${\mathrm{P}}_{\mathrm{k}}^{-}$ and electrophilic ${\mathrm{P}}_{\mathrm{k}}^{+}$ Parr functions for the neutral compound are provided by each ASD collected at different atoms of the radical anion and radical cation, respectively.

The local electrophilicity ${\mathrm{\omega }}_k$ and local nucleophilicity ${\mathrm{N}}_k$ indices were recalculated using the equations (16) and (17), respectively, based on the notions of ${\mathrm{P}}_{\mathrm{k}}^{-}$ and ${\mathrm{P}}_{\mathrm{k}}^{+}$ Parr functions:

Where the N and $\mathrm{\omega }$ were obtained by the equations (4) and (5), respectively.

3.1 Molecular alignment

Thirty-one of the 4-amino-1,2,4-triazole inhibitors were aligned using a common core, where compound 25 is the most active in databases. The superposed structures are presented in Fig. 4.

Close

Fig. 4

Overlay and alignment of thirty-one of the proposed 4-amino-1,2,4-triazoles.

Export to PPT

3.2 3D-QSAR models results

A value of cross-validated Q² ˃ 0.50 hinted at the acceptable internal predictive capability as well as the strength of the model (Ouabane et al., 2023; Khaldan et al., 2021), but it is not absolute proof. The best method for evaluating a model's predictive ability is to assess it using test set molecules, which is namely external validation and characterized by r²pred which must be ˃ 0.60 (Golbraikh and Tropsha, 2002).

Combinations between steric and electrostatic fields resulted in a CoMFA model. In contrast, the combinations between five fields, electrostatic (E), steric (S), H-bond donor (D), H-bond acceptor (A), and hydrophobic (H), led to the construction of 31 CoMSIA models, as presented in supplementary Fig. 5. The PLS statistical indicators of the CoMFA and CoMSIA models are listed in Table 2. Regarding the CoMFA model, Table 2 indicates that this model has an essential Q² value of 0.504, a significant R² value of 0.911, F value of 51.248, four optimum numbers of components, and a standard error of estimate (SEE) value of 0.046. The CoMFA model's steric and electrostatic contributions were 0.558 and 0.442, respectively, showing that the steric contribution is more significant than the electrostatic one. The external capacity of the CoMFA model was verified by an external validation characterized by the ${\mathrm{r}}^2\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}$ parameter. Its value was 0.710, i.e. higher than 0.6, which proves the high sturdiness of the CoMFA model. In the case of the CoMSIA model, many models were developed by combining five fields. The highest values of Q², R², ${\mathrm{r}}^2\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}$ , and the smallest value of SEE were considered for choosing the optimal CoMSIA model. The results in Table 2 show that the CoMSIA model with the electrostatic (E), Hydrophobic (H), and H-bond acceptor (A) fields is the best. So, the CoMSIA/EHA model present a correlation coefficient R² value of 0.938, cross-validated coefficient Q² value of 0.66, SEE value of 0.038, and F value of 75.805 in addition to four optimum number of components. Additionally, the CoMSIA/EHA model has been externally validated; the ${\mathrm{r}}^2\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}$ value was 0.781; indicating the high strength of the CoMSIA/EHA model. Furthermore, the electrostatic, hydrophobic, and H-bond acceptor contributions of the CoMSIA/EHA model were 0.330, 0.349, and 0.320, respectively, demonstrating that all three fields will have the same impact on the model. The observed α-glucosidase activity of the 31 4-amino-1,2,4-triazole molecules and their calculated pIC₅₀ using the advised models are clarified in Table 3. The residual between the two activities is less than 1, indicating a strong correlation between them (Table 3 and Fig. 6).

Close

Fig. 5

Q² and R² values of the 31 potential combinations fields of CoMSIA model (S: steric, E: electrostatic, H: hydrophobic, D: H-bond donor, A: H-bond acceptor).

Export to PPT

Close

Fig. 6

Predicted and observed pIC₅₀ plots of the investigated molecules along with their residuals using the CoMFA and CoMSIA/ EHA models.

Export to PPT

The graphs of the actual and calculated α-glucosidase activity values of the all molecules involved in creating the CoMFA and CoMSIA/EHA models are exhibited in Fig. 6. The training set molecules are denoted by blue circles, while red triangles represent the test set compounds. Fig. 6 shows that the solid points of the blue circles and red triangles are near to the Y = X line, advising that the observed and calculated α-glucosidase inhibitory activity of the 31 studied inhibitors exhibit a strong linear relationship.

3.3 CoMFA and CoMSIA/EHA models validation

The established CoMFA and CoMSIA/EHA models were subjected to further external validations using Golbraikh, Tropsha, and Roy criteria, as well as six compounds to evaluate the molded models' predictability and dependability power. The obtained outcomes are registered in Table 4.

The Golbraikh and Tropsha test outcomes show that CoMFA and CoMSIA/EHA models had $r_0^2$ and ${r_0^{\prime }}^2$ values of 0.945, 0.950, and 0.968, 0.966, respectively (Table 4). They also have a $r^2$ pred value of 0.710 and 0.781, respectively. The plot slope values $k$ and $k^{\prime }$ of the elaborated models are in the reliable range. The gap between $r^2$ and $r_0^2$ values divided by $r^2$ is too a criterion of Golbraikh and Tropsha, and its value was found of −0.343 for CoMFA and −0.216 for CoMSIA/EHA models. The identical procedure was again repeated by changing the $r_0^2$ value to the ${r_0^{\prime }}^2$ value and the generated results were still in a reliable range (<0.1). Therefore, the Golbraikh and Tropsha results demonstrated the good quality of the recommended CoMFA and CoMSIA/EHA models. On the other side, the recommended CoMFA and CoMSIA/EHA models met some Roy criteria, and the proposed CoMSIA/EHA model remains the best (Table 4). The derived outcomes confirm the good stability and predictability of the CoMSIA/EHA model for predicting the novel α-glucosidase activity of the proposed 4-amino-1,2,4-triazole inhibitors.

3.4 Randomization approach result

The Y-randomization approach was implemented to assess the models' strength and capability to predict the α-glucosidase activity of innovative 4-amino-1,2,4-triazole compounds. Table 5 demonstrates that, after applying the Y-randomization method, the Q² and R² values for the CoMFA and CoMSIA/EHA models are significantly lower than the original models, proving that the advised models are trustworthy and could be approved to predict the α-glucosidase activity of new title inhibitors.

3.5 3D-QSAR contour maps

CoMFA and CoMSIA/EHA contour plots were created to describe and identify the impact of electrostatic, steric, hydrophobic, acceptor and donor hydrogen bonds on α-glucosidase activity of the investigated compounds. Therefore, identify the activity's preferable and unfavorable substitutes. The obtained outcomes are displayed in Figs. 7 and 8.

Close

Fig. 7

Contour maps results of the CoMFA elaborated model. Steric (a), Electrostatic (b).

Export to PPT

Close

Fig. 8

Contour maps results of CoMSIA/EHA model. Electrostatic (a), Hydrophobic (b), and H-bond acceptor (c).

Export to PPT

3.6 CoMFA contour maps results

In CoMFA steric contour plots (Fig. 7a), the green outlines were found around phenyl’s meta and para positions and near the 1,2,4-triazole moiety, suggesting the bulky groups were favored in these regions. Both structures of molecule C10 (pIC₅₀ = 5.260) and molecule C26 (pIC₅₀ = 5.360) have bulky substituents at meta position of the phenyl moiety. On the other side, the smaller yellow contour around the oxygen atom of nitrogen dioxide moiety hints that tiny moieties enable amelioration of the activity of the molecules.

Regarding the CoMFA electrostatic contour maps (Fig. 7b), the large blue contours covering the phenyl and imine groups near the 12,4-triazole moiety indicate that the electro-donating substituents in these places are preferred to achieve improved activity of the molecules. The molecule C14 (pIC₅₀ = 5.604) and molecule C25 (pIC₅₀ = 5.644) that possess electron-donating moieties (Cl and OCH₃) at ortho, meta and para sites of phenyls group are two examples that can explain this observation. Whereas the tiny red contour near oxygen atom of nitro group clearly showed that electron-withdrawing groups at this position would be beneficial for α-glucosidase activity.

3.7 CoMSIA/EHA contour maps results

The results of CoMSIA electrostatic contour plots are depicted in Fig. 8a, where red color (20 % contribution) indicates the positions where electro-attracting groups are favored. However, blue color (80 % contribution) indicates the positions favored by the electro-donor groups. Hence, red-colored portions around ortho, meta, and para positions of phenyl groups and the (–NO₂) group demonstrated that introducing electro-attracting groups could ameliorate the α-glucosidase activity. Blue-colored portion situated in the phenyl moiety's meta position demonstrates that this site is exclusively available to electro-donor groups to increase the potency of 4-amino-1,2,4-triazole compounds (Fig. 8a).

CoMSIA hydrophobic contour maps (depicted by yellow (80 % contribution) and gray (20 % contribution)) were displayed in Fig. 8b. The yellow outlines indicate the regions of the preferred hydrophobic substituents and gray colored contours shown the regions where hydrophobic moieties are unpreferred. Back to the results of Fig. 8b, a yellow contour was observed around para position of phenyl group, which indicates significance of hydrophobic groups on α-glucosidase activity in this region. Similarly, gray contours around (–NO₂) group and meta position of phenyls group demonstrate that hydrophilic entities are needed to improve α-glucosidase activity (Fig. 8b).

CoMSIA of H-bond acceptor contour plots are displayed in Fig. 8c, where magenta-colored contour (80 % contribution) of oxygen atom of the nitrogen dioxide moiety shows that the use of H-bond acceptor substituents is required to enhance the potency of the new title molecules. On the contrary, a red contour (20 % contribution) near ortho position of phenyl groups and oxygen atom of (–NO₂) group suggested that hydrogen bond accepting moieties in these positions are unfavored (Fig. 8c).

3.8 CoMFA and CoMSIA/EHA contour maps summary

The structural features removed from CoMFA and CoMSIA/EHA contour cards were investigated to determine the structural necessities for each position and the specifications for a chemical group to increase the bioactivity of proposed compounds as α-glucosidase inhibitors, as shown in Fig. 9. A bulky and electro-donating substituent such as halogens groups (Cl, Br, and F), substituted alkanes (isopropyl, isobutyl) or amines groups (methanamine, dimethylamine, and trimethylamine) is required for the R1 position in order to ameliorate α-glucosidase inhibitory activity. The R2 group is requested only to be small, electron-withdrawing, hydrophilic and H-bond donor moieties as nitrile group (–CN) or nitrosyl hydride (NHO). Bulky, electro-donating/electron-withdrawing, and hydrophilic groups such as nitrogen dioxide (–NO₂), formaldehyde (–CH₂O) or formamide (–CONH₂) are beneficial for the R3 moiety (Fig. 9). Therefore, all these recommendations will guide to introduce the correct group in the proper position to design new potent α-glucosidase inhibitors with good predicted activity.

Close

Fig. 9

Structure-activity relationship representation.

Export to PPT

3.9 Design of novel 4-amino-1,2,4-triazole inhibitors

The current study aims to design new α-glucosidase inhibitors by investing the different structural features extracted from the recommended CoMFA and CoMSIA/EHA models. In fact, we substituted the nitro group (–NO₂) by nitric oxide (–NO) group in R2 position and we added the proper substituents in R1 and R3 positions according to the contour maps (Fig. 10). Hence, four new α-glucosidase inhibitors were suggested and their α-glucosidase activity were computed by using the proposed CoMFA and CoMSIA/EHA models as displayed in Table 6. Returning to Table 6 findings, the α-glucosidase activity of the recommended new molecules was better than that of the reported 4-amino-1,2,4-triazole analogues. Fig. 10 shows the chemical structures of the new proposed compounds.

Close

Fig. 10

Chemical structures of the new proposed 4-amino-1,2,4-triazole molecules.

Export to PPT

3.10 Pharmacokinetic properties and ADMET results of the new proposed α-glucosidase inhibitors

Drug likeness is described as a complex equilibrium of numerous molecular properties and structural characteristics that decide if a scaffold is a medication or not (Zaki et al., 2018; El-mernissi et al., 2023). Lipinski’s rule of five is probably the most commonly used method, which suggests that molecules with logP < 5, hydrogen bond acceptor number < 10, molecular weight < 500, and hydrogen bond donor number < 5 are deemed the best oral absorption and permeation (Lipinski et al., 1997). In this direction, we calculated numerous molecular properties of the four new proposed molecules using SwissADME (Daina et al., 2017) and pkCSM (Pires et al., 2015) online servers, as clarified in Table 7. The obtained outcomes demonstrate that the title inhibitors have acceptable oral absorption and permeation. Continuously, excellent bioavailability and increased flexibility for effective interaction with a specific binding pocket are both characteristics of a molecule with a total polar surface area (TPSA) of less than 140 Å² and a number of rotatable bonds (nrotb) of less than 10 (El-mernissi et al., 2021).

Interestingly, the four new proposed molecules present good bioavailability and are more flexible (Table 7). On the other hand, the new proposed 4-amino-1,2,4-triazole inhibitors were screened for their synthetic accessibility to identify the ease with which these compounds can be synthesized. Table 7 reveals that all the compounds under study have synthetic accessibility (SA) values near 1 and below 10, indicating that they are all capable of being synthesized (Hasan et al., 2022). Further pharmacokinetic properties of the recommended scaffolds, ADMET properties were determined using SwissADME (Daina et al., 2017) online server. The findings are displayed in Tables 8, 9, 10 and supplementary Fig. 11.

Close

Fig. 11

BOILED-Egg model.

Export to PPT

3.11 Molecular docking results

The molecular docking method defined intermolecular interactions between 4-amino-1,2,4-triazole ligands and target 3a4a receptor. The crystal structure of the studied receptor (PDB: 3a4a) is shown in Fig. 12. In this study, the active pocket of the 3a4a receptor was docked with compound M1, voglibose (α-glucosidase inhibitor used in clinics), and the best-recommended molecule. The results indicate that molecule C25 is the most potent inhibitor in the database. These findings suggest that C25 has the potential to be developed as a novel therapeutic agent for diseases associated with 3a4a receptor inhibition. The binding affinities of the studied molecules, obtained from Autodock Vina program, are listed in Table 11. The summary of inter-molecular interactions between C25/M1, and voglibose molecules with the target protein are shown in Table 12, Figs. 13, 14, and 15, respectively. Autodock Vina results provided nine different poses; the affinity value for the best binding modes for C25, M1, and voglibose molecules were obtained as 6.600, −7.600, and −6.400 kcal/mol, respectively, as shown in Table 11. Hence, M1 molecule might have more diabetes mellitus inhibitory potential than C25 and voglibose molecules. Additionally, M1 compound provided a lower inhibition constant (Ki) than the other studied compounds, demonstrating that a smaller amount of the M1 inhibitor can inhibit diabetes mellitus (Table 11). Whereas, the docking result of C25 compound interactions, a conventional hydrogen bonding interaction was formed between (N) atom of the Arg442 residue and (N) atom of the imine group. The distance for this interaction was obtained as 3.14 Å with ≡N ···N notation. Similarly, a carbon-hydrogen bond was established between = O atom of the Asp352 residue and the aromatic CH of 1,2,4-triazole ring with = O···C–H notation at 3.64 Å value of distance. The Pi-Pi-T –shaped interactions was observed between π electrons of the phenyl group in residues Tyr72 and Tyr158 with π electrons of 1,2,4-triazole group at 4.69 and 4.93 Å values of distances, respectively. The compound C25 provided a pi-anion interaction between = O atom in residue Asp69 and π electrons of 1,2,4-triazole group at a distance of 3.89 Å. On the other hand, the docking interactions of voglibose show more conventional hydrogen bonds with different residues and distances, making it more stable in the active site of protein 3a4a. In the case of molecular docking outcomes of M1 Compound showed a pi-alkyl and five conventional hydrogen bond interactions. These critical interactions make the compound M1 pharmacologically very important compared to the molecule C25 since hydrogen bonding significantly impacts the pharmacological action of ligands. In more detail, a pi-alkyl interaction was observed between Alkyl group in residue Arg315 and π electrons of phenyl group at 4.12 Å value of distance. Two conventional hydrogen bond interactions were made between ≡N atom in residue Asn415 and = O atom of amide group with ≡N···O = notation at 2.92 Å value of distance, and = O atom in the same residue and (–NH) group of amide with = O ···HN notation at 2.77 Å value of distance, respectively. Additionally, the hydroxyl (–OH) group situated at ortho site of phenyl group formed a conventional hydrogen bonding interaction with Tyr158 residue at 1.87 Å value of distance. Two other conventional hydrogen bonding interactions were formed between ≡N atom of the Arg315 residue and = O atom of (NO) with the notation ≡N···O = at a distance value of 3.00 Å, and the = O atom of the Pro312 residue and the = O atom of nitric oxide with the notation = O···O = at a distance value of 3.40 Å, respectively.

Close

Fig. 12

The crystal structure of the studied receptor (PDB: 3a4a).

Export to PPT

Close

Fig. 13

Representation of the 3D (left) and 2D (right) interactions of the molecule C25 with 3a4a receptor.

Export to PPT

Close

Fig. 14

Representation of the 3D (left) and 2D (right) interactions of the voglibose drug with 3a4a receptor.

Export to PPT

Close

Fig. 15

Representation of the 3D (left) and 2D (right) interactions of the molecule M1 with 3a4a receptor.

Export to PPT

The reliability and effectiveness of molecular docking procedures have been confirmed through the re-docking process. As shown in Fig. 16, the re-docked conformation (in magenta) overlaps with the original ligand (in red) with an RMSD (Root Mean Square Deviation) value of 0.35 Å, which lies within the satisfactory range of 2 Å (Kramer et al., 1999). The successful superimposition of the two compounds and the RMSD value below 2 Å confirm the Autodock Vina algorithms' capability to conduct molecular docking protocols with confidence. Additionally, Fig. 16 also displays the result of the original ligand's molecular docking to the active site of the 3a4a receptor. The original ligand formed multiple interactions with varying distance values, including six conventional hydrogen bond interactions, validating its stability in the active pocket of the receptor ().

Close

Fig. 16

A) re-docking pose with the rmsd value of 0.35 Å (Red = Original, Magenta = Docked), b) Docking interactions of the co-crystalized ligand.

Export to PPT

3.12 MD simulations results

3.12.1

3.12.1 RMSD

RMSD was identified for the complex on the basis of 'Backbone' atoms by means of the VMD software. The average RMSD value for C25 and M1 complexes are 4.252 ± 0.9477 Å and 1.856 ± 0.326 Å. The protein complex's RMSD graph demonstrates that, with slight fluctuation in the range of ∼ 1 Å, the structure remained constant during the simulation time, which is a usual behavior of globular protein (Fig. 17). The ligand's RMSD was computed using the VMD algorithm on the basis of the atoms of the ligand. The average RMSD value for C25 and M1 ligands are 1.1616 ± 0.360 Å and 1.8491 ± 0.323 Å. The ligand's RMSD of C25 remained below 2 Å at all times. The fluctuations of amino acids in the protein, both in the presence and absence of the ligand, were indistinguishable until the 42 ns mark. Beyond 42 ns, it was observed that the ligand changed its position or conformation within the same active site, explaining the increased vibrational intensity of the complex up to 85 ns. After this point, the ligand's conformation remained stable but converged to around 4 Å RMSD. Generally, the complex's RMSD exceeded 2 Å, indicating the conformational destabilization of the ligand throughout the dynamic simulation. On the other hand, RMSD of M1 fluctuated until 62 ns, it stabilized at 2 Å, M1 ligand remained bound throughout the simulation (Fig. 17).

Close

Fig. 17

RMSD (left), RMSF (middle) and Rg (right) of the investigated complexes in a MD simulation time of 100 ns. C25 (top) and M1 (bottom).

Export to PPT

3.12.4

3.12.4 Protein-ligand hydrogen bonds

The total number of hydrogen bonds established between C25 and M1 ligands and 3a4a receptor throughout the simulation time’s 100 ns is depicted in Fig. 18. C25 and M1 ligands exhibit a consistent changing number of H-bonds with the protein and remains in bound form. The fluctuations show that the conformation of the C25 and M1 ligands inside the pocket is changing (Fig. 18). Hydrogen bond occupancy values for each 3a4a-C25/M1 complex are shown in Table 13. Hydrogen bond occupancy is the ratio of conformations where a particular residue engages in hydrogen bonds over of 1000 conformations of the 3a4a-C25/M1 complex. Each complex's 1000 conformations were computed using the 100 ns molecular dynamics pathway. The residues with the highest hydrogen bonding interactions are shown in Table 13.

Close

Fig. 18

Protein-ligand hydrogen bonds (left) and average distance between protein and ligand (right) of the complexes studied during the 100 ns simulation. C25 (top) and M1 (bottom).

Export to PPT

Table 13 Details and occupancy of hydrogen bonds.

Complex 3a4a-C25			Complex 3a4a-M1
Donor	Acceptor	Occupancy (%)	Donor	Acceptor	Occupancy (%)
ARG213-Side-NH1	UNK0-Side-N1	1.00 %	ARG442-Side-NH1	UNK0-Side-N2	8.80 %
ARG213-Side-NH2	UNK0-Side-N2	1.40 %	UNK0-Side-N20	SER311-Main-O	7.20 %
GLN279-Side-NE2	UNK0-Side-O15	1.70 %	THR306-Side-OG1	UNK0-Side-O17	16.50 %
HSD280-Main-N	UNK0-Side-O15	1.00 %	UNK0-Side-N20	GLN353-Side-OE1	9.70 %
ARG213-Side-NH1	UNK0-Side-N2	0.80 %	UNK0-Side-O15	ASP307-Side-OD2	8.50 %
ASN414-Side-ND2	UNK0-Side-N2	0.90 %	UNK0-Side-O15	ASP307-Side-OD1	9.90 %
LYS156-Side-NZ	UNK0-Side-N2	0.70 %	UNK0-Side-O15	SER311-Main-O	7.60 %
LYS156-Side-NZ	UNK0-Side-N1	0.60 %	UNK0-Side-N20	PHE303-Main-O	7.90 %
GLN171-Side-NE2	UNK0-Side-O15	0.40 %	ARG442-Side-NH1	UNK0-Side-N1	15.10 %
LYS503-Side-NZ	UNK0-Side-O15	0.30 %	UNK0-Side-N14	THR306-Side-OG1	19.20 %

3.12.6

3.12.6 Contact frequency (CF) analysis

A contact frequency (CF) study was conducted using the contactFreq.tcl module in VMD with a cutoff of 4 to further examine the binding between 3a4a and the C25 and M1 ligands tested; the findings are shown in Fig. 19. The residues with higher CF throughout the simulation time for C25 were Tyr72, Arg213, PHE303, Asp107, Glu277, Gln279 and Ser304. While for M1, the residues were Phe303, Thr306, Gln353, Tyr158, Asp307, Arg315, Tyr347 and Asn350.

Close

Fig. 19

Contact frequency (CF) analysis.

Export to PPT

3.12.7

3.12.7 Principal component analysis (PCA)

PCA of 3a4a-C25/M1 complex was computed using Bio3D software in R as shown in Fig. 20. All three PCs captured 70.1 % for C25 and 47 % for M1 of the structural variance in 3a4a protein.

Close

Fig. 20

PCA (left) and DCCM (right) of the compounds C25 (top) and M1 (bottom).

Export to PPT

3.12.9

3.12.9 Temperature, pressure and potential energy

The NAMD log files were used to generate the temperature, pressure and potential energy of the system throughout the MD simulation time, as shown in Fig. 21. The graph shows converged potential energy, pressure, and temperature.

Close

Fig. 21

Temperature (A), pressure (B) and potential energy (C) throughout the entire simulation of 100 ns.

Export to PPT

3.13 MM/GBSA binding energy

The MM/GBSA approach was put into practice to reevaluate the complexes because it is the quickest force field-based technique to identify free energy in comparison to other methods for calculating free energy, as the thermodynamic integration (TI) or energy perturbation (FEP) approaches. Comparative investigations revealed that the MM/GBSA technique performed better than the MM/PBSA (Molecular Mechanics/Poisson Boltzmann Surface Area) (Hou et al., 2011). MM/GBSA computations were put into effect using MolAICal program (Bai et al., 2021). Exploiting the results of MD simulations performed by the NAMD software, MolAICal provides a convenient way for computing MM/GBSA (Phillips et al., 2005). Table 14 illustrates the computed binding free energies. Table 15 shows the MD outcomes and MMGBSA relative binding energies for the leading docking molecules.

3.14 Global properties of 4-amino-1,2,4-triazole molecules

The global reactivity indices of C25 and M1 molecules were determined as shown in Table 16 to determine the electrophilic and nucleophilic compounds. The electrophilicity ω index and nucleophilicity N index values of molecule C25 are 2.532 and 2.498, respectively. They are both electrophiles and nucleophiles of sufficient strength. On the other side, the electrophilicity index ω and nucleophilicity index N of compound M1 are 4.118 and 2.622, respectively, demonstrating that this compound is considered as a strong electrophile according to the electrophilic scale (Jaramillo et al., 2008) and a good nucleophile according to the nucleophilicity scale (Domingo et al., 2002). Therefore, the title compound M1 is favored over compound C25 concerning reactivity.

Furthermore, the Frontier Molecular Orbitals (FMOs) of molecules C25 and M1 have been analyzed at B3LYP/6–311++G (d,p) level. Fig. 22 shows the FMOs' corresponding energy levels for the compounds investigated. Because it accounts for the final charge transfer interaction in the compound, the HOMO-LUMO energy gap is crucial in determining molecular electrical transport properties (Miar et al., 2021). A compound with a large frontier orbital gap (HOMO-LUMO energy gap), also known as the ${\mathrm{\Delta }\mathrm{E}}_{\mathrm{g}\mathrm{a}\mathrm{p}}$ , has low chemical reactivity since it is energetically inefficient to add an electron to the high-lying LUMO to get an electron from the low-lying HOMO (Aihara, 1999; Ruiz-Morales, 2002). The derived values of ${\mathrm{\Delta }\mathrm{E}}_{\mathrm{g}\mathrm{a}\mathrm{p}}$ indicated significant chemical reactivity of M1 title molecule.

Close

Fig. 22

The geometries of the HOMO and LUMO orbitals, along with the value of the energy gap values. a) compound C25, b) compound M1.

Export to PPT

3.15 Molecular electrostatic potential (MEP)

MEP maps are logically used to obtain information about the areas of nucleophilic and electrophilic reactivity of a molecule. The charge distributions of the molecules were viewed in three dimensions using the MEP surface. Additionally, the binding of a molecule to a receptor's active site highly depends on the three-dimensional distribution of the electrostatic potential (Han et al., 2005). The B3LYP/6–311++G (d,p) basis in the Gaussian software was used to create the MEP surfaces of the investigated compounds. Table 17 displays the MEP maps of the compounds C25 and M1. In MEP contours, the blue color shows the positive (electron-poor) region; the light blue indicates the slightly electron-deficient region, the neutral region is denoted by green, and the red color demonstrates the negative (electron-rich) part, and the yellow color represents the part slightly rich in electrons. Table 17 shows that compound C25 has a high negative potential (red color) at the level of the N1 and N2 nitrogen atoms of the 1,2,4-triazole group and a low negative potential (yellowish color) at the same level. As a result, the overall electron density of compound C25 is localized on the N1 and N2 atoms, which have the highest negative potential and hence are the most attractive targets for electrophilic attack. On the other hand, the compound M1 presents an acceptable positive potential (light blue color) at the level of the atoms C11 and C13 of the phenyl group, indicating that these centers are the most attractive targets for nucleophilic attack.

Table 17 Optimized structure and electrostatic potential surface of the molecules C25 and M1.

N°	Optimized structure	Electrostatic potential surface
C25
M1

3.16 Local properties of 4-amino-1,2,4-triazole inhibitors

The ${\mathrm{P}}_{\mathrm{k}}^{+}$ and ${\mathrm{P}}_{\mathrm{k}}^{-}$ Parr functions for the compounds C25 and M1 were calculated to find out which centers are the most nucleophilic and electrophilic. The outcomes are shown in supplementary Table S18.

Study of the electrophilic ${\mathbf{P}}_{\mathbf{k}}^{+}$ Parr function of the molecule C25 shows that the C8 and C9 atoms are the most electrophilic centers because they have the highest local electrophilicity index values ( ${\mathrm{\omega }}_k$ ). On the other side, the result of the nucleophilic ${\mathrm{P}}_{\mathrm{k}}^{-}$ Parr functions of the same molecule illustrates that N1 and N2 are centers considered as the most electrophilic due to their the highest local nucleophilicity index values ( ${\mathrm{N}}_k$ ). In the case of the title molecule M1, the C11 and C13 atoms with the largest local electrophilicity index values ( ${\mathrm{\omega }}_k$ ) were identified as the most electrophilic centers based on the electrophile ${\mathrm{P}}_{\mathrm{k}}^{+}$ Parr function analysis. While the investigation of nucleophilic ${\mathrm{P}}_{\mathrm{k}}^{-}$ Parr function of the molecule M1 shows that N14 and O15 atoms are the most electrophilic centers since they have the highest nucleophilicity index values ( ${\mathrm{N}}_k$ ).

The most favoured two-center electrophilic/ nucleophilic interaction will therefore occur between the C8 and C9 atoms of compound C25, as well as C11 and C13 atoms of compound M1, and the most nucleophilic centers of the α-glucosidase receptor site.