The studied derivatives of S-substituted 4-R₁-5-(((3-(pyridine-4-yl)-1H-1,2,4-triazole-5-yl)thio)methyl)-4H-1,2,4-triazole-3-thiols, 5-(thiophene-3-ylmethyl)-4-R₁-1,2,4-triazole-3-thiol, 3-(ethylthio)-9-methylpyrazolo[1,5-d][1,2,4]triazolo[3,4-F][1,2,4]triazine-6-alkylthio, 4-amino-5-(3-methyl-1H-pyrazol-5-yl)-4h-1,2,4-triazole-3-alkylthio, 9-methylpyrazolo[1,5-d][1,2,4]triazolo[3,4-F][1,2,4]triazine-3-alkylthio, which are presented in Scheme 1, they have been synthesized by us earlier, and their physico-chemical parameters and structure’s evidence are presented in the works (Karpun and Polishchuk 2020, 2021, Fedotov and Hotsulia 2021; Karpun 2021; Karpun et al. 2021; Khilkovets 2021a, b). Two-dimensional structures of 1,2,4-trazole derivatives were created by means of ChemAxon Marvin Sketch (Csizmadia 1999), optimized at pH 7.4, minimized using the MMFF94 force field (Halgren 1996), and converted to sdf format via OpenBabel-2.4.1 (O’Boyle et al. 2011).

Scheme 1. 

Scheme of the studied derivatives of 1,2,4-triazoles used in computable investigations.

In this study, 6 enzymes were selected from the RCSB protein database that are somehow responsible for controlling oxidative stress in the body: peroxiredoxin (peroxidase) (PDB: 3MNG), – is associated with antioxidant protection; NO-synthase (PDB: :6NGJ); NADPH – oxidase (PDB: 2CDU); tyrosinase: limits the rate of melanogenesis (PDB: 3NM8); NMDA receptor (PDB: 4KFQ); hemoxygenase, is an important stress protein which takes part in cell defense, antioxidant, and anti-inflammatory activity (PDB: 1N3). The receptors were prepared using OpenBabel-2.4.1 and Molegro Molecular Viewer (Thomsen and Christensen 2006), and structural water molecules, cofactors, and ligands were not taken into account. The missing amino acid residues in the structure of 6NGJ (Cys297, Pro298, Ser339, Gln340, His341, Thr342, Arg343, Lys344, Pro345, Glu346, Asp347, Lys717, Gly718) and 4KFQ (Gly1, Met2, Ser3) were modeled by means of Modeller 10.2 (Webb and Sali 2016). It should be noted that these residues are present in the FASTA file, and although they were not a part of the binding site, without their addition to the protein structure, the topology of the complex would be incorrect. The Chimera 1.16 program was used to minimize the structure and add missing hydrogens (Pettersen et al. 2004).

The ADME properties of the selected compounds were predicted by the freely available SwissADME software (Daina et al. 2017). Various molecular parameters such as molecular weight, number of hydrogen bond acceptors, number of hydrogen bond donors and LogP values were analyzed according to Lipinski’s rule of five.

