Imidazole and nitroimidazole derivatives as NADH-fumarate reductase inhibitors: Density functional theory studies, homology modeling, and molecular docking.

Chagas disease is caused by Trypanosoma cruzi. Benznidazole and nifurtimox are drugs used for its therapy; nevertheless, they have collateral effects. NADH-fumarate (FUM) reductase is a potential pharmacological target since it is essential for survival of parasite and is not found in humans. The objectives are to design and characterize the electronic structure of imidazole and nitroimidazole derivatives at DFT-M06-2X level in aqueous solution; also, to model the NADH-FUM reductase and analyze its intermolecular interactions by molecular docking. Quantum-chemical descriptors allowed to select the molecules with the best physicochemical properties and lowest toxicity. A high-quality three-dimensional structure of NADH-FUM reductase was obtained by homology modeling. Water molecules do not have influence in the interaction between FUM and NADH-FUM reductase. The main hydrogen-binding interactions for FUM were identified in NADH, Lys172, and Arg89; while hydrophobic interactions in Phe479, Thr174, Met63. The molecules S3-8, S2-8, and S1-8 could be inhibitors of NADH-FUM reductase.

INTRODUCTION

American Trypanosomiasis, also called Chagas disease, is a parasitic disease caused by the hemoflagellate parasite Trypanosoma cruzi. Currently, the World Health Organization (WHO) classifies it as one of the most important neglected tropical diseases worldwide. It is estimated that approximately 6–7 million people are infected and less than 1% receive treatment, , , , causing around 14,000 deaths annually in endemic regions. At present, Chagas disease has an incidence of 30,000 new cases per year in the Americas, with an economic impact in 2015 of almost USD 450 million in Latin America. This disease is transmitted mainly by the bite of triatomines, but can also be transmitted congenitally, by transplants, blood transfusions, among other forms. , The only drugs approved by the WHO for its treatment are two nitroaromatic compounds: nifurtimox and benznidazole (BNZ), which are effective if they are administered at the onset of the infection; however, it has been reported that their effectiveness decreases and their toxicity increases in long‐term treatments, which makes them banned in many countries. , , , , Such toxicity relates to the formation of free radicals generated by the reduction of the nitro group, which puts the parasite under oxidative stress and causes its death, but also causes damage to human cells. Nevertheless, in the literature it has been reported that the antichagasic effect of BNZ does not depend on the production of free radicals and can act against the parasite in different ways, so its main mechanism of action is still unknown. , , ,

Turrens et al., analyzed the trypanocidal activity of BNZ and other nitroimidazoles, as well as of some anthelmintic benzimidazoles. The results of the study showed that all the molecules studied decreased the parasite proliferation of T. cruzi and Trypanosoma brucei, and inhibited the reducing action of the enzyme NADH‐fumarate (FUM) reductase in both parasites; however, its activity decreased when the imidazole ring was replaced by a thiazole ring. These findings led the researchers to conclude that the imidazole ring was responsible of the trypanocidal activity. Furthermore, this analysis allowed the enzyme NADH‐FUM reductase to be identified for the first time as a potential and specific pharmacological target of BNZ since it was not found in humans or mammalian hosts.

The NADH‐FUM reductase is an enzyme present in the mitochondria of all morphological stages of trypanosomatid parasites (such as T. cruzi, Leishmania sp., and T. brucei), because of this, it is considered one of the most important enzymes in their metabolism. Likewise, it catalyzes the reduction of FUM to succinate, which is essential for them to carry out their electron transport chain, maintain redox balance and obtain the necessary energy for their survival. , ,

Previous studies of the electronic structure, at the ab initio level and at the density functional theory (DFT) level, allowed us to establish a theoretical model and a pharmacophore group, , which were used to design a family of new antichagasic molecules derived from nitroimidazoles. These molecules were experimentally analyzed to evaluate their biological activity and cytotoxicity. The results showed that the structures proposed by the theoretical calculations (M3 and ISO3) inhibited parasitic growth in vitro, in the CL‐Brener and CL‐14 strains of T. cruzi and did not present cytotoxicity in macrophages.

Given that the American trypanosomiasis is an important public health problem and in the absence of drugs to treat it, it is necessary to find new therapeutic strategies based on specific pharmacological targets, which leads to a rational design of drugs with better effectiveness and lower toxicity. The objectives of this study were: (a) to design and characterize the electronic structure, the physicochemical properties and the chemical reactivity of imidazole and nitroimidazole derivatives at the DFT‐M06‐2X level in an aqueous solution; (b) to model the three‐dimensional structure of the enzyme NADH‐FUM reductase, using homology modeling; and (c) to visualize and analyze their intermolecular interactions and binding energy via molecular docking studies, with the purpose of identifying new molecules with potential antichagasic therapeutic effect and a reduced chemical toxicity.

COMPUTATIONAL DETAILS

Pharmacokinetic parameters and toxicological prediction

The initial design of new molecules derived from imidazole and nitroimidazole was carried out, considering the theoretical model established in previous studies, , and molecules reported as active against Chagas disease were added as reference and head of series, , , , , , , as shown in Figure 1. The mutagenic, tumorigenic, irritability, and reproductive effects of these structures were analyzed using the virtual database of Osiris Molecular. Similarly, the prediction of its pharmacokinetic properties was carried out employing the Molinspiration database and taking into account Lipinski's Rule.

FIGURE 1

Head structures of series ,

Electronic structure calculations

According to pharmacokinetic parameters and toxicological predictions results, a group of 31 imidazole and nitroimidazole derivatives was selected (see Figures 2 and 3), with the objective to characterize the physicochemical properties and the chemical reactivity.

FIGURE 2

Selected structures of the 2‐nitroimidazole and 5‐nitroimidazole derivatives, according to their pharmacokinetic parameters and toxicological predictions results. For SESE calculations: X is a transmitting moiety (imidazole ring linked to R1 or R2 groups, respectively) and the R9 is the substituent group

FIGURE 3

Selected structures of the imidazole derivatives according to their pharmacokinetic parameters and toxicological predictions results

Electronic structure calculations were carried out using the Gaussian 09 program suite. All neutral structures were optimized employing M06‐2X functional , with a 6‐311+G(d,p) basis set, while protonated (cationic) structures were optimized at restricted RM06‐2X levels, all these with a 6‐311+G(d,p) basis set. Energies were corrected to include zero‐point vibrational energy at the same levels of theory. Single‐point calculations were performed on optimized structures at the M06‐2X level of theory with a 6‐311++G(2df,2p) basis set. All calculations were carried out in an aqueous solution, and the solvent effects were described through the SMD model. First, the proton affinity (PA) energies were calculated by isodesmic reactions, according to previous works, , , with the purpose to know the basic character of the N3 atom of the imidazole ring.

Then, to describe a more reliable chemical reactivity, we determined some selected global quantum chemical descriptors, such as chemical hardness (η), ionization potential (I), and electrophilicity index (ω). The electrophilicity index is defined as the ability of a molecule to accept an electronic charge and was calculated through the electronic chemical potential and the chemical hardness using the following equations. , , , ,

Parr and Pearson , proposed the chemical hardness (η) within the DFT as follows:

Using finite‐difference approximation, Equations (1) and (2) would be: and where , , and , are the energies of N, (N−1), and (N+1) electron systems; I and A are the adiabatic ionization potential and the electron affinity, respectively. The electrophilicity index of a system, in terms of its chemical potential and chemical hardness, is given by the next expression , :

