An in silico investigation of 1,2,4-triazole derivatives as potential antioxidant agents using molecular docking, MD simulations, MM-PBSA free energy calculations and ADME predictions

) An in silico investigation of 1,2


Introduction
Oxidative stress is an unavoidable metabolic process for all living organisms, in which there is an imbalance between the Reactive Oxygen Species (ROS) and Reactive Nitrogen Species (RNS) and the mechanisms of antioxidant protection. Reactive species of N and O are formed during chain reac-Tyrosinase is an enzyme that regulates the rate of melanin production, and its abnormal accumulation causes various skin diseases. Therefore, tyrosinase inhibitors are necessary to ensure normal melanin content (Pintus et al. 2017). Hemoxygenase-1 (НО-1) has a specific ability to regulate oxidative processes by forming antioxidant bilirubin during heme catabolism, and then converting to biliverdin when reacting with ROS (Kapitulnik 2004). ROS and RNS are usually generated by enzymes such as NO-synthase and NAD(P)H-oxidase, peroxidase. NA-D(P)H excessive stimulation causes the production of reactive oxygen species, leading to oxidative stress (Tarafdar and Pula 2018).
1,2,4-triazoles and their derivatives represent a huge class of heterocyclic compounds that have a wide range of pharmacological activity. They can easily bind to biotargets, since they demonstrate good stability in various media, while increasing the ability to form hydrogen bonds, Dipole-Dipole and π-stacking interactions. The presence of specific functional capabilities in combination with practical non-toxicity (Karpenko et al. 2022;Safonov et al. 2022) contributes to find in them such biological properties as: antioxidant , antitumor (Jiang et al. 2020), antimicrobial (Strzelecka and Świątek 2021), anti-inflammatory (Azim et al. 2021) and antituberculosis activity (Gotsulya et al. 2020).
Selected enzymes capable of regulating oxidative processes have already been the subject of research. For example, works (Yapati et al. 2016;Sudhamani et al. 2019) describe new thiourea derivatives of 5-hydroxytryptophan and metal complexes of benzo [d]thiazol that are able to bind to the most important amino acid residues in the active center of 3MNG -Asn24, Val94, Gly46, Leu96. Recently, carvacrol-derived sulfonamides (de Oliveira et al. 2021) have been investigated that interacted with the crystallographic structure of NO synthase (6NGJ) and the NMDA-glycine binding site (4KFQ). In the first case, the sulfonamide derivatives formed hydrogen bonds with Val677, Arg596, Ser334 and π-interaction with Trp678, and in the second case, the compounds were connected by hydrogen bonds with Gln405, Arg523, Thr518, and Ser572. The 1,2,4-trazoles themselves have shown that they are able to act as adenosine A2B receptor agonists (Turky et al. 2020), and to form hydrogen bonds with receptors that are a part of tumors (Aliabadi et al. 2016). The article (Akın et al. 2019) highlights the ability of the 4-(substituebenzyl)-2-heptyl-5-methyl-2,4,-dihydro-3H1,2,4-triazole-3-one derivative to bind to the active center of the enzyme tyrosinase (2Y9X), interacting with His244, His263, Phe264 and Val283.
In this work, our goal was to perform virtual screening of 112 previously synthesized compounds Polishchuk 2020, 2021;Fedotov and Hotsulia 2021;Karpun 2021;Khilkovets 2021a, b) with different functional groups and pharmacophores using a molecular docking strategy with 6 proteins responsible for regulating antioxidant action and subsequent prediction of absorption, distribution, metab-olism, and excretion (ADME) properties. Molecular dynamic modeling was used to evaluate the trajectories and interactions of the receptor-ligand complex with the highest affinity values. To get a better understanding of the structure-activity relationship of the studied compounds, the free binding energy (MM-PBSA) was also calculated.

Target selection and preparation
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).

ADME properties analysis
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. Scheme 1. Scheme of the studied derivatives of 1,2,4-triazoles used in computable investigations.

Evaluation of molecular docking. Docking protocol
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).
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.

Molecular dynamics (MD) simulation of top docked complexes
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 Clions 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).

MMPBSA analysis
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: In addition, the individual free energies for a complex, protein, or ligand can be calculated using the equation (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): 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/).

Docking study
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.

Determination of ADME for lead compounds
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.

Molecular dynamics (MD) simulation study
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-1yl)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, mo-lecular 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.
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.
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.
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.
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).
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 2 , while 6NGJ_ STD is 1.08 nm 2 (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 2 and 12.46 nm 2 , 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).
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.

Free binding energy estimation
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 mod-  erate 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.

Conclusions
In this study, we've used the in silico approach to identify bioactive compounds that can inhibit six proteins responsible for antioxidant regulation. The number of selected 1,2,4-trazole derivatives in 112 molecules was sufficient to obtain promising results in terms of the potential to exhibit inhibitory activity. Six molecules have shown the best Affinity value against NO-synthase (from -9.7 to -10.4 kcal/mol), exceeding the ligand reference values. And twenty-three hits had a higher docking score (from -8.4 to -10.3 kcal/mol) to the N-methyl-d-aspartate (NMDA) GluN1 receptor than the reference ligand. The ADME prediction of priority ligands has revealed high bioavailability for the analyzed molecules and only one compound 3.19, which had a violation of the Lipinski's rule, so the next compound 3.14 was investigated instead. Molecular modeling study of the two systems with the best affinity indicators has revealed the stable complexes formation throughout the simulation. In particular, analysis based on molecular dynamic characteristics (RMSD, RMSF, Rg, SASA and PCA) has shown that the interaction of the ligand 3.14 with 4KFQ was much more stabilizing and conformationaly favorable compared to the reference. Calculations of the binding free energy of the complex have shown similar values of the studied ligand with oxidoreductase, but exceeded them compared to the reference ligand. In contrast, ligand 3.14 had better results in terms of MM-PBSA binding energy, suggesting that this molecule is likely to be a good hit in the discovery of competitive inhibitors of the NMDA receptor GluN1. Thus, our theoretical results indicate that the studied molecules can be used as potential candidates for further in vitro and in vivo studies on antioxidant activity.