Comparisons between Crystallography Data and Theoretical Parameters and the Formation of Intramolecular Hydrogen Bonds: Benznidazole

The conformational preferences of benznidazole were examined through the application of DFT, PCM and QTAIM calculations, whose results were compared with crystallography data. The geometries were fully optimized with minimum potential energy surface by means of the Relaxed Potential Energy Surface Scan (RPESS) at AM1, followed by the B3LYP/6-311++G(d,p) theoretical level. As a result, the s-cis conformation (1C) was shown to be more stable (4.78 kcal¨moí1) than s-trans (1T). The Quantum Theory Atoms in Molecules (QTAIM) was applied in order to characterize the (N–H¨¨¨O=N) and (C–H¨¨¨=N) intramolecular hydrogen bonds. The simulation of solvent effect performed by means of the implicit Polarized Continuum Model (PCM) revealed great results, such as, for instance, that the conformation 1W is more stable (23.17 kcal¨moí1) in comparison to 1C. Our main goal was stressed in the topological description of intramolecular hydrogen bonds in light of the QTAIM approach, as well as in the solvent simulation to accurately obtain an important conformation of benznidazole.


Introduction
In medicinal chemistry, the studies of new compounds with high biological activities and less toxic effects have attracted the attention of many researchers around the world [1].This scientific interest concerns the potential of several compounds in treating serious diseases, ranging from "neglected" diseases such as Schistosoma mansoni [2] to disorders such as cancer [3].In the context of neglected diseases, for which big pharmaceutical companies are not investing in the search for therapeutic alternatives, Chagas disease, also called American trypanosomiasis, is considered the most important parasitic infection of Latin America [4].Chagas disease is manifested in humans and domestic animals living in extreme poverty and rural areas [5].In practice, the infection is caused by the Trypanosoma cruzi protozoan although it is transmitted by the Triatoma infestans insect.In the actual treatment, basically two drugs have been widely used: Benznidazole (Rochagan ® , Basel, Switzerland) and Nifurtimox (Lampit ® , Leverkusen, Germany).Due to the possibility to obtain new compounds to be used in the treatment against Chagas disease, a lot of studies have been conducted in order to improve the potentiality of benznidazole.Some time ago, a theoretical conformational study of benznidazole was performed [6].The results point to the existence of local minima in the potential energy surface, whose structures were compared with some active compounds against T. cruzi: derivatives of tetrahydro-β-carboline.This study was carried out in gas phase by using the AM1 semi empirical Hamiltonian, whose application was conditioned to the modeling of the Potential Energy Surface Scan (PESS).Only then, the optimized geometry in each point of minimum has been determined through the B3LYP functional jointly with the 6-311+G(d) basis set.From an experimental viewpoint, Soares Sobrinho et al. [7] have published a crystallography study about the benznidazole structure based on X-ray diffraction, whose results were very different from those in gas phase, including the formation of intermolecular hydrogen bonds on the crystal packaging.In molecular modeling [8], it is well established that the identification of hydrogen bonds is one of the most important criteria in analysis of biological compounds [9].These studies may be performed by virtual screening [10] or quantum chemical calculations [11,12] in the pursuit of achieving structure-activity relationships.Historically, the formation of intermolecular or intramolecular hydrogen bonds has been discussed by taking into account the van der Waals radii of the electron donor-acceptor.If we consider this statement [5], the distance values for typical hydrogen bonds should be exactly the same as or shorter than 2.6 Å, of course regarding the F, O and N atoms, or even the π cloud as proton-acceptor centers [13][14][15][16].In fact, the characterization of hydrogen bonds is not just dependent on structural parameters, actually, the application of quantum mechanical criteria is feasible in this regard.
The Quantum Theory of Atoms in Molecules (QTAIM) [17] represents a useful tool in studies of molecular stability and strength of chemical bonds [18], but is also widely applied in investigations of hydrogen bonds through the analysis of electronic density and its topological parameters [19].According to some publications [20], the QTAIM approach has been also useful in studies of biological systems [21], mainly to characterize the intramolecular hydrogen bonds of compounds with biological activity [22].In this sense, we are convinced that QTAIM can be decisive in our investigation not only regarding the intramolecular hydrogen bonds [23][24][25][26][27][28][29][30][31], but also to unveil the most stable conformation of benznidazole.Another interesting viewpoint is that the solvent effect may be responsible for drastic changes in several molecular properties, e.g., geometrical deformations; increase in the reaction rate; control of products along the reaction paths, for instance [32,33].If we consider that intramolecular hydrogen bonds can be formed in the minimum structure of benznidazole, the specialized literature informs us that the application of the Polarized Continuum Model (PCM) [34] to evaluate the solvent effect is recommended [35][36][37][38][39].By taking into account the importance of drug action in human organisms, that are predominantly aqueous, this article also presents a theoretical investigation of the solvent effect on the benznidazole structure through the application of PCM protocol, demonstrating the importance of solvation in the conformational study of this bioactive molecule [40].In the context of molecular stabilization, a topological description of intramolecular hydrogen bonds using the QTAIM approach is another important contribution of this work.