Furthermore, to identify a possible relationship between carcinogenicity and aromaticity of imidazole and nitroimidazole derivatives, we estimated the energetic difference between the molecular orbitals HOMO and HOMO‐1, taking the aromaticity index (Δ) as an indirect measure of the carcinogenic compounds according to Barone's rules. , Taking into account the first of Barone's rules: a molecule will be strongly carcinogenic if it contains a pyrene‐like structure and the parameter Δ is greater than 0.25β (β ≈ 2.4 eV), or if the molecule is nonpyrene and the parameter Δ is greater than 0.15β. Thus, this descriptor could be an indicator that allows the analysis of the relative aromaticity of the imidazole ring and its derivatives.

In the design of novel therapeutic molecules, it is critical to know how a molecule's lateral groups influence each other in the electronic structure, chemical reactivity, and physicochemical properties. Therefore, the homodesmotic reactions , , were calculated using the substituent effect stabilization energy (SESE) descriptor, according to the following equations and Figure 2, and where X is a transmitting moiety (imidazole ring linked to R1 or R2 groups, respectively) and the R9 is the substituent group (see Figure 2). A value of SESE < 0 descriptor indicates larger stabilization energy caused by the substituent effect, and a value of SESE > 0 indicates an increase of pi‐electron delocalization. , ,

In addition, atomic charges fitted to the electrostatic potential , were obtained to analyze the nucleophilic/electrophilic nature sites and investigate the electron‐withdrawing/electron‐donating influence of the substituents, mainly of the substituent ‐R9 on the imidazole ring.

In addition, an interesting local density functional descriptor of reactivity is the Fukui function used to analyze the chemical reactivity and sites selectivity within a molecule. , , , The Fukui function values were obtained from the Hirshfeld's population scheme and are defined by the expressions described below.

There are three types of Fukui functions that are defined depending on the electron transfer. ,

The condensed Fukui functions were calculated using a finite difference approach just as proposed by Yang and Mortimer, according to the following expressions where q , is the electronic population of atom k in a molecule.

Furthermore, we analyzed the frontier molecular orbitals, HOMO and LUMO (The highest occupied molecular orbital and lowest unoccupied molecular orbital, respectively) to investigate the substituent effect on the nucleophilic/electrophilic nature sites of the imidazole and nitroimidazole derivatives.

Template selection and sequence structural alignment

Comparative modeling is based on identifying template structures with the highest sequence similarity to carry out the construction of a three‐dimensional enzymatic model. Since, the crystallized structure of the enzyme NADH‐FUM reductase of T. cruzi is not available, the amino acids sequence was retrieved from the NCBI database (ID: XP_807320.1 and ID: PWV10286.1 ) and a search BLAST‐P was carried out against the Protein Data Bank, in order to identify homologous structures. Four template structures were identified that showed high sequence identity (PDB ID: 5GLG, PDB ID: 1KSU, PDB ID: 1QJD, and PDB ID: 1QO8 ). The sequence of the enzyme NADH‐FUM reductase of T. cruzi was aligned to each selected template using MultAlin, maintaining the default parameters. The alignments were inspected and manually adjusted to reduce gaps and insertions, and ensure that the catalytic sites were conserved.

Model building

The homology models of the enzyme were built using the MODELLER 9.21 software based on its alignments with templates and satisfying spatial restraints. The four templates with the highest sequence identity were used to perform homology modeling and incorporate flexibility in the receptor since the protein shows three different active conformations and one inactive conformation (PDB ID: 1QO8). Ten thousand models were modeled for each of the four templates. The best models were selected based on the DOPE statistical score (Discrete Optimized Protein Energy) and based on the root mean square deviation (RMSD) values against their respective templates. The DOPE statistical score was calculated by MODELLER 9.21 and the RMSD was calculated using ProFit. , Then, we evaluated their stereochemical quality using PROCHECK and selected the best models, according to the Ramachandran plots.

Model optimization

The selected models were refined using the ModRefiner software to generate a better quality structure. ModRefiner is an algorithm that performs a conformational search in the side chains and main chain of each amino acids of the protein through energy minimization. In the same way, we verified the stereochemical quality of each model through Ramachandran plots, using PROCHECK.

Receptor and ligands preparation

The optimized geometries of the ligands, the coenzymes NADH and FAD, and the reference molecules used for the molecular docking studies were taken from the theoretical calculations, at the M06‐2X/6‐311+G(d,p) level of theory in aqueous solution. The ligands preparation process was performed for the molecular docking study by assigning charges, protonation states, and flexibility using the OpenBabel software. For the process of preparing the enzymatic models: 5299 (1KSU), 2735 (1QJD), 4560 (5GLG), and 9109 (1QO8), were assigned protonation states and missing hydrogens, using PROPKA3. , In addition, the rotamers and the positions of the amino acids GLU, HIS, and ASP were optimized using the MolProbity server.

Docking studies

The molecular docking was performed employing Autodock Vina 1.1.2 software. The position of the crystallized FAD of the template structures was taken as a reference to locate the NADH's position and performed the first four couplings to form the NADH‐FRD complexes. The complexes were formed between the NADH and the model 5299 (from the PDB ID: 1KSU template), the NADH and the model 2735 (from the PDB ID: 1QJD template), the NADH and the model 4560 (from the PDB ID: 5GLG template); and the NADH and the model 9109 (from the PDB ID: 1QO8 template).

For each model, a conformational docking space search was set up. The defined search space included the entire catalytic domain, which ensured the analysis of all potential binding sites. For the model 5299 (1KSU), the search space was defined to cover a three‐dimensional analysis region of 36 Å × 21 Å × 18 Å points, centered at x:14.083, y:−2.972, z:5.917. A similar procedure was performed for all models. For the model 2735 (1QJD) the region of 22.50 Å × 21.00 Å × 35.25 Å, centered at x:26.417, y:32.861, z:64.361 was defined. For the model 9109 (1QO8) the region of 32 Å × 22 Å × 22 Å, centered at x:39.417, y:68.472, z:20.500 was defined. For the model 4560 (5GLG) the region of 24 Å × 26 Å × 28Å, centered at x:8.417, y:28.472, z:0.667 was defined.

To achieve the molecular docking between the NADH‐FRD complex and the ligands, the model 5299 (1KSU) were selected, since it is the only model that preserved the FAD‐FRD interaction site of its original template. For these couplings, a grid box with dimensions of 20 Å × 22 Å × 22 Å, centered at x:6.500, y:5.111, z:10.556 was used. In order to analyze the input and output files, the MGLTOOLS 1.5.6 was employed. , The best ligands using coupling energy and dimensional location were selected.

Validation of docking protocol (redocking)

To verify and ensure that the position of NADH is adequate, the validation of the docking protocol was carried out, where the spatial location of the crystallized FAD‐1KSU complex was compared, with the docking of FAD and the 1KSU template used for modeling. This comparison was performed based on the RMSD values obtained. In the same way, docking of the coenzyme NADH with the template 1KSU and the model 5299 was performed, and the RMSD values were compared. In all cases, the Autodock Vina 1.1.2 software and a grid with dimensions of 36x21x18 points, with a spacing of 1 Å, centered at x: 14.083, y: −2.972, z: 5.917 was used.

Evaluation of the role of water molecules in the binding pocket of NADH‐fumarate reductase