To determine the mechanism of binding and interaction of selected compounds and targets, molecular docking studies were performed using Smina (Koes et al. 2013), which is a fork of Autodock Vina (http://vina.scripps.edu/) and focuses on the evaluation’s improving and minimization. In semi-flexible docking, the receptor was considered as a rigid structure, while the ligand remained flexible to achieve the best conformation due to the receptor complex. Data on validation protocols for molecular docking can be seen in the Figs 1 and 2. According to the literature data, RMSD values expressing the relationship between the calculated and crystallographic conformations of the ligand’s complex must have a difference no more than 2.0 Å (Hevener et al. 2009; da Silva Costa et al. 2018).

Figure 1. 

Data obtained during validation of molecular docking protocols for the receptor (Peroxiredoxin).

Figure 2. 

Receptors in complexes with calculated and experimental ligand conformation: peroxiredoxin-5 bond with (4S,5S)-1,2-dityan-4,5-diol (A), hemoxygenase 1 with protoporphyrin IX (B), ionotropic glutamate receptor (NMDA 1) with 1-sulfonyl[1,2,4]triazolo[4,3-a]quinoxalin-4(5H)-one (C), NO-synthase with 6-(3-fluoro-5-(3-(methylamino)prop-1-in-1-yl)phenethyl)-4-methylpyridine-2-amine (D).

Similarity in overlapping crystallographic poses (orientation + conformation, blue) and calculated ones (yellow) was obtained by molecular docking and graphically reflects a low RMSD value, which characterizes good results according to the literature data. Active sites are the ligand coordinates in the original target protein networks, and these active target receptor binding sites were analyzed and determined by means of Chimera 1.16. PyMOL v.2.5 (Schrodinger, New York) and Discovery studio Visualizer (BIOVIA 2021) were used to create receptor-ligand figures. The coordinates of the active site and the range of ligand conformations are given in Table 1.

The ideal position of each molecule was selected according to the energy index and the best fit to the active center with the lowest RMSD index. The lower ∆G, the more significant the interaction between the receptor and the ligands.

In the given study, we’ve conducted two molecular dynamic modeling experiments to support our design concept. Two complexes were taken as a basis, where the ligands had the lowest energy conformation, which also exceeded the indicators of crystallographic comparison preparations. The first experiment was for NO-synthase (PDB: 6NGJ) in complex with 2.15, and the second for the NMDA-receptor (PDB: 4KFQ) with compound 3.14. MD simulations were performed in a GROMACS 2202.2 (Abraham et al. 2015) package with CUDA acceleration and CGenFF force field (Vanommeslaeghe et al. 2009) running on ОС Ubuntu WSL 18.04. The topology’s receptors were obtained by means of the pdb2gmx command. Ligands’ parameterization was performed on the server https://www.swissparam.ch (Zoete et al. 2011). All systems were solvated using a water Model (TIP3P) with a dodecahedron box configuration with a distance of 1 nm from the protein edges in all directions. Na⁺ and Cl- ions were added to neutralize the system, and NaCl with a concentration of 0.15 M was added to simulate a physiological buffer solution. To minimize energy and eliminate any steric collisions and unwanted contacts, a 500-step steep descent algorithm was used, in which F_max did not exceed 1000 kJ/mol nm. After that, the NVT and NPT systems were balanced, maintaining a constant pressure on 1 bar level using a Berendsen barostat and a temperature of 300 K by means of the modified Berendsen thermostat for 100 ps. MD modeling was performed during 100 ns for each of the two complexes with stable temperature and pressure with a time step of 2 fs and a long-range interaction limit of 1 nm using Particle-Mesh Ewald electrostatics (Essmann et al. 1995). Trajectories were analyzed by means of various GROMACS scripts and they were visualized using Pymol. Graphs depicting simulation results were built on the basis of the QtGrace tool (QtGrace download SourceForge.net).

The molecular Mechanic / Poisson-Boltzmann Surface Area (MM-PBSA) (Kumari et al. 2014) is a method used to determine the affinity of an inhibitor protein. Typically, this free binding energy of protein and ligand complexes can be calculated using g_mmpbsa together with MD modeling by means of the equation:

ΔG_binding = ΔG_complex – (ΔG_protein + ΔG_ligand) (1)

In addition, the individual free energies for a complex, protein, or ligand can be calculated using the equation (2):

ΔG = (E_MM) + G_solvation (2)

Where (E_MM) is the average potential energy of molecular mechanics without taking into account pressure. The average free solvation energy consists of two parts: polar and nonpolar (3):

ΔG_solvation = ΔG_polar + ΔG_nonpolar (3)

That is, the binding energy consists of three energy terms of the system’s individual components: potential energy in vacuum, polar solvation energy, and nonpolar solvation energy. A full description of the MM/PBSA protocol was obtained from the Web page (http://rashmikumari.github.io/g_mmpbsa/).

Molecular docking was used to identify the mechanism of interaction between ligands and receptors as the primary method of virtual screening. Docking studies were conducted for 112 compounds due to six enzymes responsible for regulating of the oxin process. Nine conformational positions were created for each compound and one with the best affinity value was selected, which is presented in the Table 2. Compounds that exceeded the binding affinity energy values of the crystallographic reference ligands are highlighted in bold.

Twenty-three found hits that exceeded the docking value of the comparison ligands were subjected to ADME analysis using the SwissADME server (Table 3). Several rules have been used to determine the similarity of compounds to drugs; the most popular was Lipinski’s rule of five. The ADME properties included mol. mass, topological polar surface area, number of H-bond donors and acceptors, predicted solubility in water, value of the permanent blood-brain barrier, predicted metabolism and detoxification of xenobiotic compounds due to inhibition of cytochromes P450.

The results of this study showed that almost all compounds obeyed the Lipinski’s rule. Only compound 3.19 exceeded the molecular weight of 500 g/mol and the value of consensus log P > 5. Therefore, the following compound (3.14) was chosen for further virtual screening and protein-ligand interaction analysis based on the value of the energy’s estimate.

The results of docking have shown that among the studied compounds, 6 molecules (2.14, 1.15, 2.15, 1.16, 3.16, 3.19) the binding affinity values (from -9.7 to -10.4 kcal/mol) due to NO- synthase (PDB: 6NGL) relative to the control ligand were exceeded. Compound 2.15 has formed a complex with oxidoreductase and a better docking score of -10.4 kcal/mol. The stability of the complex was characterized by the presence of hydrogen bonds between the ligand and SER413 (3.11 Å), ARG414 (2.92 Å), CYS415 (3.00 Å), TRP409 (3.02 Å). The interaction of hydrogen bonds plays a vital role in maintaining the stability of the attached Complex (Panigrahi 2008). The higher activity of the ligand was probably due to the hydrophobic interaction of the π-π stacking TRP409 (3.99 Å), PHE584 (4.25 Å) and 4-fluorophenyl radical, as well as the hydrogen bonds of the carbonyl group and CYS415. The surface of the protein binding site is mostly surrounded by other bonds: ALA412 (3.81 Å) and CYS415 (3.89 Å) has a hydrophobic π-alkyl contact with the aromatic rings of 1,2,4-triazole, the presence of pyridine-4- yl causes an electrostatic π-cation bond with the AGR414 residue (3.77 Å). Active residues that interact via hydrogen, hydrophobic, and electrostatic bonding are shown in the Fig. 3, and Fig. 3D pocket representations are also provided.

Figure 3. 

2D interaction of the 2.15 compound and NO-synthase and 3D image of the active site’s surface of the enzyme.

The docking’s result of the studied compounds with the NMDA-receptor GluN1has shown that the 1,2,4-triazole derivative 3.19 has the highest binding affinity energy of -10.3 kcal/mol, but the compound violates the Lipinski’s rule, so the next ligand 3.14 with an affinity of -10.1 kcal/mol was chosen. The binding mode is presented in (Kvist et al. 2013) and demonstrates the formation of a stable complex through hydrogen bonds with 6-(3-fluoro-5-(3-(methylamino)prop-1-yn-1-yl)phenethyl)-4 -methylpyridin-2-amine: PRO124, THR126, ARG131 and SER296. The PHE92 residue stabilizes the binding of the reference ligand by means of a π-stacking interaction. In our case, PHE92 (4.97 Å) had a stabilizing π-π stacking with a 3-fluorophenyl substituent, respectively (Fig. 4).

Figure 4. 

Interaction of 2D and 3D graphic contacts of N-methyl-d-aspartate (NMDA) receptor GluN1 and compound 3.14.

The O1 atom of the carbonyl group has formed a weak hydrogen bond with the hydroxyl group of the THR126 residue (2.96 Å), the N1 atom of the 1,2,4-triazole ring is connected with the amino group of ASN128 (3.16 Å). Other hydrophobic, electrostatic, and unconventional hydrogen bonds were inherent in our complex: the electron cloud of the aromatic ring of the phenyl radical formed a π-alkyl contact with amino acid residues VAL181 (5.27 Å), LEU (4.92 Å). The presence of a fluorine atom causes a halogen bond with the carboxyl group of ASP224 (4.96 Å). Also, the π-anion bond is a stabilizing electrostatic interaction of the ASP224 carboxyl group (4.82 Å) with the polarizable Pi- electron cloud of the 3-fluorophenyl’s aromatic ring.

The molecular dynamics method was used to test the stability of Ligand-enzyme complexes. Thus, we could obtain important information about the dynamic behavior of the studied compounds considering the long trajectories, in the 100 ns time function and in the solvated medium. The physiological salt concentration was also maintained. After a cycle of 100 ns for each of the MD complexes, the simulation was analyzed for such parameters as: RMSD, RMSF, Rg, SASA, Hbonds, PCA. Molecular dynamics study was performed for complexes: proteins with the studied compounds that had the highest value of the docking index and exceeded the docking score of crystallographic comparison preparations (6NGJ_2.15; 4KFQ_3.14); proteins with crystallographic reference ligands, which were used to confirm the suitability of the molecular docking protocol, where for NO-synthase a complex of 6-(3-fluoro-5-(3-(methylamino)prop-1-in-1-yl)phenethyl)-4-methylpyridine-2-amine (6NGJ_STD), and for the NMDA GluN1 receptor, modeling was performed with 1-sulphanyl[1,2,4]triazolo[4,3-a]quinoxalin-4(5H)-one (4KFQ_STD). Also, for comparison, molecular modeling of proteins’ native forms (6NGJ_Apo; 4KFQ_Apo) was carried out.

All changes in the protein and protein-ligand complexes configuration were investigated by means of square deviation (RMSD) of the protein backbone skeleton’s coordinates during 100 ns of MD modeling. The RMSD graph shows that the calculated RMSD for the complex of NO-synthase and ligand 2.15 (Fig. 5A) was 0.37 ± 0.05 nm and was 0.05 nm larger than the RMSD of the native protein (0.32 ± 0.02 nm). For the 6NGJ_STD complex, the RMSD was 0.25 ± 0.03 nm. The RMSD graph shows that for the no synthase complex and ligand 2.15 (Fig. 4A) the calculated RMSD was 0.37 ± 0.05 nm and was 0.05 nm greater than the native protein RMSD (0.32 ± 0.02 nm). For the 6NGJ_STD RMSD complex, the RMSD was 0.25 ± 0.03 nm. It can be noted that the crystallographic reference ligand stabilizes the synthase structure to a certain extent. The rms deviation of the atomic positions of the protein framework and the complexes reached equilibrium in the initial ∼10 ns of simulation and produced stable trajectories. Thus, compound 2.15 in interaction with the protein, when a plateau appears on the RMSD graph, can take on new poses that differ from those ones indicated in the docking simulation. Therefore, when ligand 2.15 is introduced into the conformational sphere NO-synthase it leads to insignificant destabilization of the complex. In the case of NMDA receptor apo, the form had the highest RMSD value, which was 0.45 ± 0.17 nm (Fig. 5B). For coupled systems, the stability was higher, and throughout the simulation, the average RMSD of 4KFQ_STD was 0.37 ± 0.15 nm, and 4KFQ_3.14–0.31 ± 0.05 nm, which indicates stable conformational behavior during dynamics. From the above mentioned, it can be concluded that the studied compound 3.14 behaves better in the active center of the 4KFQ protein than the reference complex.

Figure 5. 

RMSD diagram: A. apo forms of NO-synthase and selected 6NGJ_2.15 and 6NGJ_STD complexes; B. native NMDA receptor, 4KFQ_3.14 and 4KFQ_STD complexes.

The RMSD analysis shows that the MD trajectories of the two studied complexes were generally stable throughout the simulation period.

The average position fluctuation is calculated from the standard deviation (RMSF) values of each residue to test the mobility or flexibility of protein residues during modeling. This allowed us to note that for NO-synthase and the corresponding complexes, the position of amino residues do not differ much. They showed a slight fluctuation difference between the native protein and the complexes, i.e. they maintained the same position of the main scaffold of the residues. The average RMSF value for the Apo form was 0.10 ± 0.07 nm, the RMSF of 6NGJ due to the compound 2.15 was 0.14 ± 0.08 nm, and for the reference ligand 0.15 ± 0.09 nm (Fig. 6A). In the case of synthase and bound ligands, the RMSF graph’s analysis shows a single peak above 0.3 nm, which indicates a large variation of amino acids in the range of Cys326 – Agr343, which are located outside the ligand binding interface. The RMSF value for the native NMDA receptor was equal to 0.26 ± 0.13 nm, in relation to the ligand 3.14–0.16 ± 0.11 nm, and for the protein complex and the reference compound 0.25 ± 0.11 nm (Fig. 6B). The highest fluctuation values were characteristic of the amino acid residues that did not belong to the active site of the enzyme and were observed in two ranges – Asn48-His57 and Arg97-Lys104. The reduction of the average RMSF by almost half in 4KFQ_3.14 shows that the secondary conformations of the receptor remain more stable during the 100 ns simulation.

Figure 6. 

Final RMSF deviations (nm): A. RMSF of atoms of the 6NGJ protein backbone, 6NGJ_2.15 and 6NGJ_STD complexes; B. RMSF of atoms of the 4KFQ protein skeleton, 4KFQ_3.14 and 4KFQ_STD complexes.

The rotation radius indicates the stability of biomolecules by measuring their structural compactness along the molecular dynamics trajectory. The Rg value is relatively constant and have a stable curve throughout the simulation for the two 6NGJ_Apo complexes, 6NGJ_2.15 and 6NGJ_STD, and they are in average 2.29 ± 0.013 nm, 2.33 ± 0.017 nm, and 2.31 ± 0.019 nm, respectively (Fig. 7A). Comparing the Rg complexes and apo forms of the corresponding protein, it can be argued that the complexes showed relatively similar compactness behavior, and no degradation occurred. In conclusion, the higher Rg for 6NGJ_2.15 shows a more detailed structure of the protein-ligand system compared to the native protein and the protein with the reference compound. The graph shown in Fig. 7B demonstrates the average Rg apo values of the 4KFQ form, 4KFQ_3.14, and 4KFQ_STD complexes, and they were 2.11 ± 0.08 nm, 2.06 ± 0.03 nm, and 2.08 ± 0.06 nm, respectively. A lower Rg for 4KFQ_3.14 shows greater compactness of the protein-ligand complex. As a result, each complex showed relatively similar behavior of compactness and a constant amount of Rg compared to the native protein and the reference one.

Figure 7. 

Rotation radius of protein skeleton atoms: A. Apo-protein 6NGJ, of 6NGJ_2.15 and 6NGJ_STD complexes; B. of native protein 4KFQ, of 4KFQ_3.14 and 4KFQ_STD complexes.

Calculating the surface accessible solvent area (SASA) helps you to calculate the surface area of the protein and complex that is available for the solvent. The average SASA value is calculated for the native protein is 225.66 nm ± 1.50 nm, while the protein complexes with the studied compound 2.15 and the reference ligand had values of 235.87 nm ± 1.67 nm and 237.24 ± 1.52 nm, which indicates a slight difference between the studied systems (Fig. 8A). Fig. 8 B shows a graph of the SASA value’s dependence for the NMDA receptor and its associated complexes, in which the average values were very similar to 4KFQ Apo, 4KFQ_3.14, and 4KFQ_STD – 160.01 ± 1.23 nm, 160.28 ± 1.36 nm, and 160.23 ± 1.22 nm, respectively. These calculations show that all two studied complexes have a significantly similar SASA value to the reference complex.

Figure 8. 

Solvent-accessible surface area: A. 6NGJ_Aro, 6NGJ_2.15 and 6NGJ_STD complexes; B. 4KFQ_Apo, 4KFQ_3.14 and 4KFQ_STD complexes.

Hydrogen bonds are a key indicator of binding specificity between a receptor and a ligand. The average values of H-bonds after the simulation period of 100 ns formed between NO-synthase and molecule 2.15 were 2.40 ± 1.63 (Fig. 9A), and between the NMDA receptor enzyme and 3.14 molecule were 0.76 ± 0.72 (Fig. 9B).

Figure 9. 

proteins H-bonds with the studied ligands: A. 6NGJ with compound 2.15; B. 4KFQ with compound 3.14.

Principal Component Analysis (PCA) is a statistical calculation for reducing the amount of data, which in our case was based on extracting significant movements of backbone Cα atoms due to the ligand. The first 40 eigenvectors and eigenvalues of the two studied enzymes with and without ligands were selected for this study. The first 40 eigenvectors were 99.89% for 6NGJ_2.15, 99.42% for 6NGJ_STD, 99.17% for 4KFQ_3.14, and 99.99% for 4KFQ_STD. The flexibility of all systems was analyzed by calculating the trace value of the diagonalized covariance matrix, which is the sum of the Eigen values. For the 6NGJ_2.15 complex, this value is 4.08 nm², while 6NGJ_STD is 1.08 nm² (Fig. 10A). It can be argued that the complex has increased flexibility with respect to movement compared to the comparison complex, probably due to increased local movements for the complex, which was also shown by the Rg graph. On a 2D projection (Fig. 11A) it is shown that the protein structure has 5 clusters, but they occupy a smaller phase space compared to a complex that has 4 clusters. The trajectory projection for 4KFQ_3.14 and 4KFQ_STD, which is shown in Fig. 10B, demonstrates the sum of the trace values of the covariance matrix and is 1.52 nm² and 12.46 nm², respectively. The receptor complex with the test ligand 3.14 was very stable, because they occupied less space in the phase space, and the cluster was well expressed compared to the protein-reference ligand complex (Fig. 11B).

Figure 10. 

Analysis of principal components. Plots of eigenvalues with respect to the first 40 eigenvectors of the complex: A. synthase with compound 2.15 and reference ligand; B. 4KFQ_3.14 and 4KFQ_STD.

Figure 11. 

PCA scatterplots show the displacement projections of Cα atoms along the first eigenvector and the second eigenvector at each time point for: A. 6NGJ_2.15 and 6NGJ_STD; B. 4KFQ_3.14 and 4KFQ_STD.

Summarizing all the above mentioned, we can conclude that 2.15 and 3.14 compounds for which molecular dynamic characteristics were studied confirmed the docking study, since these tested molecules tended to remain bound to the selected enzymes and did not fall out of the active site throughout the simulation. It should be noted that the synthase complex with ligand 2.15, was although quite stable, but was unfortunately inferior in stability to the reference complex. The 6NGJ_STD complex occupied a smaller confirmation space, had fewer residual fluctuations, and even showed less flexibility of the protein skeleton compared to Apo synthase. In contrast, the NMDA receptor complex with 3.14 has showed a fairly significant difference from the reference complex and Apo protein. The studied ligand caused a greater compactness of the enzyme structure, the complex itself had fewer residual fluctuations, therefore it was more stable than the reference 4KFQ_STD complex.

The free binding energy of the simulated complexes was estimated by the MM-PBSA method for the last 20 ns trajectories (Table 4). The calculated ΔG value for the 6NGJ_2.15 and 4KFQ_3.14 complexes was 5.33 ± 2.79 and -24.14 ± 0.32 kcal/mol, respectively.

The binding energy of the complex showed a small difference between the studied ligand 2.15 and the reference ligand. Although the free binding energy was higher in the studied complex, similar stability of the systems was still present. Instead of, the 4KFQ complex shows better binding affinity to compound 3.14 (-24.14 ± 0.32 kcal/mol) than to the reference ligand 1.24 ± 2.77 kcal/mol. The electrostatic energy shows significant moderate values in the case of the 4KFQ_3.14 complex. The binding energy showed the stability of the complex, and we can assume that ligand 3.14 can be used as a potential inhibitor/antagonist for the selected enzyme.