Computational Procedure
The procedure of Relaxed Potential Energy Surface Scan (RPESS) was performed by taking particular care with the dihedral angles (θ 1 , θ 2 , θ 3 , and θ 4 ), which are illustrated in Figure 1.These calculations were executed by using the AM1 semiempirical level with angle variations from 0 ˝to 360 ˝, in intervals of 10 ˝.After that, the minimum conformation was fully optimized by the B3LYP/6-311++G(d,p) level of calculation without any symmetry constraint and with all geometries modeled at a minimum of potential energy because no imaginary frequencies were found.As a result, the conformation recognized as 1C was generated.compared with some active compounds against T. cruzi: derivatives of tetrahydro-β-carboline.This study was carried out in gas phase by using the AM1 semi empirical Hamiltonian, whose application was conditioned to the modeling of the Potential Energy Surface Scan (PESS).Only then, the optimized geometry in each point of minimum has been determined through the B3LYP functional jointly with the 6-311+G(d) basis set.From an experimental viewpoint, Soares Sobrinho et al. [7] have published a crystallography study about the benznidazole structure based on X-ray diffraction, whose results were very different from those in gas phase, including the formation of intermolecular hydrogen bonds on the crystal packaging.In molecular modeling [8], it is well established that the identification of hydrogen bonds is one of the most important criteria in analysis of biological compounds [9].These studies may be performed by virtual screening [10] or quantum chemical calculations [11,12] in the pursuit of achieving structure-activity relationships.Historically, the formation of intermolecular or intramolecular hydrogen bonds has been discussed by taking into account the van der Waals radii of the electron donor-acceptor.If we consider this statement [5], the distance values for typical hydrogen bonds should be exactly the same as or shorter than 2.6 Å, of course regarding the F, O and N atoms, or even the π cloud as proton-acceptor centers [13][14][15][16].In fact, the characterization of hydrogen bonds is not just dependent on structural parameters, actually, the application of quantum mechanical criteria is feasible in this regard.
The Quantum Theory of Atoms in Molecules (QTAIM) [17] represents a useful tool in studies of molecular stability and strength of chemical bonds [18], but is also widely applied in investigations of hydrogen bonds through the analysis of electronic density and its topological parameters [19].According to some publications [20], the QTAIM approach has been also useful in studies of biological systems [21], mainly to characterize the intramolecular hydrogen bonds of compounds with biological activity [22].In this sense, we are convinced that QTAIM can be decisive in our investigation not only regarding the intramolecular hydrogen bonds [23][24][25][26][27][28][29][30][31], but also to unveil the most stable conformation of benznidazole.Another interesting viewpoint is that the solvent effect may be responsible for drastic changes in several molecular properties, e.g., geometrical deformations; increase in the reaction rate; control of products along the reaction paths, for instance [32][33].If we consider that intramolecular hydrogen bonds can be formed in the minimum structure of benznidazole, the specialized literature informs us that the application of the Polarized Continuum Model (PCM) [34] to evaluate the solvent effect is recommended [35][36][37][38][39].By taking into account the importance of drug action in human organisms, that are predominantly aqueous, this article also presents a theoretical investigation of the solvent effect on the benznidazole structure through the application of PCM protocol, demonstrating the importance of solvation in the conformational study of this bioactive molecule [40].In the context of molecular stabilization, a topological description of intramolecular hydrogen bonds using the QTAIM approach is another important contribution of this work.