Rational drug design studies have reported that water mediates some ligand–protein interactions. , , , , In this sense, it is postulated that a ligand must first displace the water molecules from the area that defines its binding pocket to interact with the protein. However, some water molecules—also known as nondisplaceable or conserved—interact strongly and remain attached to the protein, thus affecting the positioning of ligands within the binding site and, consequently, influencing the ligand‐binding affinity. , , , ,

There is evidence that imidazole derivatives inhibit the enzyme NADH‐FUM reductase and that these ligands may bind at the FUM–binding pocket. Therefore, before docking imidazole and nitroimidazole derivatives into the homology model of the enzyme NADH‐FUM reductase, we investigated the potential presence of conserved (nondisplaceable) water molecules within the binding site, which could play a role in the mode of binding and affinity of ligands. To Identify such water molecules, WaterDock and WaterDock2.0 were used. These programs dock water molecules within the protein's binding pocket to identify those that are nondisplaceable or conserved.

WaterDock predicts the location conserved water molecules within the binding pocket without the explicit presence of the ligand. This program is relevant because it makes the prediction independently of the nature of the ligand. By contrast, WaterDock2.0 predicts the location of conserved water molecules with the explicit presence of the ligand and includes solvation data of the ligand functional groups. We used both programs to obtain a more complete and accurate estimation of the potential presence of conserved water molecules within the ligand‐binding site. In addition, these programs were selected because they were validated using crystallographic data, showing an extremely high‐success rates: 88% accuracy for WaterDock and 91% for WaterDock2.0.

RESULTS AND DISCUSSION

In the present discussion, the BNZ is the reference molecule, since it is used as a first‐line drug to treat Chagas disease. From BNZ, we made structural modifications with the different functional groups reported as active against T. cruzi, as previously mentioned. The structures analyzed in this study comprise the following series: first, 2‐nitroimidazoles with R2 substituent in position five. Second, 5‐nitroimidazoles with R1 substituent in position two. Third, nonnitrated imidazole derivatives with R1 substituent in position two. Fourth, nonnitrated imidazole derivatives with R2 substituent in position five, as shown in Figures 2 and 3. According to the results of in silico prediction and considering Lipinski's Rules, we selected 32 molecules that were not toxic. It is important to notice that the results of the in silico predictions showed that the BNZ and other reference structures indicated a potential reproductive risk. Similarly, metronidazole (MTZ) showed negative values in the partition coefficient (Log P), which could result in poor absorption of the drug (Tables S1 and S2, See Supporting Information).

Calculations of the electronic structure

The optimized geometry of the selected molecules for imidazole and nitroimidazole derivatives is shown in Figures 4 and 5. The geometrical parameters agree with the available X‐ray experimental data. ,

FIGURE 4

Optimized structures for 2‐nitroimidazole and 5‐nitroimidazole derivatives at the M06‐2X/6‐311+G(d,p) levels of theory. The H, C, N, S, and O atoms are represented as white, gray, blue, yellow, and red balls, respectively

FIGURE 5

Optimized structures for imidazole derivatives at the M06‐2X/6‐311+G(d,p) levels of theory. The H, C, N, S, and O atoms are represented as white, gray, blue, yellow, and red balls, respectively

On the other hand, the analysis of electron density‐based molecular descriptors in pharmacological molecules is important. In this sense, it is essential to know the interaction sites and the chemical reactivity of these molecules. Tables 1 and S3 show the calculated values for the global reactivity quantum chemical descriptors: PA energy, chemical hardness (η), ionization potential (I), electrophilicity index (ω), the SESE, and aromatic index (Δ).

TABLE 1

Proton affinity and quantum‐chemical descriptors for reference and selected molecules, at the M06‐2X/6‐311++G(2df,2p)//6‐311+G(d,p) levels of theory

	Molecule	PA (kcal/mol)	η (eV)	I (eV)	ω (eV)	SESE (kcal/mol)	∆ (eV)
Ref. Mol.	BNZ*	15.37	10.18	6.83	0.15	−0.47	0.24
	M3**	8.21	10.09	6.88	0.17	10.23	0.18
	M35m	14.67	10.18	6.88	0.16	−0.70	0.08
	M2**	12.92	10.01	6.82	0.16	0.42	0.16
	M7D	8.92	10.46	7.20	0.19	7.66	0.20
Imidazoles	S1‐1	−0.15	7.36	6.47	0.53	–	0.31
	S1‐6	0.24	7.69	6.54	0.47	–	0.07
	S1‐8	0.05	8.39	6.45	0.30	–	0.21
	S2‐1	−2.36	7.18	6.22	0.48	–	0.68
	S2‐6	−2.40	7.39	6.24	0.44	–	0.15
	S2‐8	−2.29	8.16	6.22	0.28	–	0.08
	S3‐1	−1.91	7.15	6.19	0.48	–	0.71
	S3‐6	−0.58	7.39	6.21	0.43	–	0.21
	S3‐8	−1.30	8.12	6.19	0.28	–	0.10
	S4‐1	1.80	6.96	6.02	0.46	–	1.01
	S4‐8	1.80	7.95	6.02	0.26	–	0.30
Nitroimidazoles	2NS5‐2	16.86	10.78	7.38	0.18	−2.02	0.08
	2NS1‐6	16.60	9.89	6.55	0.13	2.40	0.11
	2NS5‐3	16.69	10.40	7.02	0.16	−2.03	0.32
	2NS3‐6	16.11	9.85	6.54	0.13	0.88	0.11
	2NS3‐7	14.94	10.51	7.22	0.18	−1.44	0.15
	5NS2‐6	8.47	9.75	6.53	0.14	10.44	0.07
	5NS1‐1	11.17	10.17	6.91	0.16	11.20	0.16
	5NS5‐4	14.18	10.23	6.97	0.17	−0.27	0.25

Abbreviations: eV, electron Volts; I, ionization potential; PA, proton affinity energy; SESE, substituent effect stabilization energy; Δ, aromatic index; η, chemical hardness; ω, electrophilicity index.

Values obtained from Reference [16].

Values obtained from Reference [15].

The PA values of the reference molecules and their nitrated and nonnitrated analogs are shown in Tables 1 and S3; a larger value of the protonation energy implies higher basicity. Hence, we analyzed the basicity features of the nitrogen atom (N3) of the imidazole group with respect to the R9‐substituent groups and their relationship with the position of the nitro group. The results indicated that the basicity of the N3 atoms is larger for 2‐nitroimidazole derivatives, 2NS5‐2, 2NS5‐3, 2NS1‐6, 2NS3‐6, BNZ, 2NS3‐7, and M35m than the 5‐nitroimidazoles derivatives, 5NS5‐4, M2, 5NS1‐1, M7D, 5NS2‐6, and M3. This behavior could be due to the proximity of the nitro group to the N3 atom and the influence exerted by their inductive and resonance effects. In addition, the imidazole derivatives without the nitro group have the lowest PA values, indicating that the presence of the nitro group is crucial for preserving the basicity of the N3 atom. These results are consistent with previously reported studies.