Computational Procedure
The procedure of Relaxed Potential Energy Surface Scan (RPESS) was performed by taking particular care with the dihedral angles (θ1, θ2, θ3, and θ4), which are illustrated in Figure 1.These calculations were executed by using the AM1 semiempirical level with angle variations from 0° to 360°, in intervals of 10°.After that, the minimum conformation was fully optimized by the B3LYP/6-311++G(d,p) level of calculation without any symmetry constraint and with all geometries modeled at a minimum of potential energy because no imaginary frequencies were found.As a result, the conformation recognized as 1C was generated.The crystal structure (1K) was optimized using the same level of calculation presented above.The conformation 1T was obtained from an arbitrary variation in the geometry of 1C, corresponding to a previously reported geometry [6], which remains as a minimum conformation (less stable than 1C) after B3LYP/6-311++G(d,p) optimization.The solvent effect was evaluated on both 1K and 1C through the arguments of the Polarized Continuum Model (PCM) both providing 1W.All these calculations were carried out by GAUSSIAN 03W quantum suite of codes [41].Both the nature of electronic density and the characterization of intramolecular hydrogen bonds were investigated through the QTAIM formalism, whose calculations were performed by the AIM 2000 1.0 software package [42].

Minimum Structure, Bond Lenghts and Vibration Modes
Particularly, the study of dihedral θ 3 (Figure 1) is important due to the high energetic barrier of amide group (HN-C=O) [43].In line with this, it is worthy to assume that two conformations (s-trans and s-cis) must coexist (Figure 2) when the energetic interconversion between them is about 20-22 kcal¨mol ´1.The terminology of s-cis and s-trans (s means to conformation into N-C single bond) was used to describe the relative positioning between the R (RNH) and C=O groups when they are in the same direction (s-cis) or even in an opposite alignment (s-trans), respectively.From RPESS protocol, the lowest energy conformation was selected, and after undergoing an optimization procedure at the B3LYP/6-311++G(d,p) level of theory, the conformation symbolized as 1C was obtained (Figure 3).As such, 1C shows both RN and C=O groups in s-cis position.Moreover, the conformation 1T (Figure 4) was also obtained from the geometry described previously [6] through the optimization at B3LYP/6-311++G(d,p) level of theory.Note that the groups (R-N and C=O) are in s-trans position in 1T.
The crystal structure (1K) was optimized using the same level of calculation presented above.The conformation 1T was obtained from an arbitrary variation in the geometry of 1C, corresponding to a previously reported geometry [6], which remains as a minimum conformation (less stable than 1C) after B3LYP/6-311++G(d,p) optimization.The solvent effect was evaluated on both 1K and 1C through the arguments of the Polarized Continuum Model (PCM) both providing 1W.All these calculations were carried out by GAUSSIAN 03W quantum suite of codes [41].Both the nature of electronic density and the characterization of intramolecular hydrogen bonds were investigated through the QTAIM formalism, whose calculations were performed by the AIM 2000 1.0 software package [42].