Chemical hardness (η) is a descriptor of global reactivity, in which larger energy values are related to more stable molecules. In Tables 1 and S3, it is observed that the nitrated molecules have larger chemical hardness values compared to their respective nonnitrated imidazole analogs, which suggests that these molecules are less reactive. Moreover, it is observed that the nonnitrated molecules that possess a phthalazine group in the R9 substituent (S1‐8, S2‐8, S3‐8 and S4‐8) have higher values of chemical hardness than its analogs with benzimidazole substituent (S1‐6, S2‐6, S3‐6) or benzyl groups (S1‐1, S2‐1, S3‐1, S4‐1) as shown in Tables 1, S3 and Figures 3 and 5. These results indicate that the presence of the phthalazine group in the R9 substituent confers more stability to nonnitrated imidazole derivatives. In addition, a similar behavior was observed for phthalazine and benzimidazole substituents in the 5‐nitroimidazole derivatives (Figures 2 and 4).

The ionization potential (I) values are shown in Tables 1 and S3. It is observed that the 2‐nitroimidazole and 5‐nitroimidazoles derivatives, as well as molecules of series 1 (S1‐1, S1‐6, and S1‐8), have higher ionization potential values than nonnitrated imidazole derivatives. These results suggest a lower oxidative effect on nonnitrated imidazole derivatives with R1= H and nitroimidazoles with electron‐donor substituents in R9, indicating an increase in its structural stability.

Furthermore, to know how the substituent groups influence these molecules' stability and chemical reactivity, we investigated and evaluated the effect of the R9 substituents and the nitro group on the nitroimidazole derivatives via the (SESE) descriptor, where more negative values indicate greater energy stability. In Tables 1 and S3, we can see that the 2‐nitroimidazoles that contain mercaptopyridine, benzyl, and phthalazine functional groups in the R9 substituent, respectively (2NS5‐2, 2NS5‐3, and 2NS3‐7—as well as BNZ, M35m)—have negative values of the SESE descriptor. Therefore, they are the most energetically stable molecules, which agrees with the results obtained regarding chemical hardness and ionization potential. However, the stability decreases when the nitro group is in position 5, as it can be observed in the M7D, M3, 5NS2‐6, and 5NS1‐1 molecules, which indicates that the 5‐nitroimidazole molecules are less stable than the 2‐nitroimidazoles.

Since toxicity is a critical parameter in drug design, the descriptors electrophilicity index (ω) and aromaticity index (Δ) were determined. In the last decades, these descriptors have been related to mutagenicity or carcinogenicity of organic molecules, respectively. In this context, Roy et al., , have related high values of electrophilicity index to biological toxicity due to organic molecules with high values could form covalent bonds with DNA, which might cause mutations. In this respect, it was observed in Tables 1 and S3 that nonnitrated molecules S1‐1, S1‐6, S2‐1, S2‐6, S3‐1, S3‐6, and S4‐1, have the largest electrophilicity index values, therefore, they could be considered as the molecules with the highest probability of causing mutagenicity; on the other hand, nitrated structures and nonnitrated molecules S1‐8, S2‐8, S3‐8, and S4‐8, possess the lowest values, and they could be considered the least mutagenic structures in this study. In addition, the electrophilicity index (ω) values showed that the presence of the nitro group decreases the electrophilicity index of nitroimidazole derivatives, that is, BNZ, reference molecules and nitrated molecules have lower electrophilicity index values than their respective nonnitrated analogs, therefore, nitrated structures are less susceptible to accepting electronic charge. This behavior may be related to the global reactivity of the molecules, which is consistent with the results of chemical hardness and ionization potential.

The aromaticity index (Δ) is a global descriptor that indirectly predicts the carcinogenicity of a compound when it has aromatic rings. , According to the Barone indexes, a molecule is considered carcinogenic if it presents values above 0.36 eV. Tables 1 and S3 show that the molecules S2‐1, S3‐1, and S4‐1 have highest aromaticity values, therefore, these molecules could be considered carcinogenic. By contrast, it is observed that the molecules 5NS2‐6, S1‐6, 2NS5‐2, S2‐8, S3‐8, M35m, 2NS1‐6, and 2NS3‐6 possess the lowest values, which indicates that they can be considered as the molecules with the least probability of causing carcinogenicity.

On the other hand, the Tables 2 and S4 show the atomic charge values (fitted to ESP) of nonnitrated molecules and nitroimidazole derivatives for the following atoms: N1, C2, N3, C4, C5, and the substituent (R9), as well as NO2 groups, in order to estimate which atomic or group charges could be related to the active sites. It is observed that the negative charges were located in N3 atom and NO2 group, whereas the positive charges were located in C2, C4 atoms, and R9 substituent groups, except for S3‐6 and 5NS1‐1. Similarly, the results suggest that N3 atoms have larger negative values of the atomic charges than N1 atoms in all molecules, which agrees with the PAs and the amphoteric properties of imidazole derivatives reported in previous studies.

TABLE 2

Atomic charges values for imidazole and nitroimidazole derivatives at the M06‐2X/6‐311++G(2df,2p)//6‐311+G(d,p) levels

Molecule	N1	C2	N3	C4	C5	O7	O8	R9	NO2
BNZ*	0.314	0.201	−0.595	0.221	−0.347	−0.511	−0.528	0.071	−0.212
M3**	0.257	0.444	−0.683	0.269	−0.357	−0.558	−0.564	0.050	−0.169
M35m	−0.228	0.432	−0.604	0.059	0.143	−0.489	−0.453	0.308	−0.309
M2**	0.075	0.706	−0.649	0.207	−0.316	−0.569	−0.581	0.100	−0.191
M7D	−0.021	0.545	−0.683	0.244	−0.243	−0.547	−0.556	0.167	−0.183
S1‐1	0.291	0.167	−0.685	0.220	−0.471	–	–	0.028	–
S1‐6	0.300	0.147	−0.680	0.199	−0.419	–	–	0.030	–
S1‐8	0.172	0.211	−0.694	0.212	−0.388	–	–	0.098	–
S2‐1	−0.039	0.525	−0.738	0.171	−0.395	–	–	0.200	–
S2‐6	0.229	0.408	−0.727	0.152	−0.449	–	–	0.047	–
S2‐8	0.121	0.446	−0.732	0.152	−0.398	–	–	0.105	–
S3‐1	−0.074	0.235	−0.689	0.060	0.060	–	–	0.182	–
S3‐6	0.250	0.135	−0.678	0.073	−0.106	–	–	−0.002	–
S3‐8	0.082	0.189	−0.677	0.046	0.002	–	–	0.101	–
S4‐1	0.027	0.682	−0.677	0.052	−0.378	–	–	0.116	–
S4‐8	−0.068	0.743	−0.695	0.075	−0.365	–	–	0.172	–
S5‐1	0.129	0.454	−0.724	0.191	−0.402	–	–	0.058	–
2NS5‐2	0.280	0.242	−0.575	0.201	−0.273	−0.459	−0.419	0.062	−0.244
2NS1‐6	0.367	0.178	−0.585	0.222	−0.355	−0.486	−0.477	0.022	−0.194
2NS5‐3	0.026	0.339	−0.579	0.078	0.059	−0.467	−0.426	0.132	−0.279
2NS3‐6	0.226	0.177	−0.567	0.080	−0.002	−0.504	−0.483	0.059	−0.205
2NS3‐7	0.158	0.183	−0.560	0.062	0.030	−0.525	−0.531	0.122	−0.230
5NS2‐6	0.167	0.461	−0.681	0.236	−0.271	−0.487	−0.546	0.068	−0.186
5NS1‐1	0.378	0.228	−0.630	0.259	−0.328	−0.540	−0.509	−0.023	−0.148
5NS5‐4	−0.045	0.778	−0.683	0.209	−0.209	−0.524	−0.453	0.177	−0.236