Minimum Structure, Bond Lenghts and Vibration Modes
Particularly, the study of dihedral θ3 (Figure 1) is important due to the high energetic barrier of amide group (HN-C=O) [43].In line with this, it is worthy to assume that two conformations (s-trans and s-cis) must coexist (Figure 2) when the energetic interconversion between them is about 20-22 kcal•mol −1 .The terminology of s-cis and s-trans (s means to conformation into N-C single bond) was used to describe the relative positioning between the R (RNH) and C=O groups when they are in the same direction (s-cis) or even in an opposite alignment (s-trans), respectively.From RPESS protocol, the lowest energy conformation was selected, and after undergoing an optimization procedure at the B3LYP/6-311++G(d,p) level of theory, the conformation symbolized as 1C was obtained (Figure 3).As such, 1C shows both RN and C=O groups in s-cis position.Moreover, the conformation 1T (Figure 4) was also obtained from the geometry described previously [6] through the optimization at B3LYP/6-311++G(d,p) level of theory.Note that the groups (R-N and C=O) are in s-trans position in 1T.Although it should be expected that the most stable conformation must be the solvated structure [48,49], geometrically the strength of the intramolecular N=O b •••H c -C hydrogen bond is not sufficient evidence for that.The strongest statement regarding the hydrogen bonds' formation, either intra or intermolecular, actually has its basis in the deformations of the bond lengths that compose it.Although it should be expected that the most stable conformation must be the solvated structure [48,49], geometrically the strength of the intramolecular N=O b ¨¨¨H c -C hydrogen bond is not sufficient evidence for that.The strongest statement regarding the hydrogen bonds' formation, either intra or intermolecular, actually has its basis in the deformations of the bond lengths that compose it.Besides the intramolecular hydrogen bonds, Table 1 also enumerates the bond length results of the donors X-H (N-H a , C-H c , C-H e and C-H d ) and acceptors Y (O b =N) of protons.Although a precise measurement of the bond length variation is not allowed because the monomer configuration of the proton donors cannot be attained, it can be seen that the HBond distances are not in agreement with the bond length variations of X-H as well as of Y, whose values are quite similar.Regarding the characterization of the infrared vibration modes as one of the preconditions to recognize the most stable conformation, the values of stretch frequencies and absorption intensities are invariable, even though a slight difference of 97.3 (1C) and 107.6 km•mol −1 (1K) for IN-H a has been observed.Nevertheless, the new vibration modes could not be unveiled, although by taking into account the long values of the HBond distances beyond 2.000 Ǻ, their stretch frequencies and absorption intensities must be active in an infrared region lower than 50 cm −1 .Besides the intramolecular hydrogen bonds, Table 1 also enumerates the bond length results of the donors X-H (N-H a , C-H c , C-H e and C-H d ) and acceptors Y (O b =N) of protons.Although a precise measurement of the bond length variation is not allowed because the monomer configuration of the proton donors cannot be attained, it can be seen that the HBond distances are not in agreement with the bond length variations of X-H as well as of Y, whose values are quite similar.Regarding the characterization of the infrared vibration modes as one of the preconditions to recognize the most stable conformation, the values of stretch frequencies and absorption intensities are invariable, even though a slight difference of 97.3 (1C) and 107.6 km•mol −1 (1K) for IN-H a has been observed.Nevertheless, the new vibration modes could not be unveiled, although by taking into account the long values of the HBond distances beyond 2.000 Ǻ, their stretch frequencies and absorption intensities must be active in an infrared region lower than 50 cm −1 .Besides the intramolecular hydrogen bonds, Table 1 also enumerates the bond length results of the donors X-H (N-H a , C-H c , C-H e and C-H d ) and acceptors Y (O b =N) of protons.Although a precise measurement of the bond length variation is not allowed because the monomer configuration of the proton donors cannot be attained, it can be seen that the HBond distances are not in agreement with the bond length variations of X-H as well as of Y, whose values are quite similar.Regarding the characterization of the infrared vibration modes as one of the preconditions to recognize the most stable conformation, the values of stretch frequencies and absorption intensities are invariable, even though a slight difference of 97.3 (1C) and 107.6 km¨mol ´1 (1K) for IN-H a has been observed.Nevertheless, the new vibration modes could not be unveiled, although by taking into account the long values of the HBond distances beyond 2.000 Å, their stretch frequencies and absorption intensities must be active in an infrared region lower than 50 cm ´1.

PCM Calculations
The stabilization of 1C, 1T and 1W can be discussed on the basis of thermodynamic properties, whose results organized in Table 2 represent the sum of electronic and zero-point energies (ε 0 + ZPE), sum of electronic and thermal free energies (ε 0 + G corr ) and ∆G, obtained from the difference between (e 0 + G corr ) for conformation and the corresponding value for 1C.The PCM calculations were performed on 1C, although this procedure has been based on the X-Ray crystal structure (1K), which both culminated with the development of the 1W conformation.Regarding Figure 3, there is less tendency to form an intramolecular hydrogen bond, which can be explained by the solvent effect in the stabilization of benznidazole, although the structural interaction strength recommended by the RC-H c ¨¨¨O b =N contact points out that 1K is slightly more strongly bonded than 1C.The value of ´23.17 kcal¨mol ´1 may mean that 1W is the most stable structure.However, this HBond presents a distance of 2.468 Å, which, being one of the longest interactions, which seems to conflict with the stabilization statement presented above.It is important to point out that the solvent effect should reveal electronic energy values that may differ from a hypersurface investigated in the gas phase.Thus, the displayed geometry is consistent with the weakening of the intramolecular interaction favoring the influence of the external environment (aqueous).Really, the solvent effect should be carefully used to explain the additional molecular stabilization, depending on the type of system studied [50], including drug molecules that act in predominantly aqueous environments.