Note: Atomic charges of atoms and groups are in a.u.

Values obtained from Reference [16].

Values obtained from Reference [15]

In addition, it is observed that the atomic charges of the N1, N3, C5 atoms, and the R9‐substituent are modified due to the presence of the nitro group and of its substituents R1= –CH3 and R2= –OCH3, as shown for the BNZ, S1‐1, M35m, S3‐1, M3, S2‐1, M2, and S4‐1 molecules (see Tables 2 and S4). It is interesting to note that N3 atoms for no‐nitrated analogs have larger negative values of the atomic charges than nitroimidazole derivatives; these results show that the electron‐withdrawing effect of the NO2 group is stronger than the R9‐substituent effect. In addition, it was observed that for nitroimidazoles, the atomic charges of C2 and C5 atoms are affected by the presence of the nitro group in position two, while the atomic charges of C4 atoms are affected mainly by the presence of the nitro group in position five of the imidazole ring, respectively (see Figures 4 and 5).

The condensed Fukui function allows the analysis of the interaction sites in specific regions of a molecule and the analysis of which atoms are potential reactive sites for electrophilic, nucleophilic, or radical attacks, respectively. Tables 3 and S5 collect the Fukui function values for the following atoms: C2, N3, C4, C5, N6, O7, and O8 for the reference molecules, nitroimidazole derivatives, and their respective nonnitrated analogs.

TABLE 3

Condensed Fukui function values for imidazole and nitroimidazole derivatives at the M06‐2X/6‐311++G(2df,2p)//6‐311+G(d,p) levels

Mol	f −					f +				f °
	C2	N3	C4	C5	R9	N6	O7	O8	R9	C4	C5	O7	O8
BNZ*	–	–	–	–	0.145 C21 0.181 C29	0.164	0.245	0.212	–	0.027	0.034	0.123	0.106
M3**	–	–	–	–	0.148 C16 0.181 C24	0.159	0.229	0.192	–	0.051	0.013	0.115	0.097
M35m	0.132	0.065	0.188	0.145	–	0.163	0.239	0.196	–	0.121	0.099	0.158	0.129
M2**	0.078	0.087	0.115	0.141	–	0.156	0.223	0.191	–	0.113	0.085	0.159	0.129
M7D	0.109	0.085	0.135	0.162	–	0.160	0.230	0.198	–	0.118	0.094	0.164	0.133
S1‐1	0.161	0.090	0.202	0.216	–	–	–	–	0.153 C20 0.152 C24	0.102	0.108	–	–
S1‐6	–	–	–	–	0.123 C22 0.136 C24	–	–	–	0.129 C20 0.133 C24	–	–	–	–
S1‐8	–	–	–	–	0.142 N32 0.167 N33	–	–	–	0.109 C21 0.108 C25	–	–	–	–
S2‐1	0.123	0.087	0.197	0.208	–	–	–	–	0.141 C19 0.145 C23	0.099	0.105	–	–
S2‐6	0.122	0.088	0.196	0.210	–	–	–	–	0.129 C19 0.133 C23	0.099	0.105	–	–
S2‐8	0.123	0.088	0.197	0.210	–	–	–	–	0.109 C20 0.108 C24	0.099	0.105	–	–
S3‐1	0.158	0.089	0.199	0.167	–	–	–	–	0.142 C20 0.145 C21	0.101	0.084	–	–
S3‐6	0.160	0.089	0.199	0.168	–	–	–	–	0.123 C19 0.131 C23	0.100	0.084	–	–
S3‐8	0.158	0.089	0.199	0.167	–	–	–	–	0.109 C20 0.108 C24	0.100	0.084	–	–
S4‐1	0.094	0.091	0.183	0.197	–	–	–	–	0.147 C19 0.148 C23	0.092	0.099	–	–
S4‐8	0.094	0.091	0.183	0.198	–	–	–	–	0.109 C20 0.108 C24	0.092	0.099	–	–
S5‐1	0.116	0.087	0.194	0.202	–	–	–	–	0.144 C17 0.141 N23	0.100	0.102	–	–
2NS5‐2	–	–	–	–	0.169 S11 0.124 C17	0.162	0.244	0.206	–	0.029	0.031	0.129	0.107
2NS1‐6	–	–	–	–	0.118 C21 0.133 C23	0.155	0.225	0.193	–	0.027	0.036	0.116	0.098
2NS5‐3	0.125	–	0.187	0.141	–	0.160	0.240	0.201	–	0.069	0.061	0.159	0.131
2NS3‐6	–	–	–	–	0.122 C20 0.129 C22	0.153	0.222	0.190	–	0.028	0.028	0.115	0.096
2NS3‐7	0.058	0.034	0.083	0.067	0.078 C19 0.078 C23	0.156	0.226	0.191	–	–	–	0.127	0.110
5NS2‐6	–	–	–	–	0.131 C18 0.145 C22	0.159	0.227	0.184	–	0.051	0.012	0.115	0.092
5NS1‐1	–	–	–	–	0.143 C16 0.179 C23	0.156	0.221	0.197	–	0.049	0.011	0.114	0.101
5NS5‐4	0.070	0.083	0.107	0.134	0.099 O10	0.154	0.225	0.196	–	0.106	0.077	0.161	0.132

Values obtained from Reference [16].

Values obtained from Reference [15].

From these results, it can be noted that for nonnitrated imidazole derivatives, the negative function is found mainly on the C2, N3 C4, and C5 atoms; likewise, the Fukui functions, for nitroimidazole derivatives, as well as for the S1‐6 and S1‐8 molecules are mainly located on the R9‐substituent. Therefore, these results indicate that these sites can undergo possible electrophilic attacks. The positive function for nitroimidazole derivatives is principally located on the N6, O7, and O8 atoms, and for nonnitrated imidazole derivatives, is located on the R9‐substituent; these results indicate that these sites could suffer a nucleophilic attack. The radical function values for nonnitrated molecules are located on the C4 and C5 atoms, and, for nitroimidazole derivatives on the O7 and O8 atoms, suggesting possible radical attacks. It is interesting to mention that, in some cases, the regions susceptible to nucleophilic, electrophilic, and radical attacks change their position when the nitro group is not present in the structure, that is, the inductive effect of the nitro group in position two or five modifies the reactivity sites in the nitroimidazole molecules. These results agree with the reported in previous studies.

The frontier molecular orbitals allow to identify the most important interaction sites susceptible to electrophilic and nucleophilic attacks of a molecule of pharmacological interest. Figure 6 shows the molecular orbitals (HOMO and LUMO) of the reference molecules, nitroimidazole derivatives and their nonnitrated analogs, respectively. Hence, the HOMO orbitals for BNZ and M3 are located on the R9 substituent, while for M7D, M2, and M35m, they are mainly located on the imidazole ring, which indicates that these regions are susceptible to electrophilic attacks. On the other hand, the LUMO orbitals are located on the imidazole ring and on the nitro group for BNZ, M3, M35m, M7D, and M2, which indicates that these regions could undergo nucleophilic attacks, respectively. On the contrary, in nonnitrated molecules, the HOMO orbitals are located on the imidazole ring, while the LUMO orbitals are located on the R9 substituent. Interestingly, in all molecules, the HOMO and LUMO molecular orbitals are located in the same areas as the Fukui functions and the atomic charges previously analyzed, confirming the most important interaction sites for these structures.

FIGURE 6