QTAIM Parameters and Intramolecular Hydrogen Bonds
Bond pathways, Bond Critical Points (BCP), and values of ρ followed by ∇ 2 ρ were characterized in terms of the QTAIM topological integrations of the electronic density.We can see, in Figures 5  and 6 the presence of a Bond Critical Point (BCP) between (N-H a ¨¨¨O b =N), indicating the formation of intramolecular hydrogen bonds, corroborating therefore the formulated hypothesis in this work.The same conformation shows a BCP between (C-H c ¨¨¨O b =N) groups, forming a ring of six members [9].The participation of the C-H group in intramolecular hydrogen bonds has been established [51].In our current case, the simultaneous positions of the α-carbonyl and α-nitrogen-hydrogen can contribute to a peculiar acidity of C-H and also lead to the formation of an intramolecular hydrogen bond.The values of ρ and ∇ 2 ρ that characterize these interactions can be visualized in Table 3.By taking into account the virial theorem [17], G and U idealize the kinetic (always positive) and potential (always negative) energy density functions.Then, positive values of ∇ 2 ρ indicate a depletion of electronic potential energy density at BCP in favor of the kinetic energy because it outweighs the potential one, showing a concentration of electronic charge in separated nuclei [51][52][53].

QTAIM Parameters and Intramolecular Hydrogen Bonds
Bond pathways, Bond Critical Points (BCP), and values of ρ followed by ∇ 2 ρ were characterized in terms of the QTAIM topological integrations of the electronic density.We can see, in Figures 5 and 6, the presence of a Bond Critical Point (BCP) between (N-H a •••O b =N), indicating the formation of intramolecular hydrogen bonds, corroborating therefore the formulated hypothesis in this work.The same conformation shows a BCP between (C-H c •••O b =N) groups, forming a ring of six members [9].The participation of the C-H group in intramolecular hydrogen bonds has been established [51].In our current case, the simultaneous positions of the α-carbonyl and α-nitrogen-hydrogen can contribute to a peculiar acidity of C-H and also lead to the formation of an intramolecular hydrogen bond.The values of ρ and ∇ 2 ρ that characterize these interactions can be visualized in Table 3.By taking into account the virial theorem [17], G and U idealize the kinetic (always positive) and potential (always negative) energy density functions.Then, positive values of ∇ 2 ρ indicate a depletion of electronic potential energy density at BCP in favor of the kinetic energy because it outweighs the potential one, showing a concentration of electronic charge in separated nuclei [51][52][53].

QTAIM Parameters and Intramolecular Hydrogen Bonds
Bond pathways, Bond Critical Points (BCP), and values of ρ followed by ∇ 2 ρ were characterized in terms of the QTAIM topological integrations of the electronic density.We can see, in Figures 5 and 6, the presence of a Bond Critical Point (BCP) between (N-H a •••O b =N), indicating the formation of intramolecular hydrogen bonds, corroborating therefore the formulated hypothesis in this work.The same conformation shows a BCP between (C-H c •••O b =N) groups, forming a ring of six members [9].The participation of the C-H group in intramolecular hydrogen bonds has been established [51].In our current case, the simultaneous positions of the α-carbonyl and α-nitrogen-hydrogen can contribute to a peculiar acidity of C-H and also lead to the formation of an intramolecular hydrogen bond.The values of ρ and ∇ 2 ρ that characterize these interactions can be visualized in Table 3.By taking into account the virial theorem [17], G and U idealize the kinetic (always positive) and potential (always negative) energy density functions.Then, positive values of ∇ 2 ρ indicate a depletion of electronic potential energy density at BCP in favor of the kinetic energy because it outweighs the potential one, showing a concentration of electronic charge in separated nuclei [51][52][53].These values also indicate intramolecular hydrogen bonds with a significant covalent character, as argued by Gilli et al. [54] and Siskos et al. [55].All these considerations can also be visualized in terms of the Resonance Assisted Hydrogen Bond phenomenon (RAHB) [56,57], in which the π electrons of the pirazole system are in resonance with the NO 2 group, increasing the electron density at the BCP and decreasing the bond lengths, that are considerably smaller than sum of van der Walls radii.On the other hand, it is important to be highlighted that the values of ∇ 2 ρ, ρ presented here are much smaller than the ones obtained at BCP for common covalent interactions.
Regarding the interaction strength, the ´G/U ratio is a qualitative manner to express the appearing of covalent effect [32].According to values of 0.733 and 0.738 gathered in Table 3, the C-H c ¨¨¨O b =N intramolecular hydrogen bonds behave with a trend of covalence in 1K and 1W, although the electronic density values are not the largest, and the bond lengths are not the shortedones, as demonstrated in Figure 7.This can be stated once these results are smaller than 1.0 [32,58].In regards to the remaining hydrogen bonds, all of them are non-covalent.However, besides the hydrogen bonds, 1T also reveals the formation of a dihydrogen bond [59][60][61], namely C-H d ¨¨¨H f -C, wherein by means of the values of ρ and ∇ 2 ρ as well as ´G/U, this interaction is also non-covalent.