Frontier molecular orbitals (HOMO and LUMO) for imidazole and nitroimidazole derivatives at the M06‐2X/6‐311+G(d,p) level of theory

Template selection and structural alignment

In absence of the crystal structure of the enzyme NADH‐FUM reductase of T. cruzi, , , we carried out homology modeling to obtain its three‐dimensional structure and to use this model in subsequent molecular docking studies. Through a structural quality analysis, the templates that met the criteria required for such modeling were selected (percentage of coverage >30%, percentage of identity > 30%, E‐value < 10‐4, resolution < 3.5, score > 200, R‐value < 0.2, and R‐free > R‐value).

Table 4 shows the four template structures with the highest sequence identity that were identified: PDB ID: 5GLG, chain A with 37.33% identity, crystallized ligand: succinate, coenzyme: FAD; PDB ID: 1KSU, chain A with 35.39% identity, crystallized ligand: FUM, coenzyme: FAD; PDB ID: 1QJD, chain A with 35.19% identity, crystallized ligand: malate, coenzyme: FAD; PDB ID: 1QO8, chain A with 34.79% identity, no crystallized ligand, coenzyme: FAD.

TABLE 4

Structural quality criteria

ID	% Identity	R‐value	R‐free	% Coverage	Resolution	Score	E‐value	Coenzyme	Ligand
5GLG	37.3	0.18	0.21	40	1.80*	272	5.00E‐80	FAD	Succinate
1KSU	35.4	0.17	0.24	39	2.00**	239	1.00E‐67	FAD	Fumarate
1QJD	35.2	–	0.27	39	1.80***	238	5.00E‐67	FAD	Malate
1QO8	34.8	0.23	0.28	41	2.15****	240	5.00E‐68	FAD	–

Method X‐ray diffraction (Values obtained from Reference [51]).

Method X‐ray diffraction (Values obtained from Reference [52]).

Method X‐ray diffraction (Values obtained from Reference [53]).

Method X‐ray diffraction (Values obtained from Reference [54]).

Although available templates in the Protein Data Bank database , , , showed a coverage percentage of ≤41%, they include the segment spanning from residues 315 to 950, which correspond to the catalytic site of FUM and the site of interaction for the formation of the Coenzyme‐FRD complex, as reported in literature. , , , , Figure 7 shows the multiple alignment where the residues that integrate the catalytic site of the protein with FAD are highlighted in green, whereas the residues that form the catalytic site of the protein with FUM are highlighted in yellow.

FIGURE 7

Multiple sequence alignment for template and problem sequence. The residues that integrate the catalytic site of the FAD are in green box, and the residues that integrate the catalytic site of the fumarate are in yellow box

The 10,000 models generated based on each of the four proposed templates are shown Figure 8. From these results, the models with the best values of RMSD and DOPE score were selected. Based on 1QJD template, the model 2735 was selected, with values of RMSD = 3.87 Å and DOPE = −56,774.23. For the 1KSU template, the model 5299 was selected, with values of RMSD = 3.87 Å and DOPE = −57,022.21. For the 5GLG template, the model 4560 was selected, with values of RMSD = 2.01 Å and DOPE = −49,050.61. For template 1QO8, the model 9109 was selected, with values of RMSD = 8.34 Å and DOPE = −52,790.86.

FIGURE 8

Models generated for: (A) 1QJD (purple), (B) 1KSU (blue), (C) 5GLG (red), (D) 1QO8 (green). The values of the root mean square deviation (RMSD) were calculated with ProFit, , and the DOPE statistical score (discrete optimized protein energy) was calculated by MODELLER 9.21

We visually evaluated and compared the selected models with their respective templates using the PyMOL software, as shown in Figure 9. The results illustrated that the modeled structures are highly similar to their respective templates, especially in the regions of the catalytic site, therefore, Ramachandran plots were constructed to evaluate the stereochemical quality of models, as shown in Figure 10.

FIGURE 9

Selected models (magenta) and compared with their respective templates (blue): (A) 1QJD vs the model 2735, (B) 1KSU vs the model 5299, (C) 5GLG vs the model 4560, (D) 1QO8 vs the model 9109. The models were generated by MODELLER 9.21

FIGURE 10

Ramachandran plots for the selected models: (A) model 2735 (1QJD), (B) model 5299 (1KSU), (C) model 4560 (5GLG), (D) model 9109 (1QO8)

The Ramachandran plots showed that the amino acids of each selected model were located in allowed sites. In the case of the model 2735 (1QJD), a total percentage of 98.7% was obtained; for the model 5299 (1KSU), a total percentage of 98.7% was obtained; for the model 4560 (5GLG), a total percentage of 98.3% was obtained; and for the model 9109 (1QO8), a total percentage of 98.5% was obtained. According to our results shown in the homology modeling study, was confirmed that the modeled structures have the quality required to continue with the molecular docking study.

Docking studies of coenzyme‐protein complex

We performed molecular docking between the NADH and the model 5299 (1KSU), the NADH and the model 2735 (1QJD), the NADH and the model 4560 (5GLG), as well as the NADH and the model 9109 (1QO8), to form the respective NADH‐FRD complexes. We also compared the spatial location of the coenzyme, taking as reference the position of the crystallized FAD in the template structures, as shown in Figure 11.

FIGURE 11

Molecular docking for: (A) The crystallized FAD‐1QJD (green) overlaid with NADH coupled to model 2735 (blue), (B) the crystallized FAD‐1KSU (green) overlaid with NADH coupled to model 5299 (red), (C) the crystallized FAD‐5GLG (green) overlaid with NADH coupled to model 4560 (yellow), and (D) the crystallized FAD‐1QO8 (green) overlaid with NADH coupled to model 9109 (magenta)

The results shown in Figure 11 indicated that, the NADH‐5299 (1KSU), is the only one out of the four analyzed complexes that conserved the interaction space of the FAD with its respective template, with a binding energy value of −7.1 kcal/mol. Therefore, we selected the model 5299 to continue with the molecular docking studies and validate the docking protocol (redocking).

The validation of the docking protocol was carried out by means of comparing the spatial location of the crystallized FAD‐1KSU complex, with the FAD coupled to the 1KSU template and with the FAD coupled to model 5299. Similarly, the spatial locations of NADH coenzyme coupled to the 1KSU template, with NADH coupled to model 5299, were compared. All the comparisons were carried out through the RMSD values, as shown in Figure 12.

FIGURE 12

Spatial location of the crystallized FAD (green) overlaid with FAD coupled to model 5299 (blue), the FAD coupled to 1KSU template (magenta), the NADH coupled to 1KSU template (yellow), and the NADH coupled to model 5299 (red)

Figure 12 shows that the binding positions presented by the FAD coupled to 1KSU template, versus the crystallized FAD‐1KSU complex, have a minimal difference in their spatial location, with RMSD = 0.136 Å. In addition, it shows that the location of the FAD coupled to model 5299, agrees with the position of the crystallized FAD at 1KSU, with a difference of RMSD = 1.144 Å. Likewise, the coupling positions between the complex NADH‐5299 and NADH‐1KSU are very similar, with a difference of RMSD = 1.719 Å. It is important to note that the spatial locations of all the complexes generated in the validation of the coupling agree with the position of the crystallized FAD at the 1KSU template; therefore, these results allowed us to confirm that the NADH possesses a very similar interaction site to that of the crystallized FAD.

To analyze hydrogen bonds and hydrophobic interactions, we compared the FAD‐1KSU complex with the NADH‐5299 complex. The interactions were visualized by the LigPlot and PyMOL software, as shown in Figure 13.

FIGURE 13

(A) Interactions of the FAD‐1KSU complex; (B) interactions of the NADH‐5299 complex. The green dotted lines illustrate the hydrogen bonds interactions, and red dotted lines display the hydrophobic interactions

Figure 13 shows the intermolecular interactions that correspond to the hydrogen bonds for FAD‐1KSU complex are located in the 1KSU residues: Gly278 (with N5 atom from 6‐aminopurine and N atom from the amine substituent of the 6‐aminopurine of the FAD), Asp344 (with N atom from the 6‐aminopurine of the FAD), Glu534 (with O atom from the phosphate of the flavin nucleotide of the FAD), Gly136 (with O atom from the phosphate of the flavin nucleotide of the FAD), Ala137 (with O atom from the phosphate of the flavin nucleotide of the FAD), Glu156 (with O atom from the hydroxyl ribose of the adenine nucleotide of the FAD), Ala549 (with O atom from the flavin group of the FAD), Ile550 (with O atom from the flavin group of the FAD), Gly171 (with O and N atoms from the flavin group of the FAD), Gly170 (with O atom from the flavin group of the FAD), Ala169 (with O and N atoms from the flavin group of the FAD), Ala165 (with O atom from the phosphate of the adenine nucleotide of the FAD), Asn164 (with O atom from the phosphate of the adenine nucleotide of the FAD), and Gln338 (with O atom from the hydroxyl ribose of the FAD). Similarly, it was observed that its hydrophobic interactions are located in the following residues: Val132, Arg277, Gly133, Thr313, Gly314, Gly163, Gly533, Ala312, Ile553, Gly135, Tyr505, Asn548, Arg544, Gly547, Val253, Ala168, Leu167, His504, Asn337, Met375, Thr336, Glu156, Lys157, Glu158, Thr276.

Moreover, the intermolecular interactions that correspond to the hydrogen bonds of the NADH‐5299 complex are located in the residues: Asn405 (with N atom from the nicotinamide of the NADH), His552 (with O atom from the ribose of the nicotinamide nucleotide of the NADH), Cys609 (with O atom from the hydroxyl ribose of the adenine nucleotide of the NADH), Cys139 (with N atom from the 6‐aminopurine of the NADH), Ala337 (with N8 atom from the 6‐aminopurine of the NADH), Ala142 (with O atom from the phosphate of the nicotinamide nucleotide of the NADH), Leu141 (with O atom from the phosphate of the nicotinamide nucleotide of the NADH), Thr364 (with O atom from the hydroxyl ribose of the nicotinamide nucleotide of the NADH), while their hydrophobic interactions are located in the residues: Arg600, Ala173, Pro412, Ser605, Pro261, Leu606, Ala142, Thr260, Gly140, Thr338, Trp368, Gly138, Gly339, Gly168, Gly167, Ser170, Glu590, Asn365, Ser170, Glu590, Asn365.

Validation of the water docking protocol

The water docking protocol was validated by assessing whether programs could reproduce conserved water molecules' presence or absence in an experimental 3D protein crystal structure. High‐resolution complex crystal structure 1KSU, which is homologous to enzyme NADH‐FUM reductase, served as input for the validation analysis. This protein is bound to its coenzyme, FAD, and to its native ligand, FUM. Importantly, crystallographic data show no evidence of the presence of conserved water molecules in the vicinity of the ligand (FUM)‐binding pocket of 1KSU, as shown in Figure 14.

FIGURE 14

Crystal structure of 1KSU reveals no evidence of the presence of conserved water molecules within the ligand (fumarate, FUM)‐binding pocket

When water was docked into FAD‐1KSU complex, in the absence of ligand (FUM), six water molecules were identified as potentially conserved in the environs of the ligand‐binding pocket (Figure 15A, left panel). However, only two of those water molecules were found to interact with the complex through hydrogen bonds with FAD. Hydrogen bonds show long lengths, ranging from 2.88 Å and up to 3.06 Å, suggesting that the interactions may be considerably weak, , and thus, significantly reducing their influence in the stabilization of the FAD‐1KSU complex.

FIGURE 15

Docking of water molecules to (A) FAD‐1KSU complex and to (B) FAD‐1KSU‐fumarate (FUM) complex (red spheres represent water molecules)

Interestingly, in the explicit presence of the ligand (FUM), the docking of water molecules into FAD‐1KSU‐FUM complex revealed no conserved water molecules within the ligand‐binding pocket (Figure 15B, left panel). This may confirm that the interactions of conserved water molecules, identified previously, are indeed weak and that such molecules are completely displaced by the presence of the ligand, FUM, which is strongly stabilized through a series of hydrogen bonds and hydrophobic interactions with the protein (Figure 15B, right panel). These results agree with the lack of conserved water molecules in the crystallographic structure of 1KSU, therefore, these finding indicate that our protocol may be suitable to predict conserved (nondisplaceable) waters in similar ligand–protein complexes.

Docking of water molecules into the homology model of NADH‐fumarate reductase

We used the validated protocol to investigate the potential presence of conserved water molecules within the ligand–binding pocket of the enzyme NADH‐FUM reductase. To conduct this study, in the absence of a 3D crystal structure of the enzyme to T. cruzi, we employed a homology model (model 5299, see section Model Building), which was built based on the structure of 1KSU. This protein model is bound to NADH (coenzyme) and to FUM (ligand).

Without the explicit presence of FUM, the docking of water molecules into the homology model found six conserved water molecules (Figure 16, left panel). However, these water molecules appear to be somewhat displaced toward the vicinity of the NADH‐binding pocket, where four of such molecules interact with the coenzyme by forming hydrogen bonds, Figure 16(A) (right panel). This suggests that the influence of these water molecules may be insignificant within the FUM‐binding site.

FIGURE 16

Docking of water molecules to (A) NADH‐5299 complex and (B) NADH‐5299‐fumarate (FUM) complex (red spheres represent water molecules)

In the explicit presence of FUM (ligand), water docking into the NADH‐5299‐FUM complex identified eight potential nondisplaceable water molecules (Figure 16B, left panel). But only four of such molecules form hydrogen bonds with FUM (Figure 16B, right panel). However, the interaction distances between FUM and water molecules are significantly long, with bond lengths of 2.86, 2.91, 3.00, and 3.00 Å, which, according to the literature, result in weak interactions. ,

Thus, our in silico docking assessment suggest the absence of conserved water molecules in the vicinity of the FUM‐binding site of the enzyme NADH‐FUM reductase. These findings are consistent with experimental crystallographic data. High‐resolution crystal structure of 1KSU reveals no conserved water molecules, and crucially, it is an evolutionary homolog to the enzyme NADH‐FUM reductase. Hence, this reductase, similar to 1KSU, should also lack conserved water molecules, which our findings precisely suggest. Therefore, subsequent docking studies involving this enzyme may be conducted without incorporating explicit water molecules.

Docking studies of coenzyme–protein–ligand complex

Once the NADH‐FRD (NADH‐5299) complex was obtained and validated, the molecular docking of the complex with FUM, reference molecules and selected ligands were carried out. In addition, to corroborate that the interactions are found in the internal part of the protein and not on the surface, a visual analysis was carried out of the spatial location of the designed ligands, as shown in Figure 17.

FIGURE 17