*
All values are given in electronic units (e.u.).
These values also indicate intramolecular hydrogen bonds with a significant covalent character, as argued by Gilli et al. [54] and Siskos et al. [55].All these considerations can also be visualized in terms of the Resonance Assisted Hydrogen Bond phenomenon (RAHB) [56,57], in which the π electrons of the pirazole system are in resonance with the NO2 group, increasing the electron density at the BCP and decreasing the bond lengths, that are considerably smaller than sum of van der Walls radii.On the other hand, it is important to be highlighted that the values of ∇ 2 ρ, ρ presented here are much smaller than the ones obtained at BCP for common covalent interactions.
Regarding the interaction strength, the −G/U ratio is a qualitative manner to express the appearing of covalent effect [32].According to values of 0.733 and 0.738 gathered in Table 3, the C-H c •••O b =N intramolecular hydrogen bonds behave with a trend of covalence in 1K and 1W, although the electronic density values are not the largest, and the bond lengths are not the shortedones, as demonstrated in Figure 7.This can be stated once these results are smaller than 1.0 [32,58].In regards to the remaining hydrogen bonds, all of them are non-covalent.However, besides the hydrogen bonds, 1T also reveals the formation of a dihydrogen bond [59][60][61], namely C-H d •••H f -C, wherein by means of the values of ρ and ∇ 2 ρ as well as −G/U, this interaction is also non-covalent.

Conclusions
In this work, we observed that intramolecular hydrogen bonds can control the conformations of benznidazole.The structural parameters embodied as bond lengths and intermolecular distances revealed unsystematic tendencies in comparison with the stabilization criteria.Also, the new vibration modes of the intramolecular hydrogen bonds could not be unveiled, and in addition, the difficulties in identifying the frequency shifts of the proton donors prevent us from correlating this parameter with the stabilization criteria.The QTAIM calculations were able to characterize classical and non-classical intramolecular hydrogen bonds by means of the values of electronic density and Laplacian.Also, the appearance of a partial covalent character was observed in the C-H c ¨¨¨O b =N links of 1K and 1W.Approaches considering solvent or gas phase effect can lead to deviations of the crystal structure, but they are all s-cys and not s-trans.This must be taken into account in studies of the structure-activity relationships of benznidazole.

Figure 1 .
Figure 1.Structure of benznidazole and dihedral angles examined in RPESS procedure.

Figure 1 .
Figure 1.Structure of benznidazole and dihedral angles examined in RPESS procedure.

Figure 7 .
Figure 7. Relationship between the distance values of intramolecular hydrogen bonds and QTAIM electronic density amounts.

Figure 7 .
Figure 7. Relationship between the distance values of intramolecular hydrogen bonds and QTAIM electronic density amounts.

Table 1 .
Values of HBond distances and bond lengths of donors and acceptors of protons in the 1C, 1T, 1K and 1W conformations of benznidazole.
* Values of R (Intramolecular hydrogen bonds) and r (bond lengths) are given in Å; Values of υ (stretch frequencies) and I (absorption intensity) are given in cm ´1 and km¨mol
a 1C is the reference minimum.
* All values are given in electronic units (e.u.).