Molecular docking of the NADH‐5299 complex with fumarate, selected ligands and reference molecules that conserve the catalytic site of the fumarate

The structures shown in Figure 17 correspond to the molecules S1‐1, S1‐6, S1‐8, S2‐8, S3‐6, S3‐8, S4‐8, 5NS1‐1, and 2NS1‐6 which were found in the internalized region of the protein and conserved the catalytic site of FUM with the NADH‐5299 complex. By contrast, the structures BNZ, M3, M35m, M2, M7D, S2‐6, S4‐6, 5NS2‐6, 2NS3‐6, 2NS3‐7, 2NS5‐2, 2NS5‐3, and 5NS5‐4 were located on the superficial region of the protein, which hinders their interaction with the catalytic site of FUM.

In order to analyze the binding energies and intermolecular interactions among the NADH‐5299 complex and the selected ligands, in Table 5, we formatted against an italic the ligands that conserved the FUM catalytic site; so, it can be observed that the main intermolecular interactions of FUM for the hydrogen bonds are located on NADH, Lys172, Tyr480, Arg89; while for the hydrophobic interactions, the residues are located on Asp258, Phe479, Thr174, Met63. It is important to mention that the residues: Lys172 and Arg89 for the hydrogen bonds and the residues: Phe479, Thr174, Met63 for the hydrophobic interactions were localized in almost all coupled ligands.

TABLE 5

Ligand interactions with the NADH‐5299 complex

Ligands	Binding energy (kcal/mol)	Hydrophobic interactions	Hydrogen bonds
FUM	−4.7	Asp258, Phe479, Thr174, Met63	NAD, Lys172, Tyr480, Arg89
BNZ	−6.2	Ala62, Arg89, Ile57, Thr483, Met63, Thr174, Tyr480, Tyr64, Phe479, NAD	Lys172
M3	−6.4	Ala62, Leu85, Ile57, Arg89, Asp258, Lys172, Thr174, Phe479, Lys58, Met63, Tyr64, Ser61, NAD
M35m	−6.0	Ile57, Ala62, Leu85, Phe479, Met63, Tyr64, Asp258, Thr174, Arg89, Lys172, NAD
M2	−6.7	Arg89, Lys172, Asp258, Thr174, Phe479, Met63, Tyr64, Ala62, Lys58, Ser61, Ile57, Leu85
M7D	−5.7	Ser56, Lys58, Met63, Leu85, Ala62, Ser61, Thr483, Phe479	Asn482
S1‐1	−6.0	Tyr64, Ile57, Thr483, Leu85, Phe479, Met63, Thr174, NAD	Arg89, Lys172
S1‐6	−6.0	Tyr480, Thr174, Phe479, Met63, Ala62, Ile57, Leu85	Arg89, Lys172
S1‐8	−7.2	Ile57, Ala62, Met63, Tyr480, Thr174, Phe479, Leu85	Arg89, Lys172
S2‐6	−6.7	Leu85, Ala62, Met63, Ser56, Phe479, Asn482, Ile57, Thr483, Lys58	Arg89, Ser61
S2‐8	−7.5	Leu85, Ile57, Tyr64, Tyr480, Phe479, Thr174, Met63	Arg89, Lys172
S3‐6	−6.1	Met63, Tyr64, Ile57, Leu85, Arg89, Phe479, Thr174, Tyr480, NAD	Lys172
S3‐8	−8.0	Ala62, Met63, Tyr480, Thr174, Phe479, Tyr64, Leu85, Ile57, NAD	Arg89, Lys172, Ser61
S4‐6	−6.4	Thr483, Asn482, Met63, Phe479, Ser56, Arg89, Leu85, Ile57, Lys58, Ala62	Ser61
S4‐8	−6.8	Leu85, Ile57, Thr174, Tyr480, Phe479, Met63	Arg89, Lys172
2NS5‐2	−6.1	Ser61, Thr483, Ser56, Phe479, Met63, Tyr64, Ala62, Leu85, Ile57	Arg89, Lys58
2NS1‐6	−6.8	Thr174, Met63, Arg89, Ile57, Leu85, Tyr64, Phe479	Tyr480, Lys172, NAD
2NS5‐3	−5.8	Phe479, Thr483, Ser56, Met63, Leu85, Ala62, Ser61	Asn482
2NS3‐6	−7.0	Phe479, Lys58, Ile57, Arg89, Ala62, Leu85, Met63, Thr483	Asn482, Ser56
2NS3‐7	−7.3	Phe479, Ser56, Lys58, Ile57, Met63, Ala62, Leu85, Arg89, Thr483	Asn482
5NS2‐6	−6.9	Phe479, Lys58, Met63, Ile57, Ala62, Leu85, Arg89, Thr483	Asn482, Ser56
5NS1‐1	−6.5	Leu85, Tyr64, Phe479, Ser61, Ala62, Met63, Arg89, Lys172, Thr174, NAD
5NS5‐4	−5.8	Ser56, Phe479, Thr483, Ile57, Arg89, Leu85, Ala62	Lys58

The molecular docking studies indicated that the structures S3‐8 (−8.0 kcal/mol), S2‐8 (−7.5 kcal/mol), and S1‐8 (−7.2 kcal/mol) have the lowest values of binding energy.

CONCLUSIONS

A computational study was undertaken, at the DFT‐M06‐2X level in aqueous solution, to characterize and analyze the electronic structure and evaluate the importance of the nitro group and the substituent effects on the chemical reactivity and the physicochemical features of new imidazole and nitroimidazole derivatives. The virtual platforms indicated that the selected molecules were not toxic and had good pharmacokinetic properties.

The results of the quantum chemical descriptors allowed us to identify the structures: S2‐8, S3‐8, 2NS5‐2, 2NS3‐7, 2NS5‐3, and 5NS1‐1, which showed the best physicochemical properties and the lowest toxicity. In addition, the results showed that, the presence of the phthalazine group in the R9 substituent confers more stability to nonnitrated imidazole derivatives.

Using homology modeling and applying a thorough refinement protocol, we constructed a high‐quality 3D structure of the enzyme NADH‐FUM reductase. The homology model was built based on the experimental crystallographic structure from Shewanella frigidimarina (PDB ID:1KSU. It showed the highest stereochemical robustness and structural consistency as demonstrated by the low RMSD, 3.87 Å, when superimposed to the template protein.

This model was then used to investigate the potential presence of conserved water molecules within the ligand–binding pocket of the enzyme. Our in silico docking approach, complemented with experimental crystallographic data, revealed the absence of water molecules in the molecular interactions of NADH‐FUM reductase with ligands. This finding allowed us to perform subsequent ligand–docking studies without incorporating explicit water molecules in the simulations.

The selected molecules and evaluated via DFT studies were tested as ligands in molecular docking. The results showed that the ligands S1‐1, S1‐6, S1‐8, S2‐8, S3‐6, S3‐8, S4‐8, 5NS1‐1, and 2NS1‐6 interacted at the catalytic site of the FUM with similar intermolecular interactions.

Finally, the computational studies showed that the molecules S3‐8 (−8.0 kcal/mol), S2‐8 (−7.5 kcal/mol), and S1‐8 (−7.2 kcal/mol), could be considered as potential inhibitors of the enzyme NADH‐FUM reductase; however, future studies of molecular dynamics simulations and experimental studies of enzymatic inhibition are contemplated.

APPENDIX S1 Supporting information

Click here for additional data file.